rungeefeito, R - rextester

rungeefeito

f<-function(s){1/(1+s^2)}
curve(f,-5,5)

x=seq(-5,5,by=1); x
y=f(x);y

points(x,y,col="red")  # 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,-5,5)
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+a_5s^5 é 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)
u=TSup(R,v)
#---------------------------------------------------------------------------------------

print("Soulução") ; u

poli<-function(s){ #  p(s)=a_0+a_1s+a_2s^2+a_3s^3+a_4s^4+a_5s^5+a_6s^6
q=u[1]
for (i in 2:n){q=q+u[i]*s^(i-1)}
q
}

Intpoli<-function(s){ #  p(s)=a_0s+a_1s^2/2+a_2s^3/3+a_3s^4/4+a_4s^5/5+a_5s^6/6+a_6s^7/7
q=u[1]*s
for (i in 2:n){q=q+u[i]*s^(i)/i}
q}

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

Integral=Intpoli(5)-Intpoli(-5); Integral

Trap<-function(f,v){
n=length(v)
trap=0
for ( i in 1:(n-1)){trap=trap+(f(v[i])+f(v[i+1]))*(v[i+1]-v[i])/2}
trap
}

Trap(f,seq(-5,5,by=0.05))


run \| edit \| history \| help	0