26-08-2020AjusteCurvaLinearizab, R

26-08-2020AjusteCurvaLinearizab

# Ajuste de curvas

x=seq(-5,5,by=0.5)
print("x");x

n=length(x)
y=c(76.7717771, 66.3145240, 56.4127785, 47.2264003, 38.9028982, 31.3924261, 24.3614203, 17.8432769, 11.0583345, 3.7080725, -0.1512441, 4.0841671, 11.2307420, 17.5647328, 24.0860429, 31.2664318, 39.0262100, 47.3342189, 56.5095134, 66.1088441, 76.7198095)

print("y");y

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

#---------- Mudança de variável z=e^(20x^2/y)=ax^4+bx^2+c

z=exp(20*x^2/y); print("z");z

plot(x,z,col="red")  # Teste de alinhamento

A=matrix(0,n,3); v=0*x

for ( i in 1:n){A[i,]=c(x[i]^4,x[i]^2,1);v[i]=z[i]}

print("A"); A

glim<-function(u){  #  |Au-v|^2 / 2
w=A%*%u-v
t(w)%*%w/2}

# Método dos gradientes conjugados 
B=t(A)%*%A; print("B");B; b=t(A)%*%v;print("b");b      # Adaptando o sistema para A*Au=A*v, ou Bu=b, em que B=A*A é positiva definida e b=A*v.

m=5                              # Em tese bastam 3 iterações
u=c(2,2,2)                       # u0
w=B%*%u-b                        # w0=graglim(u0)
d=w                              # d1
h= t(w)%*%w /(t(B%*%w) %*%w)     # Escolha de h0
u=u-h[1]*w                       # u1
w=B%*%u-b

for ( j in 1:(m-1)){
a= t(B%*%d) %*%w/(t(B%*%d) %*%d) # Escolha de alpha
d=w-a[1]*d                       # w - projeção ortogonal de w na direção de d
bt= t(d) %*%w/(t(B%*%d) %*%d)    # Escolha de beta
u=u-bt[1]*d                       # (quase) Método de Euler u_{n+1}=u_n-h (w_n-proj_d(w_n))
w=B%*%u-b                        # gradg(u_n)   
}
print("u"); u                                # Aproximação para u(t), sendo t a soma dos passos h.
print("Erro na solução do sistema Au-v");A%*%u-v

print("Teste de erro no ajuste linearizado"); glim(u)

fapmv<-function(s){u[1]*s^4+u[2]*s^2+u[3]} # Aproximação para f linearizada

print("Grafico ajuste linearizado")
curve(fapmv,-5,5)
points(x,z,col="red")

fap<-function(s){20*s^2/(log(u[1]*s^4+u[2]*s^2+u[3]))} # Aproximação para f

print("Grafico ajuste des-linearizado")
curve(fap,-5,5,ylim=c(0,78))
points(x,y,col="red")


run \| edit \| history \| help	0