//
// equation du pendule simple non linéaire
// schémas symplectique vs Euler exp. et Heun (ordre 2)
//
T=100;
h=0.01;nt=int(T/h);
//
t=0:h:T;gsurl=1;
thetae=zeros(nt+1,1);
pe=zeros(nt+1,1);
thetas=zeros(nt+1,1);
ps=zeros(nt+1,1);
//
thetae(1)=%pi/2;pe(1)=0;
thetas(1)=%pi/2;ps(1)=0;
for i=1:nt
    thetae(i+1)=thetae(i)+h*pe(i);
    thetas(i+1)=thetas(i)+h*ps(i);
    pe(i+1)=pe(i)-h*gsurl*sin(thetae(i)) //Euler exp.
    ps(i+1)=ps(i)-h*gsurl*sin(thetas(i+1)) // symplectique
end
//
deff('ydot=f(t,y)','ydot=[y(2);-gsurl*sin(y(1))]')
//
yh=[%pi/2;0];
for i=1:nt
yend=yh(:,$)    
yh=[yh,yend+h/2*(f(t(i),yend)+f(t(i+1),yend+h*f(t(i),yend)))];
end
//
figure(1)
clf()
plot2d(t,thetas,1)
plot2d(t,thetae,2)
plot2d(t,yh(1,:),3)
legends(['symp','Euler exp.','Heun'],[1,2,3])
//
figure(2)
clf()
plot2d(thetas,ps,1)
plot2d(thetae,pe,2)
plot2d(yh(1,:),yh(2,:),3)
legends(['symp','Euler exp.','Heun'],[1,2,3])
