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 : Sélectionner tout - Visualiser dans une fenêtre à part
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.