Ex11-03-21EdsonAmanda, R

Language: Layout:

# Considere a função 
f<-function(s){sin(s)^3/log(3*s)}

curve(f,3,7)
# Queremos cakcular a integral de $f(s), no intervalo [3,7].
# Para isso vamos usar interpolação polinomial para aproximar f(s) em um conjunto de pontos e observar se o erro cometido na 
# aproximação e aceitável. Se assim o for faremos então a integral do polinomio e usaremos isso como se fosse a integral de f(s).

# ---- Para determinar o polinômio interpolador vamos resolver um sistema linear pelo método dos gradientes conjugados
#
#--- Inicio da função para o método dos gradientes conjugados
gradconj<-function(A,b,m){      
B=t(A)%*%A; bb=t(A)%*%b           # Adaptando o sistema para A*Au=A*v, ou Bu=b, em que B=A*A é positiva definida e b=A*v.
u=0*b
w=B%*%u-bb                        # w0
v=w                              # d1
alpha= t(w)%*%w /(t(B%*%w) %*%w); alpha=alpha[1]     # Escolha de h0
u=u-alpha*w                       # u1

for ( j in 1:m){
w=(B%*%u-bb)
beta=t(B%*%w) %*%v/(t(B%*%v) %*%v); beta=beta[1]
v=w-beta*v                   # Método de Euler
alpha= t(v) %*%w/(t(B%*%v) %*%v);alpha=alpha[1]
u=u-alpha*v                      # Método quase Euler (integral)
}
return(u)
} 
 
#----- Fim do método dos gradientes conjugados.

# Para reduzir erros na interpolação vamos dividir o intervalo [3,7] em dois outros intervalos, no caso [3,5] e [5,7].

# Aproximação polinomial em [3,5], como polinômio de grau 5.
curve(f,3,5)

h=2/5; x=seq(3,5,by=h); y=f(x)
points(x,y,col="red")

# Determinando os coeficientes do polinômio interpolador como solução do 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)
}
}

print("Aproximação para os coeficientes do polinômio interpolador em [3,5]");u=gradconj(A,y,50000);u

print("Erro na determinação dos coeficentes"); A%*%u-y
#---------------------------------------------------------------------------------------

poli<-function(s){   #  Polinômio interpolador
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")

Erro<-function(s){f(s)-poli(s)}
curve(Erro,3,4)

# Aproximação polinomial em [5,7], como polinômio de grau 5.
curve(f,5,7)

h=2/5; x=seq(5,7,by=h); y=f(x)
points(x,y,col="red")