16-09-2020ExVitorc, 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=c(2,5,8,11,14,17,27,31,41,44)
n=length(x)

y=c(94.8,89.7,81.3,74.9,68.7,64.0,49.3,44.0,39.1,31.6)

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

# Suponha que y=f(x) (aproximado) e que f(s)= s/(a+bs). Descobrir aproximações para 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).
# 
# 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]
p=0
for (i in 1:n){
p=p+( y[i]-1/(a+b*x[i]) )^2}
p/2
}

gradg<-function(u){   # função gradiente de g(u).
p=0*u
a=u[1];b=u[2]
m=length(u)
for ( i in 1:n){p=p+( y[i]-1/(a+b*x[i]) )*c(1/(a+b*x[i])^2,x[i]/(a+b*x[i]))}
p
}

#---- 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*gradg(w)}
w
}

# Teste de aproximação de solução do problema de minimização. 
u0=c(1,2); t=10; m=10000000
print("Erro na condição inicial") ; g(u0)                     # Teste de g(u0).

u=ZeroEuler(u0,t,m) ; print("Solução"); u     # Aproximadamente u(t)
print("Erro na aproximação"); g(u)                        # Teste de minímo aproximado
gradg(u)

# Definição da aproximação de f(u).

fap<-function(s){ 1/(u[1]+u[2]*s)}

curve(fap,2,44)
points(x,y,col="red")        # Teste de ajuste.

# Se quiser pode aumentar t e n ou escolher outra condição inicial 
# para se convencer do resultado obtido.

# Ajuste linearizado 
#
# Mudança de variável z=1/y=~ 1/(f(x)) = a+bx. Determinar a e b.

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

# Resolver o sistema linear A^*(Au-z)=A^*Au-A^*z=0 pelo método A=QR.

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

# Agora é conhecido que o gradiente de g(u) é dado por A^*(Au-v).
#
# Como o problema é linear, vamos resolver pelo Método de decomposição A=QR.
# Isso nos leva a resolver a equação Ru=t(Q)v.

# ---------------Processo de ortogonalização de Gram-Schmdit

pe<-function(u,v){ # Dados os vetores u e v, somo o produto de cada coordenada de u com a respectiva cooordenada de v.
n=length(u)
p=u[1]*v[1]
for ( i in 2:n){p=p+u[i]*v[i]}
p
}

mod<-function(u){sqrt(pe(u,u))}  # módulo de u

proj<-function(u,v){pe(u,v)*v}  # Projeção de u sobre v (ortogonal), com norma de v igual 1

GramS<-function(A){  # Ortogonliza as colunas de A
n=length(A[1,])

Q=0*A

Q[,1]=A[,1]/mod(A[,1])

for (i in 2:n){
p=A[,i]
for (j in 1:(i-1)){
p=p-proj(A[,i],Q[,j])
}
Q[,i]=p/mod(p)
}

Q
}

Q=GramS(A); print("Matriz Q"); Q

t(Q)%*%Q      # Teste de ortogonalização

R=t(Q)%*%A; print("Matriz R"); R

#------ Resolução de sistema triangular Ru=t(Q)v.

TSup<-function(A,b){
p=0*b; n=length(b)
p[n]=b[n]/A[n,n]
for (i in (n-1):1){
q=b[i]
for (j in n:(i+1)){
q=q-A[i,j]*p[j]
}
p[i]=q/A[i,i]
}
p
}

v=t(Q)%*%(1/y);v
u=TSup(R,v)
#---------------------------------------------------------------------------------------

print("Soulução") ; u                                # Aproximação para u que minimisa |Au-v|^2.

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

curve(fapmv,2,44)                           # Teste de ajuste
points(x,1/y,col="red")

[ + ] Show input
Absolute running time: 15.08 sec, cpu time: 22.91 sec, memory peak: 54 Mb, absolute service time: 15,32 sec