n=4 //nombre d'oscillateurs
d=1.0/(n+1) // distance moyenne      
function [dq]=dK(p)
    dq=p
endfunction
		        

function [dp]=dV(q)
    dp(1)=(q(1)-d)+(q(1)-q(2)+d)
    for k=2:n-1,
        dp(k)=(q(k)-q(k-1)-d)+(q(k)-q(k+1)+d)
    end;
    dp(n)=(q(n)-q(n-1)-d)+(q(n)-1.0+d)
endfunction
		        
function e=energie(q,p)
    e=0.5*p(1)^2+0.5*(q(1)-d)^2;
    for k=2:n-1,
        e=e+0.5*p(k)^2+0.5*(q(k)-q(k-1)-d)^2;
    end;
    e=e+0.5*p(n)^2+0.5*(q(n)-q(n-1)-d)^2+0.5*(q(n)-1.0+d)^2;
endfunction

function [q,p]=pasVerlet(dK,dV,h,qn,pn)
    pi=pn-0.5*h*dV(qn)
    q=qn+h*dK(pi)
    p=pi-0.5*h*dV(q)
endfunction
		        

function [mq,mp,e,t]=verlet(q0,p0,T,n1,n2,dK,dV,energie)
    h=T/(n1*n2);
    s=size(q0);
    nv=s(1);
    mq = zeros(nv,n2+1);
    mp = zeros(nv,n2+1);
    e = zeros(1,n2+1);
    t=zeros(1,n2+1);
    mq(:,1)=q0;
    mp(:,1)=p0;
    q=q0; 
    p=p0;
    e(1)=energie(q0,p0);
    t(1)=0;
    for k2=1:n2,
        for k1=1:n1,
            [q,p]=pasVerlet(dK,dV,h,q,p)
        end;
        mq(:,k2+1)=q;
        mp(:,k2+1)=p;
        e(k2+1)=energie(q,p);
        t(k2+1)=t(k2)+n1*h;
    end;
endfunction
		    

function v=vitesse()
    v=(rand()-0.5)*0.1;
endfunction
q0=zeros(n,1);
p0=zeros(n,1);
for k=1:n,
    q0(k)=d*k;
    p0(k)=vitesse();
end;
n1=1;
n2=1000;
T=400;
[mq,mp,e,t]=verlet(q0,p0,T,n1,n2,dK,dV,energie);
		        
figA=scf();
plot2d(t,mq(1,:),style=color("black"),rect=[0,0,T,1]);
plot2d(t,mq(2,:),style=color("red"));
plot2d(t,mq(3,:),style=color("brown"));
plot2d(t,mq(4,:),style=color("orange"));
xtitle('','t','q');
		        
figB=scf();
plot2d(t,e);
xtitle('','t','Em');

function ydot=systeme(t,y)
    k=1;
    ydot(k)=y(k+n);
    ydot(k+n)=-(y(k)-d)-(y(k)-y(k+1)+d);
    for k=2:n-1,
        ydot(k)=y(k+n);
        ydot(k+n)=-(y(k)-y(k-1)-d)-(y(k)-y(k+1)+d);
    end;
    k=n;
    ydot(k)=y(k+n);
    ydot(k+n)=-(y(k)-y(k-1)-d)-(y(k)-1.0+d);
endfunction
y0=[q0;p0];
t0=0;
t=0:T/n2:T;
rtol=1e-5;
atol=1e-6;
[y,w,iw]=ode("adams",y0,t0,t,rtol,atol,systeme);
            
figC=scf();
plot2d(t,y(1,:),style=color("black"),rect=[0,0,T,1]);
plot2d(t,y(2,:),style=color("red"));
plot2d(t,y(3,:),style=color("brown"));
plot2d(t,y(4,:),style=color("orange"));
xtitle('','t','q');
		        
function e=energie2(y)
s=size(y);
e=zeros(1,s(2));
for j=1:s(2),
    e(j) = 0.5*y(1+n,j)^2+0.5*(y(1,j)-d)^2;
    for k=2:n-1,
        e(j)=e(j)+0.5*y(k+n,j)^2+0.5*(y(k,j)-y(k-1,j)-d)^2;
    end;
    e(j)=e(j)+0.5*y(n+n,j)^2+0.5*(y(n,j)-y(n-1,j)-d)^2+0.5*(y(n,j)-1.0+d)^2;
end;
endfunction
Em2=energie2(y);
		        
figD=scf();
plot2d(t,Em2);
xtitle('','t','Em');
iw(12)
n1=50;
n2=1000;
T=20000;
[mq,mp,e,t]=verlet(q0,p0,T,n1,n2,dK,dV,energie);
		    
figE=scf();
plot2d(t,e);
xtitle('','t','Em');
t=0:T/n2:T;
[y,w,iw]=ode('adams',y0,t0,t,rtol,atol,systeme);
Em2=energie2(y);

figF=scf();
plot2d(t,Em2);
xtitle('','t','Em');
