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,40,200000)
#---------------------------------------------------------------------------------------

print("Soulução"); u                                # Aproximação para u que minimisa |Au-v|^2.
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
fapmv<-function(s){a/(b+c*exp(-a*(s-1930)/10))}  # Aproximação da função linearizada.

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

print("População estimada em 2021");fapmv(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))

t0=0  # tempo inicial
tf=1 # t final 
e0=u  # condição inicial 
n=4
h=(tf-t0)/n  # Tamanho do passo
Y=matrix(0,3,n)
Y[,1]=e0
for ( i in 1:(n-1)){
w=solve(Hessng(Y[,i]),-h*gradng(Y[,i]))
Y[,i+1]=Y[,i]+w
}

print("Aproximação para ponto de mínimo por Newton"); u=Y[,n];u
print("Aproximação para o valor mínimo"); g(Y[,n])

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

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

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

[ + ] Show input
Absolute running time: 4.17 sec, cpu time: 6.16 sec, memory peak: 59 Mb, absolute service time: 4,42 sec