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

# Para resolver sistemas lineares

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

#-------Interpolação de f(s) no intervalo [0,1]
print("Intervalo [0,1]")

x=seq(0,1,by=0.1);x
y=f(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.

# 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)
A1=matrix(0,n,n)

for ( i in 1:n){
for (j in 1:n){
A1[i,j]=x[i]^(j-1)
}
}

A1

print(" Aproximação para u");u1=gradconj(A1,y,50);u1

print("Erro em Au-y"); A1%*%u1-y
#---------------------------------------------------------------------------------------

                         

poli1<-function(s){
u=u1
n=length(u)
q=u[1]
for (i in 2:n){q=q+u[i]*s^(i-1)}
q
}

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

Erro1<-function(s){f(s)-poli1(s)}
#curve(Erro1,min(x),max(x))


Intpoli1<-function(s){
u=u1
n=length(u)
q=u[1]*s
for (i in 2:n){q=q+u[i]*s^i/i}
q
}

#plot(Intpoli1,min(x),max(x))

#-------Interpolação de f(s) no intervalo [1,2]

print("Intervalo [1,2]")

x=seq(1,2,by=0.1);x; 
y=f(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.

# 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)
A2=matrix(0,n,n)

for ( i in 1:n){
for (j in 1:n){
A2[i,j]=x[i]^(j-1)
}
}

A2

print(" Aproximação para u");u2=gradconj(A2,y,50);u2

print("Erro em Au-y"); A2%*%u2-y
#---------------------------------------------------------------------------------------
                   
poli2<-function(s){
u=u2
n=length(u)
q=u[1]
for (i in 2:n){q=q+u[i]*s^(i-1)}
q
}

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

Erro2<-function(s){f(s)-poli2(s)}
#curve(Erro2,min(x),max(x))


Intpoli2<-function(s){
u=u2
n=length(u)
q=u[1]*s
for (i in 2:n){q=q+u[i]*s^i/i}
q
}

#plot(Intpoli2,min(x),max(x))



#-------Interpolação de f(s) no intervalo [2,2.5]-------------------------------------------

print("Intervalo [2,2.5]")

x=seq(2,2.5,by=0.05);x
y=f(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.

# 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)
A3=matrix(0,n,n)

for ( i in 1:n){
for (j in 1:n){
A3[i,j]=x[i]^(j-1)
}
}

A3

print(" Aproximação para u");u3=gradconj(A3,y,50);u3

print("Erro em Au-y"); A3%*%u3-y
#---------------------------------------------------------------------------------------
                   
poli3<-function(s){
u=u3
n=length(u)
q=u[1]
for (i in 2:n){q=q+u[i]*s^(i-1)}
q
}

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

Erro3<-function(s){f(s)-poli3(s)}
#curve(Erro3,min(x),max(x))


Intpoli3<-function(s){
u=u3
n=length(u)
q=u[1]*s
for (i in 2:n){q=q+u[i]*s^i/i}
q
}

#plot(Intpoli3,min(x),max(x))

#-------------------------------------------------------
#-------------Aproximação polinomial por partes
fap<-function(s){
p=poli1(s)
if ( s>1){p=poli2(s)}
if (s>2){ p=poli3(s)}
p
}

fap(0.2)
f(0.2)

fap(1.2)
f(1.2)

fap(2.2)
f(2.2)

#-----# Teste de ajuste polinomial por partes
x=seq(0,2.5,by=0.1)
y=f(x)
xx=seq(0,2.5,by=0.01)
Y=0*xx
for ( i in 1:length(Y)){Y[i]=fap(xx[i])}
plot(xx,Y,"l")             
points(x,y,col="red")


#---Integral da aproximação polinomial, por partes, no intervalo [0,2.5].
# Definição aproxomada para g(s).
Inapf<-function(s){
if (s<=1){p=Intpoli1(s)}
if (s>1){
if (s<=2){p=Intpoli1(1)+Intpoli2(s)-Intpoli2(1)}
else {p=Intpoli1(1)+Intpoli2(2)-Intpoli2(1)+Intpoli3(s)-Intpoli3(2)}
}
p
}

Inapf(2.5)


# Teste da estimativa da integral, para estimar o t da equação g(t)=1

xx=seq(0,2.5,by=0.01)
Y=0*xx
for ( i in 1:length(Y)){Y[i]=Inapf(xx[i])}
#plot(xx,Y,"l")                          

#-------Resolvendo equação g(s)=1, g(t)=int_0^tf(s)ds usando a aproximação polinomial.
# Defina G(s)=|g(s)-1|^2/2 a ser minimizada. 
# Isso nos dá grad G(s)=(g(s)-1)f(s)
G<-function(s){(Inapf(s)-1)^2/2}

xx=seq(0,2.5,by=0.01)
g=0*xx
for ( i in 1:length(Y)){g[i]=G(xx[i])}
plot(xx,g,"l") 

GradG<-function(s){(Inapf(s)-1)*f(s)}

# Método de Euler para resolver a equação u'(t)=-GradG(u(t)), u(0)=2.

u=2 # condição incial
h=0.1
it=100
for (i in 1:it){
u=u-h*GradG(u)
}

print("Solução da equação  g(s)=1, g(t)=int_0^tf(s)ds");u
G(u)
GradG(u)
Inapf(u)