# Help me solve for the equations of motion for a double spherical pendulum

I've been working on a project to simulate the movement of a double spherical pendulum through Lagrangian mechanics. I found this link, which has the equations of motion in. I need to solve for the second time derivative of theta1, phi1, theta2, and phi2.

What I did was change all the time derivative symbols (') and replace them with a d meaning that theta1' is now theta1d etc. I'm certain that this is probably wrong but I'm not sure how I would do it otherwise

vars = {theta1dd, phi1dd, theta2dd, phi2dd}
Equations = {equations of movement}

Solve[equations, vars]


Is this incorrect syntax? If so, what should I do?

P.S. I'm only 16 so i'm sorry for my ignorance

• You are trying to solve a differential equation, i.e. an equation in which the unknown is a function, given some conditions on its values and the values of its derivatives. Start by looking up DSolve and perusing How To Solve a Differential Equation. May 13 '20 at 1:42
• Here is a simple example: DSolve[{y[x] == y'[x], y[0] == 1}, y[x], x]. Here we ask, what is the function $y(x)$ that is everywhere equal to its own first derivative $y'(x)$, and that has a value of 1 when its argument is 0? In this case, the answer is $y(x)=e^x$, as the above DSolve expression will tell you. It's a pretty deep rabbit hole... May 13 '20 at 1:46
• Here is a treatment of a spherical pendulum as part of an introductory college course on Mathematica: MATHEMATICA TUTORIAL for the Second Course. Part III: Spherical pendulum from Prof. Vladimir Dobrushkin at Brown. May 13 '20 at 1:55
• @MarcoB thanks soooo much. The paper from the professor is PERFECT, and I believe the final bit of code is just what I need to solve. :D thanks so much dude May 13 '20 at 2:55
• I’m very glad it helped! May 13 '20 at 3:48

## 1 Answer

I'm often a little bit unsure about how to go about entering everything I need into DSolve and NDSolve, so I usually like to start with the simplest example I can, and then slowly work my way up to what I actually want to do.

I would strongly recommend trying to work through these on your own as much as possible if you want to improve your understanding. But if you get stuck, I've added quite a bit of code here. I find this sort of simulation really interesting, so I couldn't help but look into it. There are some really nice answers for a 2D coupled pendulum on this question, so I hope that my answer here can help with the 3D case.

The VariationalMethods package has a nice function EulerEquations which automatically calculates the Euler-Lagrange equation for each variable and saves some extra work, so I'll be using it here.

Simple Pendulum:

Needs["VariationalMethods"]
x[t_] := Sin[θ[t]]
y[t_] := -Cos[θ[t]]
L = 1/2 m l^2 (x'[t]^2 + y'[t]^2) - m g l y[t] // FullSimplify
ee = EulerEquations[L, θ[t], t]


$$\frac{1}{2} l m \left(2 g \cos (\theta (t))+l \theta '(t)^2\right)$$

$$-l m \left(g \sin (\theta (t))+l \theta ''(t)\right)=0$$

Here I'm importing the VariationalMethods package, and then defining my Cartesian coordinates x[t] and y[t]. The Lagrangian is just the kinetic energy ($$1/2mv^2$$) minus the potential energy ($$mgy$$). Then, I ask EulerEquations to provide the Euler-Lagrange equations for the Lagrangian with respect to the coordinate $$\theta(t)$$ and independent variable $$t$$.

While I believe there is a closed form for the simple pendulum that relies on non-elementary functions, it's difficult to find analytical expressions for differential equations. Since the double spherical pendulum certainly won't have an analytical expression, I'll start using NDSolve here which provides a numerical result.

sol = First@NDSolve[{
ee /. {m -> 1, l -> 1, g -> 9.81},
θ'[0] == 0,
θ[0] == π/8
},
θ[t],
{t, 0, 20}
];


I'm replacing the mass $$m$$, the length $$l$$, and acceleration due to gravity $$g$$ in the Euler-Lagrange equation (using /.) before I ask it to solve the equation. There are a number of ways that you could specify these values, including just defining global variables m = 1; l = 1; g = 9.81 or making the functions accept these as arguments, but either way these should have numerical values by the time you call NDSolve.

Then I add in my initial conditions where I've set the angular velocity $$\theta'(0)$$ to 0, and the initial angle $$\theta(0)$$ to $$\pi/8$$. I'm asking it to solve for $$\theta(t)$$ for $$t$$ ranging from 0 to 20. It's unitless here, but if we assume $$m$$, $$l$$, and $$g$$ were in base SI units, we can read this as 0 seconds to 20 seconds.

Next, I want to plot this result to see what happened. I'm going to plot it 2 ways: first I'll plot $$\theta(t)$$ against $$t$$ to make sure it looks sinusoidal (I started with a small angle so it should be pretty close). Second, I want to see the motion of the pendulum.

Plot[
θ[t] /. sol,
{t, 0, 20},
AxesLabel -> {"t", "θ(t)"},
PlotRange -> {-π, π}
]
ParametricPlot[
{x[t], y[t]} /. sol,
{t, 0, 10},
AxesLabel -> {"x", "y"}
]


The second graph doesn't look all that interesting, but shows us the expected motion of a pendulum.

Spherical Pendulum:

I think I've explained most of the steps for the simple pendulum, so I'll include less explanation for these next cases.

Needs["VariationalMethods"]
x[t_] := Sin[θ[t]] Cos[ϕ[t]]
y[t_] := Sin[θ[t]] Sin[ϕ[t]]
z[t_] := -Cos[θ[t]]
L = m l^2 (x'[t]^2 + y'[t]^2 + z'[t]^2)/2 - m g l z[t] // FullSimplify;
ee = EulerEquations[L, {ϕ[t], θ[t]}, t];
sol = First@NDSolve[{
Splice[ee/.{m -> 1, l -> 1, g -> 9.81}],
ϕ'[0] == 0.5,
θ'[0] == 0,
ϕ[0] == 0,
θ[0] == π/8
},
{ϕ[t], θ[t]},
{t, 0, 100}
];
ParametricPlot3D[{x[t], y[t], z[t]} /. sol, {t, 0, 100}]


For a different set of starting conditions ($$\theta(0) = \pi/2$$ and only going up to a maximum time of 50), I get:

Double Spherical Pendulum:

Now that we understand a bit more about how NDSolve works and how to specify the arguments, we can try the toughest one. Notice that I defined the lengths l1 and l2 here. This helped me to keep the definitions of the Cartesian coordinates and the Lagrangian relatively short. This isn't my favourite way of doing it, but I haven't been able to figure out a good way to keep the definitions simple and not have the Cartesian coordinates include the lengths.

Needs["VariationalMethods"]
l1 = 1;
l2 = 1;
x1[t_] := l1 Sin[θ1[t]] Cos[ϕ1[t]]
y1[t_] := l1 Sin[θ1[t]] Sin[ϕ1[t]]
z1[t_] := -l1 Cos[θ1[t]]
x2[t_] := x1[t] + l2 Sin[θ2[t]] Cos[ϕ2[t]]
y2[t_] := y1[t] + l2 Sin[θ2[t]] Sin[ϕ2[t]]
z2[t_] := z1[t] - l2 Cos[θ2[t]]
L = m1 (x1'[t]^2 + y1'[t]^2 + z1'[t]^2)/2 +
m2 (x2'[t]^2 + y2'[t]^2 + z2'[t]^2)/2 - m1 g  z1[t] -
m2 g  z2[t] // FullSimplify;
ee = EulerEquations[
L, {ϕ1[t], θ1[t], ϕ2[t], θ2[t]}, t];
sol = First@NDSolve[{
Splice[ee /. {m1 -> 1, m2 -> 1, g -> 9.81}],
ϕ1'[0] == 0.75,
ϕ2'[0] == -0.215,
θ1'[0] == 0.2,
θ2'[0] == -0.09,
ϕ1[0] == 0.5,
ϕ2[0] == 0,
θ1[0] == 4 π/8,
θ2[0] == π/8
},
{ϕ1[t], θ1[t], ϕ2[t], θ2[t]},
{t, 0, 100},
Method -> {"EquationSimplification" -> "Residual"}
];
ParametricPlot3D[
Evaluate[{{x1[t], y1[t], z1[t]}, {x2[t], y2[t], z2[t]}} /. sol], {t,
0, 10}]


We can see the path of the first pendulum in blue, and the second in yellow.

Animation:

Because I couldn't stop myself, I decided to make an animation of what this might look like.

pendulum1[t_] := Evaluate[{x1[t], y1[t], z1[t]} /. sol]
pendulum2[t_] := Evaluate[{x2[t], y2[t], z2[t]} /. sol]
frames = Table[
Show[
ParametricPlot3D[
{x1[t], y1[t], z1[t]} /. sol,
{t, Max[0, time - 5], time},
ColorFunction -> (Directive[Red, Opacity[#4]] &)
],
ParametricPlot3D[
{x2[t], y2[t], z2[t]} /. sol,
{t, Max[0, time - 5], time},
ColorFunction -> (Directive[Blue, Opacity[#4]] &)
],
Graphics3D[{
Black,
Ball[{0, 0, 0}, 0.02],
Line[{{0, 0, 0}, pendulum1[time]}],
Line[{pendulum1[time], pendulum2[time]}],
Red,
Ball[pendulum1[time], 0.1],
Blue,
Ball[pendulum2[time], 0.1]
}
],
Axes -> True,
AxesOrigin -> {0, 0, 0},
Boxed -> False,
PlotRange -> {{-2, 2}, {-2, 2}, {-2, 2}},
ImageSize -> 500,
ViewAngle -> 17 Degree
],
{time, 0.01, 10, 0.05}
];
Export["~/Desktop/sphericalPendulum.gif", frames,
"DisplayDurations" -> 0.05]


(Actually, I had to decrease the resolution and number of frames in order to make the GIF small enough for upload.) Due to the "DisplayDurations" option, this should play at approximately real speed, i.e. 1 "unit" of time passes in the simulation for ever real second that passes.

EDIT:

It seems like I misunderstood the question in your post, sorry about that. Your method should work. I haven't tried it with the equations you found because I'm way too lazy to type out the one million characters necessary, but we can adapt some code I've already used. I switched the symbol names from $$\phi$$ and $$\theta$$ to phi and theta since you probably can't input symbols in Java. I'm also replacing all of the derivatives with your d/dd notation and removing any [t]s.

Needs["VariationalMethods"]
x1[t_] := l1 Sin[theta1[t]] Cos[phi1[t]]
y1[t_] := l1 Sin[theta1[t]] Sin[phi1[t]]
z1[t_] := -l1 Cos[theta1[t]]
x2[t_] := x1[t] + l2 Sin[theta2[t]] Cos[phi2[t]]
y2[t_] := y1[t] + l2 Sin[theta2[t]] Sin[phi2[t]]
z2[t_] := z1[t] - l2 Cos[theta2[t]]
L = m1 (x1'[t]^2 + y1'[t]^2 + z1'[t]^2)/2 +
m2 (x2'[t]^2 + y2'[t]^2 + z2'[t]^2)/2 - m1 g z1[t] - m2 g z2[t] //
FullSimplify;
ee = EulerEquations[L, {phi1[t], theta1[t], phi2[t], theta2[t]}, t] //
FullSimplify;
eqns = ee /. {
Derivative[1][theta1][t] -> theta1d,
Derivative[1][theta2][t] -> theta2d,
Derivative[1][phi1][t] -> phi1d,
Derivative[1][phi2][t] -> phi2d,
Derivative[2][theta1][t] -> theta1dd,
Derivative[2][theta2][t] -> theta2dd,
Derivative[2][phi1][t] -> phi1dd,
Derivative[2][phi2][t] -> phi2dd,
a_[t] :> a
};
Solve[eqns, {theta1dd, theta2dd, phi1dd, phi2dd}]


I'm afraid the output is long and ugly. I'm not sure if there's a simpler form. You could try another FullSimplify, but it would probably require you to manually rearrange things for it to be simpler. If it's possible, I would still recommend sticking to the Lagrangian method I show in my above examples, but if you can just copy and paste the functions, it might not be too much work to use the acceleration method. Since they're all elemental functions I think it will still run fairly quickly despite being so long.

• bruh, I can't believe u can do this type of stuff in mathematica, wtf. Since I'm programming in java, and I don't have access to the Euler-Lagrange equation solver, do you think there is anyway to slightly modify your code so that it could spit out an equation that directly represents the acceleration. this link has the equivalent equation for a 2D double pendulum. (it's second from the bottom). Do you know how I could find this otherwise? May 13 '20 at 8:41
• @Samuel, did you try inspecting the result of EulerEquations[L, {ϕ1[t], θ1[t], ϕ2[t], θ2[t]}, t]? MassDefect, this is a lovely tutorial, and as I told Samuel in another thread, this is how he should be attacking this: start simple, and slowly work yourself up to a more complicated problem. May 13 '20 at 9:52
• @Samuel, they are supposed to be, considering they're what MassDefect subsequently plugs into NDSolve[] to generate the animations. May 13 '20 at 14:19
• @SamuelCobb If you have access to Mathematica or the Wolfram Engine, the easiest way is to get it to calculate the equations for you. You can also find them with $\frac{d}{dt}\left(\frac{\partial L}{\partial q'}\right) - \frac{\partial L}{\partial q} = 0$ where q is each of the 4 coordinates (so you do it 4 times to get 4 equations equal to 0). EulerEquations just automates this process. Yeah, MMA is a powerful tool and I'm always amazed what you can do in a few dozen lines of code. I also think animations and plots are a great way to help convey physics–equations alone don't do it for me. May 13 '20 at 16:03
• @SamuelCobb I see, I think I misunderstood your question. You mentioned you wanted to use Lagrangian mechanics and I thought you meant you wanted to solve the pendulum using Lagrange mechanics, but you actually want the angular accelerations $\ddot{\phi_1}$, $\ddot{\phi_2}$, $\ddot{\theta_1}$, and $\ddot{\theta_2}$. You did say that, but I misunderstood. It can certainly be done, but there's a reason we usually use Lagrangian mechanics for this sort of thing. I checked the way you did it in your question (Solve[eqns, vars]), and it works fine. The only issue is that the output is quite long. May 14 '20 at 2:03