Gab1Bdicussão, R - rextester

Language: Layout:

# Ajuste de curvas por mínimos quadrados.

# Dados os vetores x e y, supomos que y_i=f(x_i)+erro_i, como erro_i um erro "pequeno".
x=seq(1930,2010,by=10);x 
x=(x-x[1])/10 
n=length(x)
y=c(0.63, 0.76 ,0.92 ,1.1, 1.2, 1.3, 1.5, 1.7 ,1.9)

plot(x,y,col="red")

# Seguem os comando para as definições de funções necessárias.

g<-function(u){    # função g(u)
a=u[1];b=u[2];c=u[3]
p=0
for (i in 1:n){
p=p+( y[i]- a/(b+c*exp(-a*x[i]) ) )^2}
p/2
}

# Teste do código para g(u).
g(c(1,2,3))

gradg<-function(u){   # função gradiente de g(u).
a=u[1];b=u[2];c=u[3]
hh=0.0001
p=0*u
v=u; v[1]=u[1]+hh; w=u; w[1]=u[1]-hh; p[1]=(g(v)-g(w))/(2*hh) #dg/da

v=u; v[2]=u[2]+hh; w=u; w[2]=u[2]-hh; p[2]=(g(v)-g(w))/(2*hh) #dg/db

v=u; v[3]=u[3]+hh; w=u; w[3]=u[3]-hh; p[3]=(g(v)-g(w))/(2*hh) #dg/dc

}

# Teste da função gradg(u)
gradg(c(1,2,3))

#---- Início da rotin para o método de Euler.

ZeroEuler<-function(u0,t,m){
w=u0
h=t/m
for (i in 1:m){w=w-h*gradg(w)}
w
}

# Teste de aproximação de solução do problema de minimização. 
u0=c(0.5,0.5,0.5);print("Condição inicial");u0; 
t=10; m=100000
print("Teste de g(u0)");g(u0)                       # Teste de g(u0).

u=ZeroEuler(u0,t,m) ; print("Aproximação"); u     # Aproximadamente u(t)
print("Teste de minímo aproximado");g(u)                        # Teste de minímo aproximado

# Definição da aproximação de f(u).

fap<-function(s){a=u[1];b=u[2];c=u[3]
 a/(b+c*exp(-a*s) ) }

curve(fap,x[1],x[9])
points(x,y,col="red")        # Teste de ajuste.

plot(x,y-fap(x),'l')

print(" População estimada para 2021"); fap((2021-1930)/10)

# Se quiser pode aumentar t e n ou escolher outra condição inicial 
# para se convencer do resultado obtido.

hh=0.000001
Hessg<-function(u){          # aproximação de f''(x,y) ou da hessiana de f(x,y).
A=matrix(0,3,3)
v=u; v[1]=u[1]+hh;  w=u;w[1]=u[1]-hh
A[1,]=(gradg(v)-gradg(w))/(2*hh)
v=u; v[2]=u[2]+hh;  w=u;w[2]=u[2]-hh
A[2,]=(gradg(v)-gradg(w))/(2*hh)
v=u; v[3]=u[3]+hh;  w=u;w[3]=u[3]-hh
A[3,]=(gradg(v)-gradg(w))/(2*hh)
return(A)
}

Hessg(u)
det(Hessg(u))

u=u # condição inicial
m=20# número de iterações
h=1 # Tamanho do passo

for ( i in 1:m){
w=solve(Hessg(u),-h*gradg(u))
u=u+w
}

print("Aproximação par aponto de mínimo"); u
print("Aproximação para o valor mínimo"); g(u)

# Definição da aproximação de f(u).

fap<-function(s){a=u[1];b=u[2];c=u[3]
 a/(b+c*exp(-a*s) ) }

curve(fap,x[1],x[9])
points(x,y,col="red")        # Teste de ajuste.

plot(x,y-fap(x),'l')

print(" População estimada para 2021"); fap((2021-1930)/10)

[ + ] Show input
Absolute running time: 4.07 sec, cpu time: 5.44 sec, memory peak: 59 Mb, absolute service time: 4,3 sec