Discussion :

[Ginette123]

Bonjour,

Je débute dans la résolution numérique avec l’algorithme de Runge-Kutta. Voilà, j'aimerais simuler le mouvement d'un pendule avec les équations de Hamilton (avec les équations (172) du cette page ).

Comme je n'ai pas toujours accès au logiciel Matlab, j'utilise le logiciel Octave en ligne, avec le code suivant:

Code matlab :

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
 
function pendule()
% ---------- %
% Parametres %
% ---------- %
r = 5; % rayon du pendule 
m = 2; % masse du pendule
L = 6; % L = P_{\phi} = cst = moment cinetique
dt = 0.002;
duree = 7;
nbpts = duree/dt; % NomBre de PoinTS
% -------------- %
% Initialisation %
% -------------- %
tab = zeros(duree/dt,6); % [t,theta, phi, P_{\theta}, x, y]
tab(1,1) = 0;
tab(1,2) = 0.2*pi; % theta
tab(1,3) = 0; % dans le plan vertical passant pas l'axe des x
tab(1,4) = 0; % P_{\theta} = m*r^2*\dot{\theta} =0 : vitesse angulaire \dot{\theta}=0
tab(1,5) = r*sin(tab(1,2))*cos(tab(1,3)); % x
tab(1,6) = r*sin(tab(1,2))*sin(tab(1,3)); % y
% ---------------------------- %
% Resolution par Runge-Kutta 4 %
% ---------------------------- %
for k = 1:nbpts
    a1 = dt * f(tab(k,4), m, r);
    b1 = dt * g(tab(k,2), L, m, r);
    c1 = dt * h(tab(k,2), L, m, r);
    a2 = dt * f(tab(k,4) + a1/2, m, r);
    b2 = dt * g(tab(k,2) + b1/2, L, m, r);
    c2 = dt * h(tab(k,2) + c1/2, L, m, r);
    a3 = dt * f(tab(k,4) + a2/2, m, r);
    b3 = dt * g(tab(k,2) + b2/2, L, m, r);
    c3 = dt * h(tab(k,2) + c2/2, L, m, r);
    a4 = dt * f(tab(k,4) + a3, m, r);
    b4 = dt * g(tab(k,2) + b3, L, m, r);
    c4 = dt * h(tab(k,2) + c3, L, m, r);
    tab(k+1,1) = tab(k,1) + dt;
    tab(k+1,2) = tab(k,2) + (a1 + 2*a2 + 2*a3 + a4)/6; % theta
    tab(k+1,3) = tab(k,3) + (b1 + 2*b2 + 2*b3 + b4)/6; % phi
    tab(k+1,4) = tab(k,4) + (c1 + 2*c2 + 2*c3 + c4)/6; % P_{\theta}
    tab(k+1,5) = r*sin(tab(k+1,2))*cos(tab(k+1,3)); % x
    tab(k+1,6) = r*sin(tab(k+1,2))*sin(tab(k+1,3)); % y
end;
% --------- %
% Graphique %
% --------- %
plot(tab(:,5),tab(:,6),'r-');
% --------- %
% Fonctions %
% --------- %
function dottheta = f(pth,m,r)
dottheta = pth/(m*r*r);
 
function dotphi = g(theta,L,m,r)
dotphi = L/(m*r*r*sin(theta)*sin(theta));
 
function dotp = h(theta,L,m,r)
dotp = (L*L*cos(theta)/(m*r*r*sin(theta)*sin(theta)*sin(theta))) - (m*9.81*r*sin(theta));

D'après la théorie, je devrais obtenir une rosace. Or, la solution que j'obtiens est différente (voir figure obtenue). Il semble qu'il y a une répulsion par l'origine. Je ne comprends pas pourquoi. Qu'est-ce qu'il y a de faux ? Il y a déjà un moment que j'essaie de changer les paramètres, mais sinon, je ne vois pas comment résoudre ce problème. Lorsque je change les paramètres d'un tout petit peu, la solution devient très différente (par exemple avec dt = 0.002 et duree = 10). J'ai l'impression qu'il y a quelque chose de faux dans la programmation de la méthode, mais je ne sais pas où.

Merci beaucoup d'avance pour vos réponses.

Ginette

[dourouc05]

Ça ne résoudra pas ton problème, mais il serait nettement plus propre d'avoir une vraie fonction RK qui effectue juste l'intégration, en prenant en entrée les conditions initiales et une fonction dérivée (ce que fait déjà MATLAB, donc Octave : http://octave.sourceforge.net/odepkg...ion/ode45.html). Tu aurais donc d'un côté le système différentiel du pendule, de l'autre ton RK.

Ensuite, le pas peut avoir une très grande importance sur la méthode : elle peut très bien se mettre à converger ou à donner une solution complètement aberrante. Par exemple, http://codeflow.org/entries/2010/aug...s-runge-kutta/. Pour éviter ces problèmes, tu peux regarder du côté des méthodes à pas variable.

(Maintenant, vu la tête du graphe, j'aurais tendance à penser à une erreur de signe à un moment… Peut-être un qui se retrouve lors de l'évaluation d'une formule, uniquement à cause d'erreurs de précision.)

[Ginette123]

Salut,

Merci pour la méthode, c'est vrai que c'est mieux de commencer à voir ce que ça donne avec la fonction de Matlab/Octave dédiée à cela.

Ginette

[Ginette123]

Re-bonjour,

J'ai essayé avec une ode45 et ode23. Je trouve aussi qu'il y a comme une répulsion depuis le centre.
Les conditions initiales sont celles de Hairer, Norsett, Wanner à la page 32 et les équations celles de la page 33.

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
 
% f = f(dottheta, dotimpulsetheta, dotphi)
L=0.5;
f = @(t,x) [x(2); (L^2*cos(x(1))/((sin(x(1)))^3)) - sin(x(1)) ; L^2/((sin(x(1)))^2)];
options = odeset('MaxStep',1000,'RelTol',1e-5,'AbsTol',[1e-5 1e-5 1e-5],'InitialStep',100);
[t,xa] = ode23(f,[0 100],[1 0 0],options); % ode23 ou ode45 donnent le même résultat
plot(sin(xa(:,1)).*cos(xa(:,3)),sin(xa(:,1)).*sin(xa(:,3)));

Je pense que j'ai commis une erreur quelque part, mais où ...