21-09-2020Intpoli, R

21-09-2020Intpoli

Language: Layout:

x=c(9,8,1,6,10,5,3);x

y=c(3,2,1,5,12,7,13);y

plot(x,y)  # Supomos y=f(x) e que f(x_i)=p(x_i), em que p(s)=a_0+a_1s+a_2s^2+a_3s^3+a_4s^4 é um polinômio.

# Podemos resolver esse problema resolvendo o sistema linear 
# Au=y
# em que a matriz A tem coluna i dada pelo vetor (1,x_i,x_i^2,...,x_i^(n-1)).

n=length(x)
A=matrix(0,n,n)

for ( i in 1:n){
for (j in 1:n){
A[i,j]=x[i]^(j-1)
}
}

# A=QR, então Au=y torna-se Ru=t(Q)y
# ---------------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)%*%(y);v
u=TSup(R,v)
#---------------------------------------------------------------------------------------

print("Soulução") ; u

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

poli<-function(s){
q=u[1]
for (i in 2:n){q=q+u[i]*s^(i-1)}
q
}

curve(poli,min(x),max(x))                           # Teste de ajuste
points(x,y,col="red")

[ + ] Show input
Absolute running time: 0.85 sec, cpu time: 0.93 sec, memory peak: 45 Mb, absolute service time: 1,11 sec

[1] 9 8 1 6 10 5 3 [1] 3 2 1 5 12 7 13 [,1] [,2] [,3] [,4] [,5] [,6] [,7] [1,] 1 9 81 729 6561 59049 531441 [2,] 1 8 64 512 4096 32768 262144 [3,] 1 1 1 1 1 1 1 [4,] 1 6 36 216 1296 7776 46656 [5,] 1 10 100 1000 10000 100000 1000000 [6,] 1 5 25 125 625 3125 15625 [7,] 1 3 9 27 81 243 729 [1] "Matriz Q" [,1] [,2] [,3] [,4] [,5] [,6] [1,] 0.3779645 3.750000e-01 0.1152731 -0.2389979 -0.4178417 -0.57304218 [2,] 0.3779645 2.500000e-01 -0.1667690 -0.4895528 -0.2643172 0.43368617 [3,] 0.3779645 -6.250000e-01 0.5617140 -0.3417617 0.1796821 -0.04232741 [4,] 0.3779645 1.110223e-16 -0.4412698 -0.1455581 0.4183840 0.31583339 [5,] 0.3779645 5.000000e-01 0.4938430 0.4642860 0.3228672 0.18968111 [6,] 0.3779645 -1.250000e-01 -0.4337286 0.1924501 0.3363891 -0.54505026 [7,] 0.3779645 -3.750000e-01 -0.1290628 0.5591344 -0.5751635 0.22121918 [,7] [1,] -0.37833966 [2,] 0.51886584 [3,] 0.01080971 [4,] -0.60534352 [5,] 0.08647763 [6,] 0.45400765 [7,] -0.08647765 [,1] [,2] [,3] [,4] [,5] [1,] 1.000000e+00 3.056405e-16 1.339827e-16 -1.715786e-15 1.447636e-14 [2,] 3.056405e-16 1.000000e+00 -1.277208e-15 4.891547e-15 -1.098567e-14 [3,] 1.339827e-16 -1.277208e-15 1.000000e+00 1.992157e-14 -9.920168e-14 [4,] -1.715786e-15 4.891547e-15 1.992157e-14 1.000000e+00 -9.859770e-13 [5,] 1.447636e-14 -1.098567e-14 -9.920168e-14 -9.859770e-13 1.000000e+00 [6,] -2.216378e-13 -5.383513e-14 6.063473e-13 1.288494e-11 1.350645e-10 [7,] 2.562204e-12 1.212543e-12 -3.455923e-12 -1.112806e-10 -2.087968e-09 [,6] [,7] [1,] -2.216378e-13 2.562204e-12 [2,] -5.383513e-14 1.212543e-12 [3,] 6.063473e-13 -3.455923e-12 [4,] 1.288494e-11 -1.112806e-10 [5,] 1.350645e-10 -2.087968e-09 [6,] 1.000000e+00 -4.722976e-08 [7,] -4.722976e-08 1.000000e+00 [1] "Matriz R" [,1] [,2] [,3] [,4] [,5] [1,] 2.645751e+00 1.587451e+01 1.194368e+02 9.864873e+02 8.564675e+03 [2,] 7.771561e-16 8.000000e+00 8.925000e+01 8.750000e+02 8.375250e+03 [3,] 3.885781e-16 -8.493206e-15 2.071943e+01 3.400380e+02 4.158792e+03 [4,] -4.662937e-15 1.199041e-14 6.481482e-13 4.677608e+01 1.046173e+03 [5,] 3.830269e-14 1.419975e-13 -1.303735e-12 -7.518397e-11 1.106292e+02 [6,] -5.864198e-13 -3.948952e-12 -1.871614e-11 5.431633e-10 2.859463e-08 [7,] 6.778980e-12 5.037469e-11 3.426395e-10 -2.791828e-09 -3.296816e-07 [,6] [,7] [1,] 7.671243e+04 701727.3287 [2,] 7.985300e+04 762599.2500 [3,] 4.590886e+04 483927.9137 [4,] 1.587942e+04 205562.5214 [5,] 3.117420e+03 55876.8630 [6,] 1.479241e+02 5211.5678 [7,] -1.512784e-05 217.9233 [,1] [1,] 16.25247234 [2,] 1.25000000 [3,] -0.42015434 [4,] 11.42167876 [5,] -0.75855259 [6,] 2.02175619 [7,] -0.02161951 [1] "Soulução" [,1] [1,] -3.728572e+01 [2,] 5.852977e+01 [3,] -2.452382e+01 [4,] 4.703077e+00 [5,] -4.403773e-01 [6,] 1.716272e-02 [7,] -9.920696e-05