Exintegral-19-03-2021, 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                      # u1

for ( j in 1:m){
w=(B%*%u-bb)
                 # Método de Euler
alpha= t(w) %*%w/(t(B%*%w) %*%w);alpha=alpha[1]
u=u-alpha*w                      # 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,4)

h=1/5; x=seq(3,4,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=solve(A,y);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)
points(x,Erro(x),col="red")

Intpoli<-function(s){   #  Polinômio interpolador
q=u[1]*s
for (i in 2:n){q=q+u[i]*s^(i)/i}
q
}
Int34=Intpoli(4)-Intpoli(3)

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

h=1/5; x=seq(4,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 [5,7]");u=solve(A,y);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,4,5)
points(x,Erro(x),col="red")

Intpoli<-function(s){   #  Polinômio interpolador
q=u[1]*s
for (i in 2:n){q=q+u[i]*s^(i)/i}
q
}
Int45=Intpoli(5)-Intpoli(4)

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

h=1/5; x=seq(5,6,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 [5,7]");u=solve(A,y);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,5,6)
points(x,Erro(x),col="red")

Intpoli<-function(s){   #  Polinômio interpolador
q=u[1]*s
for (i in 2:n){q=q+u[i]*s^(i)/i}
q
}
Int56=Intpoli(6)-Intpoli(5)

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

h=1/5; x=seq(6,7,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 [5,7]");u=solve(A,y);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,6,7)
points(x,Erro(x),col="red")

Intpoli<-function(s){   #  Polinômio interpolador
q=u[1]*s
for (i in 2:n){q=q+u[i]*s^(i)/i}
q
}
Int67=Intpoli(7)-Intpoli(6)

Integral =Int34+Int45+Int56+Int67  # Aproximação da integral
print("Aproximação da integral de f(s) de 3 a 7"); Integral

[ + ] Show input
Absolute running time: 8.6 sec, cpu time: 11.68 sec, memory peak: 57 Mb, absolute service time: 8,88 sec