xintf-rk4-4x4.tex

%% ========================================================================== %%
%% Méthode de Runge Kutta à l'ordre 4 avec adaptation du pas                  %%
%% ========================================================================== %%
%% Résolution d'un système différentiel 4x4 par la méthode de Runge-Kutta à   %%
%% pas variable à l'ordre 4 avec finalisation à l'ordre 5                     %%
%% ========================================================================== %%
 
\def\rkSeuils{
  \xintdeffloatvar epsilon := 0.0001;%
  \xintdeffloatvar seuil := (4/0.9)^(-5);%
  \xintdeffloatvar stop := 0.00001;%
}
 
\def\rkHillOrbit{
  \xintdeffloatvar epsilon := 0.000001;%
  \xintdeffloatvar consistance := (4/0.9)^(-5);%
  \xintdeffloatvar stop := 0.00001;%
  %
  \xintdeffloatvar m := 1/82.45;%
  \xintdeffloatvar M := 1-m;%
  %
  \xintdeffloatfunc r32(x,y,c) := ((x+c)^2+y^2)^(1.5);%
  \xintdeffloatfunc sm(t,qx,qy,px,py) := px+qy,py-qx,%
                        py-M(qx+m)/r32(qx,qy,m)-m(qx-M)/r32(qx,qy,-M),%
                        -px-qy*(M/r32(qx,qy,m)+m/r32(qx,qy,-M));%
}                                       
 
\def\rkBaseIV{%
  \xintdeffloatfunc mul(a,b,c,d,x) := a*x,b*x,c*x,d*x;%
  \xintdeffloatfunc div(a,b,c,d,x) := a/x,b/x,c/x,d/x;%
  \xintdeffloatfunc add(ra,rb,rc,rd,sa,sb,sc,sd) := ra+sa,rb+sb,rc+sc,rd+sd;%
  \xintdeffloatfunc ck(t,k1,k2,k3,k4,p1,p2,p3,p4) := sm(t,p1+k1,p2+k2,p3+k3,p4+k4);%
  \xintdeffloatfunc norme1(ra,rb,rc,rd,sa,sb,sc,sd) := abs(ra-sa)+abs(rb-sb)+abs(rc-sc)+abs(rd-sd);%
  \xintdeffloatvar  zero := 0,0,0,0;%
}  
 
\def\rkStep {
 % t, h, p doivent être affectés
 \xintdeffloatvar ka := mul(ck(t,zero,p),h);%
 \xintdeffloatvar kb := mul(ck(t+h/2,div(ka,2),p),h);%
 \xintdeffloatvar kc := mul(ck(t+h/2,div(kb,2),p),h);%
 \xintdeffloatvar kd := mul(ck(t+h,kc,p),h);%
 \xintdeffloatvar newp := add(p,add(div(ka,6),add(div(kb,3),add(div(kc,3),div(kd,6)))));%
 % newp est la valeur servie en retour
}
 
\def\rkProcIV {
  % ta : temps actuel, ha : pas actuel, pa : position actuelle
  \xintdeffloatvar t,h := ta,ha/2;%
  \xintdeffloatvar p   := pa;%
  \rkStep{}%
  \xintdeffloatvar t,h := t+ha/2,ha/2;%
  \xintdeffloatvar p   := newp;%
  \rkStep{}%
  \xintdeffloatvar pb  := newp;%
  \xintdeffloatvar t,h := ta,ha;%
  \xintdeffloatvar p   := pa;%
  \rkStep{}%
  \xintdeffloatvar s  := norme1(pb,newp)/epsilon;%
  \xintifboolexpr{s>1}%
    {%
      \xintdeffloatvar ha := 0.9s^(-0.25)ha;%
      \xintifboolexpr{ha<stop}{\xintdeffloatvar exit := 1;}{\rkProcIV}%
    }%
    {%
      \xintdeffloatvar  newqx := (16pb[0]-newp[0])/15;%
      \xintdeffloatvar  newqy := (16pb[1]-newp[1])/15;%
      \xintdeffloatvar  newpx := (16pb[2]-newp[2])/15;%
      \xintdeffloatvar  newpy := (16pb[3]-newp[3])/15;%
      \xintdeffloatvar  pa := newqx,newqy,newpx,newpy;%
      \xintdeffloatvar  QX  := QX,newqx;%
      \xintdeffloatvar  QY  := QY,newqy;%
      \xintdeffloatvar  PX  := PX,newpx;%
      \xintdeffloatvar  PY  := PY,newpy;%
      \xintdeffloatvar  ta  := ta+ha;%
      \xintdeffloatvar  LT  := LT,ta;%
      \xintdeffloatvar  H   := H,ha;%
      \xintdefvar N := N+1;%
      \xintdeffloatvar hmin := (ha<hmin)?{ha}{hmin};%;
      \xintdeffloatvar hmax := (ha>hmax)?{ha}{hmax};%;
      \xintifboolexpr{s>seuil}%
        {\xintdeffloatvar ha := 0.99s^(-0.2)ha;}%
        {\xintdeffloatvar ha := 4ha;}%
      \xintifboolexpr{ta<tf}{\rkProcIV}{\xintdeffloatvar exit:= 0;}%
    }%
}
 
\def\rkSolIV(#1,#2,#3,#4,#5,#6,#7,#8){
  \xintdeffloatvar ta := #1;%
  \xintdeffloatvar tf := #2;%
  \xintdeffloatvar ha := #3;%
  \xintdeffloatvar QX := #4;%
  \xintdeffloatvar QY := #5;%
  \xintdeffloatvar PX := #6;%
  \xintdeffloatvar PY := #7;%
  \xintdeffloatvar pa := QX,QY,PX,PY;%
  \xintdeffloatvar LT := ta;%
  \xintdeffloatvar H  := ha;%
  \xintdeffloatvar hmin := ha;%
  \xintdeffloatvar hmax := ha;%
  \IfFileExists{#8_par.xint}{%
    % Pour l'instant, rien à faire, mais que pourrait-on faire ?
  }{%
    \xintdefvar N := 1;%
    \rkBaseIV{}
    \rkProcIV{}
    \WriteXintValues{#1,#2,#3,#4,#5,#6,#7,N,hmin,hmax}{#8_par.xint}
    \WriteXintValues{LT}{#8_t.xint}
    \WriteXintValues{QX}{#8_qx.xint}
    \WriteXintValues{QY}{#8_qy.xint}
    \WriteXintValues{PX}{#8_px.xint}
    \WriteXintValues{PY}{#8_py.xint}
    \WriteXintValues{H}{#8_h.xint}
  }%
  % La suite se déroulera toujours avec les données enregistrées...
  % Plus de séquences mais des variables préfixées !
  \ReadXintValues{ct}{#8_par.xint}
  \ReadXintValues{QX}{#8_qx.xint}
  \ReadXintValues{QY}{#8_qy.xint}
  \ReadXintValues{PX}{#8_px.xint}
  \ReadXintValues{PY}{#8_py.xint}
  \ReadXintValues{LT}{#8_t.xint}
  \ReadXintValues{H}{#8_h.xint}
  \xintdeffloatvar N,hmin,hmax := ct8,ct9,ct10;%
}
\endinput