Gab1B-newton, 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-1930)/10
x   # Mudança de escala
n=length(x)

y=c(0.63,0.76,0.92,1.1,1.2,1.3,1.5,1.7,1.9);y

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

# Suponha que y=f(x) (aproximado) e que f(s)= ace^(as)/(1+bce^(as)) e f'(s)=(a-bf(s))f(s). Descobrir aproximações para a, b e c.
# 
# Note que f(0) nos dá c, caso conhecermos a e b. 
# 
# Para determinar os valores dos parâmetros resolvemos o problema de minimização 
#
#        min g(u)
# 
# em que g(u)= sum (y_i-f(x_i))^2, sendo u=(a,b,c(a,b)).
# 
# Para fazer isso usamos o método de Euler para resolver a equaçõ diferencial 
#
#      u'(t)=-grad g(u(t))
#      u(0)=u_0

# 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).
print(" Teste g"); g(c(1,2,1))

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

# Teste da função gradg(u)
print(" Teste gardiente de g"); gradng(c(1,2,1))

#---- 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*gradng(w)}
w
}

u0=c(1,1,1)

print("Teste condição inicial"); g(u0)
u=ZeroEuler(u0,20,20000)
#---------------------------------------------------------------------------------------

print("Soulução"); u                                
a=u[1]; b=u[2]; c=u[3]

print(" Resíduo da aproximação"); g(u)
print("Gradiente no ponto aproximado"); gradng(u)
print("Parâmetros a, b e c d encontrados"); a;b;c
fap<-function(s){a/(b+c*exp(-a*(s-1930)/10))}  # Aproximação da função.

print("Gráfico para ajuste")
curve(fap,1930,2010)                           # Teste de ajuste
points(seq(1930,2010,by=10),y,col="red")

print("População estimada em 2021");fap(2021)

# plot(tt,mf(Y[1,],Y[2,]),col = "blue",'l') # curva f(u(t),v(t)) para caminho de maximização obtido pelo método de Euler.

#-----------------------Segunda forma: método de Newton para maximizar e resolver a equação grad f(u)=0.

hh=0.000001
Hessng<-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,]=(gradng(v)-gradng(w))/(2*hh)
v=u; v[2]=u[2]+hh;  w=u;w[2]=u[2]-hh
A[2,]=(gradng(v)-gradng(w))/(2*hh)
v=u; v[3]=u[3]+hh;  w=u;w[3]=u[3]-hh
A[3,]=(gradng(v)-gradng(w))/(2*hh)
return(A)
}

Hessng(u)
det(Hessng(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(Hessng(u),-h*gradng(u))
u=u+w
}

print("Aproximação para ponto de mínimo por Newton"); u
print("Aproximação para o valor mínimo"); g(u)

a=u[1]; b=u[2]; c=u[3]

print("Gradiente no ponto aproximado"); gradng(u)
print("Hessiana e valores singulares no ponto aproximado");Hessng(u);eigen(Hessng(u))$values; eigen(solve(Hessng(u)))$values
print("Parâmetros a, b e c d encontrados"); a;b;c

fapb<-function(s){a/(b+c*exp(-a*(s-1930)/10))}  # Aproximação da função.

print("Gráfico para ajuste")
curve(fapb,1930,2010)                           # Teste de ajuste
points(seq(1930,2010,by=10),y,col="red")

print("População estimada em 2021");fapb(2021)

[ + ] Show input
Absolute running time: 1.12 sec, cpu time: 1.37 sec, memory peak: 55 Mb, absolute service time: 1,35 sec