Double Pendulum with euler (copy and past on your computer), R

Double Pendulum with euler (copy and past on your computer)

# EDO: método de Euler.

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

#-----------------------------------------------
Euler<-function(f,u0,a,b,n){

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

Y=matrix(0,m,n+1)
Y[,1]=u0

# O método Euler
for (i in 1:n)
{
 Y[,i+1]=Y[,i]+f(t[i],Y[,i])*h
}

Y
}

#--------------------------------------------------------
# Equação do pendulo duplo acoplado.

a=0

b=100

u0=c(pi/12,pi/12,2,2)

# Dados

f<-function(s,U) # Definicao de f(t,Y)=Y'.
{
a=U[1]
b=U[2]
u=U[3]
v=U[4]
dv=( (2/5)*cos(a-b)*(sin(a-b)*v^2+5*sin(a))+sin(a-b)*u^2-2*sin(b) )/(1-(2/5)*(cos(a-b))^2)
du= -(2/5)*(dv*cos(a-b)+sin(a-b)*v^2+5*sin(a))
p=c(u,v,du,dv)
p
}

n=5000

Y=Euler(f,u0,a,b,n)

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

plot(t,Y[2,],'l',col = "blue",ylim=c(-10,40)) 
points(t,Y[1,],'l',col = "red")

# x(t) e y(t) movimento do pêndulo duplo no plano.
x1=5*sin(Y[1,])
y1=10-5*cos(Y[1,])
plot(x1,y1,'l')

x2=5*sin(Y[1,])+5*sin(Y[2,])
y2=10-5*cos(Y[1,])-5*cos(Y[2,])

plot(x2[1],y2[1],'l',xlim=c(-10,10),ylim=c(0,20))

for (i in 100:n){
x2a=x2[(i-100):i]
y2a=y2[(i-100):i]
x1a=x1[(i-100):i]
y1a=y1[(i-100):i]

plot(x2a,y2a,'l',xlim=c(-10,10),ylim=c(0,20))
points(x1a,y1a,'l',col="gray")
x2b=x2[1:3]
y2b=y2[1:3]

x2b[1]=x2[i]
y2b[1]=y2[i]
x2b[2]=x1[i]
y2b[2]=y1[i]
x2b[3]=0
y2b[3]=10
points(x2b,y2b,'l',col="blue")
points(x2b,y2b,pch = 19,col="red")

}


run \| edit \| history \| help	0