Runge Kutta 4 à pas variables [Résolu]

[aaaallleex]

Bonjour à tous,
Je voudrais resoudre les équations du pendule double grâce à la méthode de Runge Kutta à pas variable.
Mais je n'arrive pas à definir et à utiliser correctement le pas pour chaque étapes.
je m'explique: j'ai d'abord créé une boucle pour un pas constant ici 0.001 qui marche très bien.
L'utilisation d'un pas variable me ferai gagner du temps de calcul qui me permettrais de résoudre l'équation sur un temps plus long.

Comment prendre en compte le pas variable?

L'idée est la suivante : après avoir obtenus om1 et om2 avec un pas h, on les compares avec om1_1 et om2_1 qui on été calculés avec un pas h/2 et qui sont donc plus précis.
Si la difference est inférieur à la précision demandée, on garde le pas, sinon on le reduit et on recommence.
Voir
théorie :http://media4.obspm.fr/public/m2r/co...on3_2_4_6.html
ou théorie et code: http://www.physique.usherbrooke.ca/p...HQ405_ipad.pdf (p61)
ou encore http://users-phys.au.dk/fedorov/Numeric/10/odes.pdf
et enfin un chapitre de Numerical recipes ici: http://www.pit.physik.uni-tuebingen....mrec/c16-2.pdf.

J'ai donc tenter de rajouter une boucle Runge Kutta plus précise dans la première, de comparer les resultats et ensuite de modifier le pas si necassaire.

Voici le code :

Code :

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
# -*- coding: utf-8 -*-
from numpy import*
from scipy import*
import matplotlib.pylab as plt
 
#Condition initiales
l1,l2,m1,m2=1,1,1,1                         #Longueur et masses des pendules 
t0,tf= 0,2                                  #Temps de l'experience
th1_0,th2_0,om1_0,om2_0=0,3,0,0             #Angles et vitesses initiales
g=10
 
#Les deux equations du pendule double
dom1=lambda t,om1,om2,th1,th2 : (-m2*l1*sin(th1-th2)*cos(th1-th2)*om1**2-m2*l2*sin(th1-th2)*om2**2-m1*g*sin(th1-th2)*cos(th2))/l1*(m1+m2*(sin(th1-th2))**2)
dom2=lambda t,om1,om2,th1,th2 : ((m1+m2)*l1*sin(th1-th2)*om1**2+m2*l2*sin(th1-th2)*cos(th1-th2)*om2**2+(m1+m2)*g*sin(th1-th2)*cos(th1))/l2*(m1+m2*(sin(th1-th2))**2)
 
hmin,h,hmax=0.01,0.1,0.5                                             #Pas initial
t,th1,th2,om1,om2 = t0,th1_0,th2_0,om1_0,om2_0      #Initialisation des variables
th1_1,th2_1,om1_1,om2_1=th1_0,th2_0,om1_0,om2_0                                                
p=0.001                                             #Precision demandée
 
#Tableaux des données
T=array([t])                                  
Th1=array([th1])
Th2=array([th2])
Om1=array([om1])
Om2=array([om2])
 
#Boucle Runge-Kutta 4
 
