Integracao-Interpolacaa, R

Integracao-Interpolacaa

# Calcule a integral da função f(s), no intervalo [0,1], com precisão de três casas decimais exatas.

f<-function(s){log(s^2*sin(s) + cos(s))}
curve(f,0,1)

x=seq(0,1,by=0.2);x
y=f(x);y

points(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.

# Vamos resolver o problema agora usando a ideia de Lagrange p(s)=y_0L_0(s)+y_1L_1(s)+...+y_(n-1)L_(n-1)(s)

n=length(x)

L<-function(i,s){
p=1
for (j in 1:n){
if ( i!=j){p=p*(s-x[j])/(x[i]-x[j])}
}
p
}

poliLag<-function(s){
p=0
for ( i in 1:n){p=p+y[i]*L(i,s)}
p
}

curve(poliLag,min(x),max(x))
points(x,y,col="red")

Erro<-function(s){f(s)-poliLag(s)}
curve(Erro,0,1)
points(x,Erro(x),col="red")

##----- Sistema linear e QR, por exemplo.

# 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

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

Intpoli-function(s){                                # Aproximação da integral 
q=u[1]*s
for (i in 2:n){q=q+u[i]*s^i/i}
q
}

print("Aproximação da integral"); Intpoli(1)-Intpoli(0)


run \| edit \| history \| help	0