Ex11-03-21EdsonAmanda, R

Ex11-03-21EdsonAmanda

f<-function(s){sin(s)^3/log(3*s)}

curve(f,3,7)

curve(f,3,5)

h=(5-3)/5;print("h");h

x=seq(3,5,by=h);print("x");x

y=f(x);print("y=f(x)");y
points(x,y,col="red")

AproxInt=sum(y)*h    # (y_1+y2+....+y41)h = aproximação da integral por soma de Riemann.
print("Aproximação integral");AproxInt # 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)
}
}

gradconj<-function(A,b,m){      # Início do método
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)
}

print(" Aproximação para u");u=gradconj(A,y,50000);u

print("Erro em Au-y"); A%*%u-y
#---------------------------------------------------------------------------------------

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

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

curve(f,5,7)

h=2/5;print("h");h

x=seq(5,7,by=h);print("x");x