while t<tf:                                  
 
    if t+h>tf: h=tf-t
 
    th1 += h*om1
    th2 += h*om2
 
    k1 = h*dom1 (t,om1,om2,th1,th2)
    j1 = h*dom2 (t,om1,om2,th1,th2)
    k2 = h*dom1 (t+h/2,om1+h*k1/2,om2+h*j1/2,th1,th2)
    j2 = h*dom2 (t+h/2,om1+h*k1/2,om2+h*j1/2,th1,th2)
    k3 = h*dom1 (t+h/2,om1+h*k2/2,om2+h*j2/2,th1,th2)
    j3 = h*dom2 (t+h/2,om1+h*k2/2,om2+h*j2/2,th1,th2)
    k4 = h*dom1 (t+h,om1+h*k3,om2+h*j2,th1,th2)
    j4 = h*dom2 (t+h,om1+h*k3,om2+h*j2,th1,th2)
 
    om1 += k1/6 + k2/3 + k3/3 + k4/6
    om2 += j1/6 + j2/3 + j3/3 + j4/6
 
    h=h/2                                  #Division du pas par 2  
 
    for k in range (2):                      #Double boucle Runge Kutta avec un pas divisé par 2
 
        th1_1 += h*om1_1                     #===> Meilleur precision 
        th2_1 += h*om2_1
 
        k1_1 = h*dom1 (t,om1_1,om2_1,th1_1,th2_1)
        j1_1 = h*dom2 (t,om1_1,om2_1,th1_1,th2_1)
        k2_1 = h*dom1 (t+h/2,om1_1+h*k1_1/2,om2_1+h*j1_1/2,th1_1,th2_1)
        j2_1 = h*dom2 (t+h/2,om1_1+h*k1_1/2,om2_1+h*j1_1/2,th1_1,th2_1)
        k3_1 = h*dom1 (t+h/2,om1_1+h*k2_1/2,om2_1+h*j2_1/2,th1_1,th2_1)
        j3_1 = h*dom2 (t+h/2,om1_1+h*k2_1/2,om2_1+h*j2_1/2,th1_1,th2_1)
        k4_1 = h*dom1 (t+h,om1_1+h*k3_1,om2_1+h*j2_1,th1_1,th2_1)
        j4_1 = h*dom2 (t+h,om1_1+h*k3_1,om2_1+h*j2_1,th1_1,th2_1)
 
        om1_1 += k1_1/6 + k2_1/3 + k3_1/3 + k4_1/6
        om2_1 += j1_1/6 + j2_1/3 + j3_1/3 + j4_1/6
 
    h=h*2                                        #On retablit le pas
    err=array([abs(om1-om1_1),abs(om2-om2_1)])   #Calcule de l'erreur
    err=max(err)
 
    if err<=p:                          #Si l'erreur est acceptable
        T=append(T,t)                #Stockage des données
        Th1=append(Th1,th1)
        Th2=append(Th2,th2)
        Om1=append(Om1,om1)
        Om2=append(Om2,om2)
                           #Augmentation du pas 
        t+=h                        #et on passe a l'etape suivante
 
    else:                               #Sinon diminution du pas 
                               #Et on recommence la meme etape
 
        h = h *min( max( 0.9 * (p/err) **0.2,0.1),4.0) 
        if h > hmax:
            h = hmax
        if h < hmin:
            h=hmin

Mais le code n'est pas bon, il s'arrête rapidement est les valeurs sont très loin de celle obtenus sans faire varier h.
J'ai l'impression que le problème vient du fait que même avec un pas plus petit, l'erreur calculée ne diminue pas, elle augmente. Dès lors, on ne se trouve plus jamais dans le cas 'err<=p' et le programme ne s'arrête plus.
Est ce que quelqu'un sait pourquoi l'erreur augment alors que le pas diminue ?
Est-ce que vous pensez que le problème est purement mathématique ou qu'il vient du programme?

Si vous avez des remarques ou des commentaires sur le code, je prend!

Merci d'avance pour vos réponses.

[Nanzilla]

Penses à factoriser ton code tu fais un copier/coller de deux fois la même portion de code ce qui n'est pas bien ^^.
Pour ton problème, si tu passes par une erreur acceptable, il faut mettre à jour om1 et om2 pour qu'ils prennent om1_1 et om2_1 comme valeur sinon cela ne sert à rien. De même pour Om1 et Om2 qui doivent prendre les valeurs om1_1 et om2_1.
Après je ne connais pas la méthode, mais d'après le document pdf, le pas de temps est modifié aussi p66 Ligne 42 du document code: http://www.physique.usherbrooke.ca/p...HQ405_ipad.pdf

[aaaallleex]

Merci pour ta réponse.

Concernant le copier coller:
Les deux boucles sont presques les même mais elles prennent des valeurs différente à chaque pas, ce sont en fait deux boucle complètement différentes. Il est vrai que j'aurais du utiliser une seule routine dans une définition et l'appeler deux fois avec des valeurs différentes.

Pour la mise à jour des valeurs:
d'après ce que j'ai compris, la deuxième boucle ne sert que de vérification, et on utilise les valeurs de la première si l'erreur est acceptable. C'est pourquoi ce sont om1 et om2 qui sont ajoutés.


Par ailleurs, j'ai demander à mon prof d'info pourquoi l'erreur augmente alors que le pas diminue qui m'a répondu: que lorsqu'on utilises une routine à pas variable il faut faire attention au fait que lorsque le pas devient très petit, l'erreur numérique tend vers zéro (en théorie), par contre l'ereur que l'on fait (due à la précision de représentation) lors des calculs augmente... il faut donc trouver un optimum.

Un algorithme est proposé dans les numerical recipes dont on peux s'inspirer.

http://apps.nrbook.com/c/index.html, (aller aux pages pages 714-722).

Finalement mise en place d'une routine à pas variable est finalement un peu trop complexe pour mon niveau de programmation, initié depuis quelques semaines...

Je n'aurais pas le temps de corriger le programme dans les prochains jours, le sujet est donc résolue.