31-08-2020-Exemplo Ajuste, R

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(-5,5,by=0.5);x  
n=length(x)

y=c( 0.1617647, 0.1701807, 0.1818182, 0.1985294, 0.2236842, 0.2638889, 0.3333333,
 0.4659091, 0.7500000,  1.3750000, 2.0000000, 1.3750000, 0.7500000, 0.4659091, 0.3333333,
 0.2638889, 0.2236842, 0.1985294, 0.1818182, 0.1701807, 0.1617647)

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

# Suponha que y=f(x) (aproximado) e que f(s)=(a+bs^2)/(c+ds^2). Descobrir aproximações para a, b, c e d.
#
# (c+ds^2)f(s)=(a+bs^2) ou f(s)=-df(s)s^2/c+a/c+bs^2/c. Mudança de variável a_1=a/c, b_1=b/c, d_1=-d/c e 
# g_1(s)=f(s)s^2, e podemos escrever
#
# f(s)=a_1+b_1s^2+d_1 g_1(s). Determinar a_1, b_1 e d_1.

# Note que 
#            f(s)=(a+bs^2)/(c+ds^2)=(a/c+bs^2/c)/(1+ds^2/c)=(a_1+b_1s^2)/(1-d_1s^2).

# Preciso costruir g_1(x_i)=y_i*x_i^2 e minimizar a função

# |Au-y|^2, em que A matriz, cuja linha i é o vetor c(1,x_i^2, y_i*x_i^2).

# Resolver o sistema linear A^*A-A^*y=0 pelo método dos gradientes conjugados.

A=matrix(0,n,3); v=0*x;
for ( i in 1:n){A[i,]=c(1,x[i]^2,y[i]*x[i]^2);v[i]=y[i]}

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

gmv(c(1,6,3))     # Teste de g(u), com mudança de variável.

# Agora é conhecido que o gradiente de g(u) é dado por A^*(Au-v).
#
# Como o problema é linear, vamos resolver pelo 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(0,0,0)                       # 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)   
}
print("Soulução") ; u                                # Aproximação para u(t), sendo t a soma dos passos h.
A%*%u-v

gmv(u)

fapmv<-function(s){(u[1]+u[2]*s^2)/(1-u[3]*s^2)}  # Aproximação da função linearizada.

curve(fapmv,-5,5)                           # Teste de ajuste
points(x,y,col="red")

[ + ] Show input
Absolute running time: 0.72 sec, cpu time: 0.89 sec, memory peak: 45 Mb, absolute service time: 0,98 sec