# EDO: método de Euler.

# Limpando a memória
rm(list=ls(all=TRUE))

#-----------------------------------------------
RungK1<-function(f,u0,a,b,n){ # O método Euler
h=(b-a)/n
t=seq(a,b,by=h)
m=length(u0)
Y=matrix(0,m,n+1)
Y[,1]=u0
for (i in 1:n){ Y[,i+1]=Y[,i]+h*f(t[i],Y[,i]) }
return(Y)
}


#-----------------------------------------------
RungK2<-function(f,u0,a,b,n){ # O método Euler
h=(b-a)/n
t=seq(a,b,by=h)
m=length(u0)
Y=matrix(0,m,n+1)
Y[,1]=u0
for (i in 1:n){ 
k1=f(t[i],Y[,i])
k2=f(t[i]+h,Y[,i]+h*k1)
Y[,i+1]=Y[,i]+h*(k1+k2)/2 }
return(Y)
}

#----------Aplicação na Equação do pendulo simples

f<-function(s,U){ # Definicao de f(t,Y)=Y'.
u=U[1]; v=U[2]
return(c(v,-9.8*sin(u)))
}

a=0; b=1  # a é o tempo inicial e b o tempo final

u0=c(pi/6,0)  # condição inicial, posição do pendulo em radianos e velocidade inial em m/s

#--- Primeira simulação para estimar Y(1).

n=10

Y=RungK1(f,u0,a,b,n); print("Estimativa com n=10");Y[,n+1]

Y2=RungK2(f,u0,a,b,n); print("Estimativa com n=10");Y2[,n+1]


h=(b-a)/n
t=seq(a,b,by=h)
plot(t,Y[1,],'l',col = "blue",ylim=c(-1,1),ylab=expression(paste(theta,"(t)"))) 
par(new=TRUE)

plot(t,Y2[1,],'l',col = "blue",ylim=c(-1,1),ylab='',lty=2) 
par(new=TRUE)
#--- Segunda simulação para estimar Y(1).

n=100

V=RungK1(f,u0,a,b,n); print("Estimativa com n=100");V[,n+1]


h=(b-a)/n
t=seq(a,b,by=h)

plot(t,V[1,],'l',col = "green",ylab='',ylim=c(-1,1)) 
par(new=TRUE)

#--- Terceira simulação para estimar Y(1).

n=1000

V=RungK1(f,u0,a,b,n); print("Estimativa com n=1000");V[,n+1]


h=(b-a)/n
t=seq(a,b,by=h)

plot(t,V[1,],'l',col = "red",ylab='',ylim=c(-1,1)) 

par(new=TRUE)

#--- Quarta simulação para estimar Y(1).

n=10000

V=RungK1(f,u0,a,b,n); print("Estimativa com n=10000");V[,n+1]


h=(b-a)/n
t=seq(a,b,by=h)

plot(t,V[1,],'l',col = "purple",ylab='',ylim=c(-1,1)) 

legend(0.75,1, legend=c("n=10", "n=100","n=1000","n=10000","RK2-n=10"),
       col=c("blue", "green","red","purple","blue"), lty=c(1,1,1,1,2), cex=0.8)