## Exemple avec Heun et Crank-Nicholson

Cette partie :
1/ illustre les résultats du cours ;
2/ sert de démonstrateur pour les exercices de la seconde partie.

l'edo à intégrer est : 

dx/dt = f(x) ; x(t=0) =x0

Fonction réalisant un pas avec la méthode de Heun

En entrée : x0, fct=f et h le pas de temps. 

UnPasHeun(x0,fct,h):=block([a21:1/3,a32:2/3,b1:1/4,b3:3/4,
                            x01,x02,x03,f01,f02,f03],
  x01:x0,f01:fct(x01),
  x02:x0+h*a21*f01,f02:fct(x02),
  x03:x0+h*a32*f02,f03:fct(x03),
  x0+h*(b1*f01+b3*f03))$

Vérification que la fonction est correcte sur un exemple simple 
On doit retouver que Heun est d'ordre 3

si f(x)=x et x0=1 alors la solution est x=exp(t)
on compare exp(h) et x1 donné par un pas de Heun.

h:'h$
f(x):=x;
x1:expand(UnPasHeun(1,f,h));
taylor(exp(h),h,0,4);

Utilisation de la méthode pour l'exemple de l'oscillateur harmonique
dx/dt=-y ; dy/dt=x
g = x^2+y^2 doit se conserver.
Et on utilise un nombre de pas sur une période relativement petit, disons 10. 

Cette fois la dimension de x est de 2, c'est le vecteur (x,y) (on note x la première composante).

On prend x0=(1,0)

F(X):=[-X[2],X[1]]$
N:10$h:float(2*%pi/N)$
X0:[1,0]$g0:X0.X0;
for n:1 thru N do block([],
 X0:UnPasHeun(X0,F,h))$
X0;
gT:X0.X0;
gT-g0;

Fonction réalisant un pas de Crank-Nicholson sur l'exemple de l'oscillateur harmonique
dx/dt=-y ; dy/dt=x et toujours x0=(1,0)
Le résultat est presque parfait (mais c'est au prix de l'implicite)

x1:'x1$y1:'y1$x0:'x0$y0:'y0$h:'h$
eq:[x1=x0-h*(y1+y0)/2,y1=y0+h*(x1+x0)/2];
sol:solve(eq,[x1,y1]);
define(UnPasCNOH(X,h),expand(float(ev([x1,y1],sol,eval,x0=X[1],y0=X[2]))));

N:10$h:float(2*%pi/N)$
X0:[1,0]$g0:X0.X0;
for n:1 thru N do block([],
 X0:UnPasCNOH(X0,h))$
X0;
gT:X0.X0;
gT-g0;

## Exercices

### Faire comme pour l'exemple de Heun avec la 3/8

C'est à dire :  
* fabriquer la fonction qui fait un pas 
* tester qu'elle est correcte sur un exemple simple
* regarder ce qu'elle donne sur l'exemple de l'oscillateur harmonique

Vérification que la fonction est correcte sur un exemple simple 
On doit retouver que 3/8 est d'ordre 4

Utilisation de la méthode pour l'exemple de l'oscillateur harmonique
dx/dt=-y ; dy/dt=x
g = x^2+y^2 doit se conserver.
Et on utilise un nombre de pas sur une période relativement petit, disons 10. 

### Même question avec la DOPRI 5-4

DOPRI 5-4 est une méthode emboitée : Pour dx/dt =f(x) et x(t=0)=x0 
elle fournit deux estimations de x(h) qui sont x1 (à l'ordre 4) et xt1 (à l'ordre 5)

Le principe de la méthode est de comparer x1 et xt1 : 
si leur différence est trop grande c'est qu'il fallait prendre un pas plus petit (h/2)
sinon on allonge un peut (1.2*h) le pas pour les pas suivants.

La DOPRI 5-4

dopri54(x0,fct,h):=
block([a21:1/5,
       a31:3/40,a32:9/40,
       a41:44/45,a42:-56/15,a43:32/9,
       a51:19372/6561,a52:-25360/2187,a53:64448/6561,a54:-212/729,
       a61:9017/3168,a62:-355/33,a63:46732/5247,a64:49/176,a65:-5103/18656,
       a71:35/384,a72:0,a73:500/1113,a74:125/192,a75:-2187/6784,a76:11/84,
       b1:35/384,b2:0,b3:500/1113,b4:125/192,b5:-2187/6784,b6:11/84,b7:0,
       bc1:5179/57600,bc2:0,bc3:7571/16695,bc4:393/640,
           bc5:-92097/339200,bc6:187/2100,bc7:1/40,
       x01,x02,x03,x04,x05,x06,x07,
       f01,f02,f03,f04,f05,f06,f07,
       x1,xt1
       ],
x01:x0,
f01:fct(x01),
x02:x0+h*a21*f01,
f02:fct(x02),
x03:x0+h*(a31*f01+a32*f02),
f03:fct(x03),
x04:x0+h*(a41*f01+a42*f02+a43*f03),
f04:fct(x04),
x05:x0+h*(a51*f01+a52*f02+a53*f03+a54*f04),
f05:fct(x05),
x06:x0+h*(a61*f01+a62*f02+a63*f03+a64*f04+a65*f05),
f06:fct(x06),
x07:x0+h*(a71*f01+a72*f02+a73*f03+a74*f04+a75*f05+a76*f06),
f07:fct(x07),
x1:x0+h*(b1*f01+b2*f02+b3*f03+b4*f04+b5*f05+b6*f06+b7*f07),
xt1:x0+h*(bc1*f01+bc2*f02+bc3*f03+bc4*f04+bc5*f05+bc6*f06+bc7*f07),
[x1,xt1]
)$

Tester DOPRI 5-4 sur dx/dt=x x(t=0)=1 

Utilisation de la méthode pour l'exemple de l'oscillateur harmonique
dx/dt=-y ; dy/dt=x
g = x^2+y^2 doit se conserver.
Mais on ne contrôle pas le nombre de pas, seulement le temps d'arrivé. Il faut fabriquer une fonction
qui gère les pas de temps. 

### Mise sous forme ad hoc pour RK d'une EDO

On considère l'EDO
z'' + z'^2 + z = cos(t)
avec z(0) =0 ; z'(0) = 0

Que vaut z(1) ?

Il faut mettre le problème sous la forme dX/dt = F(X)