Pendule sphérique - résolution numérique - Runge-Kutta 4

[triton32]

Bonjour,

J'essaie de résoudre numériquement le pendule sphérique, avec les équations de Hamilton, mais j'obtiens n'importe quoi.
Voici mon code en Matlab:

function pendulespherique()
% theta: nutation, phi:rotation
% ---------- %
% Parametres %
% ---------- %
r = 5; % rayon du pendule 
m = 2; % masse du pendule
L = 3; % L = P_{\phi} = cst = moment cinetique
dt = 0.01;
duree = 50;
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));

Je ne vois pas où sont mes erreurs: dans la méthode ou dans les paramètres ?

Je sollicite (un peu) votre aide pour éclairer ma lanterne

Merci


-----

[triton32]

Bonjour,

J'obtiens l'image suivante. Normalement, je devrait trouver des ellipses ouvertes, mais ce n'est pas droit ça ...

[triton32]

Une idée ?

[raph357]

Même pas une seule ?

[A voir en vidéo sur Futura]

[obi76]

Bonjour,

c'est quoi P dans votre cas ? Qu'il y ait un phi point et un theta point OK, mais P... ?

Et je ne vois pas bien ce que vous avez codé, c'est compliqué (et sans etre péjoratif, bordélique...).

Regardez en 2D comment vous faites, et ensuite faites la 3D en considérant que sur chaque direction c'est indépendant.