simulation of spherical pendulum

First, pleasue excuse my bad in English. I am from Germany and trying my best.

In the long run I want to simulate a spherical double pendulum. But it became obvious that I already have problems with simulating a single spherical pendulum.

Here is the code I use. The equations should ber right so far. It should run, when pasted into matlab:

%===================================================================

function Erg1 = solveitmass();

options = odeset('Mass' , @mass, 'RelTol',1e-8,'AbsTol',1e-8); 

Erg1 = ode45(@rightside, [0 10], [pi/4 0 0 0.01], options);

end

%============================Function for Massmatrix

function M  = mass( t, x)

m = 1;
l = 1;

M = [1 0  0    0;
     0 1  0    0;
     0 0  m*l^2  0;
     0 0  0  m*l^2*sin(x(1))^2]
end

%===========================Function for right hand side of equations

function dx = rightside(t, x)

l = 1;
g = 9.81;
m= 1;

dx = zeros(4, 1);

dx(1)=x(3)
dx(2)=x(4)
dx(3)=+m*l^2*sin(x(1))*cos(x(1))*x(4)^2-m*g*l*sin(x(1))
dx(4)= -2*m*l^2*sin(x(1))*cos(x(1))*x(3)*x(4)

end

%======================================================

The position of the mass is decribed by two angles. I will refer to them as azimuth-angle and elevation-angle. the first 2 positions of the initial condtion vector are the angles (Elevation first, second Azimuth). The last two positions describe the the angular rates. And here lies the problem. When starting with the angles with arbitrary values and the angular rates = 0, there ist no problem (with these values one does basicallyhave a 2D-Pendulum). But as soon as I give a value to the angular rate of the azimuth angle, the elevation angle doesn't cross zero anymore. I attached a plot, that shows the values of the elevation-angle.

http://www10.pic-upload.de/05.09.12/a746w5hckrl.jpg

Can anybody tell me what I am doing wrong? The equations should be right and as far as I can see, the implementation is also correct.

Edit

To me this doen't seem to be the natural behaviour of a pendulum. In my opinion the elevation angle should follow the course of a sinus wave. I can understand that the the angular rate of the azimuth angle gets bigger when the elevation gets smaller because the conversation of angular momentum, but I don't see a reason that the elevation cannot cross zero anymore, as soon as there is any movement in azimuth direction in the initial conditions. Instead the elevation angle really 'bounces' back again as soon as ist reaches zero. And I have no idea why this happens. I see no reason why the elevation angle just crosses zero. Since the Equations of movement are correct I was hoping that I am just using Matlab in some wrong way.

[myta]

I hope this can help http://www.youtube.com/watch?v=oPagRjbqbYA

[MARTÍN ZAPATA FERGUSON]

This might be a bit late of an answer, but what you're calculating is physically correct. Try to calculate the energy integral and see what happens when the angular momentum =/= 0. You'll find asymptotes at zero and pi, meaning your elevation can come arbitrarily close to those values but won't actually reach them. Essentially, the azimuthal spinning means you can't reach either pole.

Answer 1

[Bas Siebers]

Hi,

I have the same problem as you. Do you know a solution for this problem?