solve Lagrange eqs.

I am trying to solve the Lagrange equations for a point mass rolling (withouth friction) in a spherical bawl.

Looking at the output it seems that the solution for theta "jumps" for multiples of pi

Any idea? The code follows and thanks for any help

Renzo

clear all
clf
syms theta(t) phi(t) 
g=9.8;
R=0.1; 
t0=0;
tf = 1; %seconds
dtheta=diff(theta,t);
dphi=diff(phi,t);
%Initial Conditions
c1=  pi/2+0.1;     % phi(0)
c2= 0;     % phidot(0)
c3=0;       % theta(0)
c4=0.001;     % thetadot(0) 
y0 = [c1 c2 c3 c4];
eq1 = diff(theta,2) == -2*dtheta*dphi*cos(phi)/sin(phi);
eq2 = diff(phi,2) == (dtheta^2)*sin(phi)*cos(phi)+g*sin(phi)/R;
vars = [theta(t); phi(t)];
[V,S] = odeToVectorField([eq1,eq2]);
M = matlabFunction(V,'vars', {'t','Y'});
interval = [t0 tf];  %time interval    
ySol = ode45(M,interval,y0);
tValues = linspace(interval(1),interval(2),100000);
%yValues = deval(ySol,tValues,1); %number 1 denotes first position likewise you can mention 2 to 6 for the next 5 positions
%plot(tValues,yValues)   %if you want to plot position vs time just swap here
theta=deval(ySol,tValues,3);
phi=deval(ySol,tValues,1);
dtheta=deval(ySol,tValues,4);
dphi=deval(ySol,tValues,2);
x=R.*sin(phi).*cos(theta);
y=-R.*sin(theta).*sin(phi);
z=R.*cos(phi);
%plot3(x,y,z)
% subplot(3,1,1)
% plot(tValues,x)
% ylabel('x')
% subplot(3,1,2)
% plot(tValues,y)
% ylabel('y')
% subplot(3,1,3)
% plot(tValues,z)
% ylabel('z')
% xlabel('t')
% sgtitle(['\phi_0= ',num2str(c1),'  d\phi_0= ',num2str(c2), ...
%     '  \theta_0= ',num2str(c3),'  d\theta_0= ',num2str(c4)]) 
% T=(dtheta*sin(phi)).^2+dphi.^2;
% U=1-cos(phi);
% E=T+U;
subplot(2,2,1)
plot(tValues,theta)
% set(gca,'YTick',-2*pi:pi/2:2*pi)
% ax.YTickLabel = {'-2 \pi','- \pi','0',' \pi','2 \pi'};
ylabel('\theta')
grid on
subplot(2,2,2)
plot(tValues,phi)
ylabel('\phi')
grid on
subplot(2,2,3)
plot(tValues,dtheta)
ylabel('d\theta/dt')
grid on
subplot(2,2,4)
plot(tValues,dphi)
ylabel('d\phi/dt')
grid on

●

[Torsten]

eq1 = diff(theta,2) == -2*dtheta*dphi*cos(phi)/sin(phi);

You divide by sin(phi) which gives a division by 0 when phi is a multiple of pi.

[Lorenzo]

Yes, Iwas aware of that and in fact the initial condition on phi was set different from zero but, silly of me, I did not think to phi=pi and multiples! I have to find a way to avoid these points

Thank you very mutch!!!

Answer 1

[David Goodmanson]

Hi Lorenzo,

The sin(phi) factor in the denominator does not really matter, in the sense that if phi can be zero, then the particle is moving in only a single plane, say for example the xz plane. That's the simple pendulum example and you only need one angle to describe it. So it's safe exclude that case and to assume that sin(phi) is not zero (see below).

Your initial conditions are such that the motion, when described by polar and azimuthal angles, is very extreme. I have not plotted the motion in 3d but if the path crosses at high speed and very close to the equilibrium point at the bottom, then theta will increment by pi very quickly, even if the motion itself is benign. These changes soften up the motion:

R = 1;
tf = 4; % seconds
c1 =  pi/4;     % phi(0)
c4 = 1;     % thetadot(0)

I also changed phi to be zero at the bottom of the sphere and pi at the top, which just involved a couple of minus signs. This seems the more natural choice since when the particle has small motions about the equilibrium point at the bottom of the sphere, then the cos factor is 1 instead of -1 which makes it easier to relate the motion to a circular pendulum with small oscillations.

Another code change speeds things up since with more rounded motion you don't need 100k time points. Those changes:

eq2 = diff(phi,2) == (dtheta^2)*sin(phi)*cos(phi) - g*sin(phi)/R;
tValues = linspace(interval(1),interval(2),10000);  
z = -R.*cos(phi)

One disadvantage of using eqn. 1 in your code is that some of the physics goes missing that way. eqn. 1 arises from the Lagrangian eqn

(d/dt) mr^2 sin(phi)^2 theta' = 0

meaning that

mr^2 sin(phi)^2 theta' = constant = A  % angular momentum about the z axis

If either phi or theta' is zero at t=0 then you are back to the xz plane pendulum. So assume that neither are zero. Then since A is conserved, neither sin(phi) or theta' can be zero at any time in the motion.

Solving (3) for theta' (dtheta in your code) and plugging it into eqn. 2 results in an equation in the phi variable only and introduces the centifugal potential, which repels the particle from the origin. That eqn can be solved in terms of two state space variables instead of 4, and (3) can then be used to obtain theta' and straight itegration gives theta. PCs are so fast that this way is not necessary now, but conceptually something is missed if it's not at least looked at.

I crossed myself up a few times since in physics it's a lot more common for theta and phi to be the other way around from the way it is in this question.