26-08-2020AjusteCurvaLinearizab, R

26-08-2020AjusteCurvaLinearizab

Language: Layout:

# 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")

[ + ] Show input
Absolute running time: 0.88 sec, cpu time: 0.82 sec, memory peak: 44 Mb, absolute service time: 1,1 sec

[1] "x" [1] -5.0 -4.5 -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 [16] 2.5 3.0 3.5 4.0 4.5 5.0 [1] "y" [1] 76.7717771 66.3145240 56.4127785 47.2264003 38.9028982 31.3924261 [7] 24.3614203 17.8432769 11.0583345 3.7080725 -0.1512441 4.0841671 [13] 11.2307420 17.5647328 24.0860429 31.2664318 39.0262100 47.3342189 [19] 56.5095134 66.1088441 76.7198095 [1] "z" [1] 673.716926 449.106206 290.752869 179.069911 102.197236 53.616249 [7] 26.679108 12.452959 6.101842 3.851295 1.000000 3.401579 [13] 5.934758 12.961091 27.699812 54.483496 100.714011 176.966332 [19] 287.943225 457.721305 676.695643 [1] "A" [,1] [,2] [,3] [1,] 625.0000 25.00 1 [2,] 410.0625 20.25 1 [3,] 256.0000 16.00 1 [4,] 150.0625 12.25 1 [5,] 81.0000 9.00 1 [6,] 39.0625 6.25 1 [7,] 16.0000 4.00 1 [8,] 5.0625 2.25 1 [9,] 1.0000 1.00 1 [10,] 0.0625 0.25 1 [11,] 0.0000 0.00 1 [12,] 0.0625 0.25 1 [13,] 1.0000 1.00 1 [14,] 5.0625 2.25 1 [15,] 16.0000 4.00 1 [16,] 39.0625 6.25 1 [17,] 81.0000 9.00 1 [18,] 150.0625 12.25 1 [19,] 256.0000 16.00 1 [20,] 410.0625 20.25 1 [21,] 625.0000 25.00 1 [1] "B" [,1] [,2] [,3] [1,] 1310401.039 61825.156 3166.625 [2,] 61825.156 3166.625 192.500 [3,] 3166.625 192.500 21.000 [1] "b" [,1] [1,] 1439107.369 [2,] 68534.525 [3,] 3603.066 [1] "u" [,1] [1,] 0.9928637 [2,] 2.0988790 [3,] 2.6192714 [1] "Erro na solução do sistema Au-v" [,1] [1,] 1.91412040 [2,] 3.15152738 [3,] -0.37843152 [4,] -1.74776649 [5,] -0.26609606 [6,] 0.90475370 [7,] 0.22149820 [8,] -0.08483782 [9,] -0.39082752 [10,] -0.64524990 [11,] 1.61927136 [12,] -0.19553421 [13,] -0.22374389 [14,] -0.59296955 [15,] -0.79920574 [16,] 0.03750618 [17,] 1.21712949 [18,] 0.35581275 [19,] 2.43121219 [20,] -5.46357179 [21,] -1.06459716 [1] "Teste de erro no ajuste linearizado" [,1] [1,] 30.25674 [1] "Grafico ajuste linearizado" [1] "Grafico ajuste des-linearizado"