clear();

//==============================================================================================================
//Modélisation de l'expérience FPUT (Fermi-Pasta-Ulam-Tsingou), première simulation numérique de l'Histoire
// (1953)
//==============================================================================================================

//==============================================================================================================
//Système différentiel, modélisant N masses reliées entre elles:
//  - N: Nombre de masses
//  - T: Durée de la simulation numérique
//  - ht: Pas de subdivision en temps
//  - alpha: Paramètre d'anharmonicité
//  - delta: Exposant d'anaharmonicité (si delta = 2: anaharmonicité quadratique, si delta = 3: anharmonicité
//cubique)
//  - Nprint: Nombre d'affichages de la solution souhaité
//==============================================================================================================

N = 64;
T = 5000;
ht = 0.5;
alpha = 0.5;
delta = 2;
Nprint = 500;


//==============================================================================================================
//Construction du champ de vecteurs
//==============================================================================================================

A = 2*diag(ones(1,N-2),0)-diag(ones(1,N-3),-1)-diag(ones(1,N-3),1);
D1 = -diag(ones(1,N-2),0)+diag(ones(1,N-3),1);
D2 = diag(ones(1,N-2),0)-diag(ones(1,N-3),-1);
D= D1+D2;

function YY = R(Y)
    //Y est un vecteur colonne à 2N-4 coordonnées
    X = Y(1:N-2);
    V = Y(N-1:2*N-4);
    XX = V;
    //VV = -A*X + alpha*( (D1*X).^delta - (D2*X).^delta ); // Equation sur la page Wikipédia en français de FPUT
    VV = -0.01*(A*X).*(ones(N-2,1)+alpha*D*X) // Equation sur la page Wikipédia en anglais de FPUT
    YY = [XX;VV];
endfunction


//==============================================================================================================
//Résolution numérique du système différentiel via la méthode de Runge-Kutta 4
//==============================================================================================================
Niter = int(T/ht); //Nombre d'itérations de la méthode numérique

Xsol = zeros(N-2,Niter);
X = sin(2*%pi*(2:N-1)/(N-2))'; //Position initiale des masses
Xsol(:,1) = X;


Vsol = zeros(N-2,Niter);
V = zeros(N-2,1); //Vitesse initiale des masses
Vsol(:,1) = V;

Ysol = [Xsol;Vsol];
Y = [X;V];

for n=1:Niter-1
    k1 = R(Y);
    k2 = R(Y+(ht/2)*k1);
    k3 = R(Y+(ht/2)*k2);
    k4 = R(Y+ht*k3);
    Y = Y + (ht/6)*(k1+k2+k3+k4);
    Ysol(:,n+1) = Y
end


//==============================================================================================================
//Tracé de la solution
//==============================================================================================================

Xsol = Ysol(1:N-2,:);
Xsol = [zeros(1,Niter);Xsol;zeros(1,Niter)];

clf();
xclick();
delete();
aa=get("current_axes")
aa.data_bounds=[0 min(Ysol) ; N max(Ysol)];
set(gca(),"grid",[1 1]);
title('Evolution de la solution (valeur absolue) à t = '+string(0));
for n=1:int(Niter/Nprint):Niter
    drawlater();
    b=get("current_entity");
    delete(b.children);
    plot((1:N)',Xsol(:,n),'s-','color','blue','thickness',2);
    title('Evolution de la solution à t = '+string(int(100*(n-1)*ht)/100));
    drawnow();
end
