26-08-2020AjusteCurvaLineariza, R

26-08-2020AjusteCurvaLineariza

# Ajuste de curvas

x=seq(-5,5,by=0.5);x
n=length(x)
erro=runif(n, min = -0.25, max = 0.25)

f<-function(s){20*s^2/(log(s^4+2*s^2+3))}

curve(f,-5,5)

y=f(x)+erro;y

points(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);z

fmv<-function(s){s^4+2*s^2+3}

curve(fmv,-5,5)
points(x,z,col="red")

gmv<-function(u){  #  |Au-v|^2 / 2
a=u[1];b=u[2];c=u[3]
p=0
for (i in 1:n){
p=p+( z[i]- (a*x[i]^4+b*x[i]^2+c) )^2}
p/2
}

gmv(c(1,6,3))

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]}

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

gmv2(c(1,6,3))

# para minimizar basta igualar o gradiente de gmv a zero e vai notar que isso 
# é um sistema linear

gradgmv<-function(u){
a=u[1];b=u[2];c=u[3]
p=c(0,0,0)
for (i in 1:n){
p=p+( z[i]- (a*x[i]^4+b*x[i]^2+c) )*c(x[i]^4,x[i]^2,1)}
p
}

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

n=10
u=c(2,2,2)                       # u0
w=B%*%u-b                        # w0=grag(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:(n-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)   
}
u                                # Aproximação para u(t), sendo t a soma dos passos h.
A%*%u-v

gmv2(u)

fapmv<-function(s){u[1]*s^4+u[2]*s^2+u[3]}

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]))}

curve(fap,-5,5,ylim=c(0,78))
points(x,y,col="red")


run \| edit \| history \| help	0