# Numerical solution for coupled PDEs

I am a trying to solve a physics problem, one that resembles the double pendulum. In my case it's a spring pendulum combination. I was able to write the lagrangian for the system and use the Euler Lagrange equations, to find the differential equations, which are the eq. of motion. The set up is the following:

I want to emphasize, the what I need help with is, how does mathematica solve numerically coupled PDEs. I am writing down, what I have done up until now:

x1[t_] = lf[t]*Sin[Phi1[t]]
y1[t_] = lf[t]*(-Cos[Phi1[t]])
x2[t_] = lf[t]*Sin[Phi1[t]] + l2*Sin[Phi2[t]]
y2[t_] = lf[t]*(-Cos[Phi1[t]]) - l2*Cos[Phi2[t]]

r1 = {x1[t], y1[t]}
r2 = {x2[t], y2[t]}

v1 = D[r1, t]
v2 = D[r2, t]

T = (1/2)*(m1*v1 . v1 + m2*v2 . v2);
U = m1*g*lf[t]*(1 - Cos[Phi1[t]]) + (1/2)*k1*(lf[t] - l1)^2 +
m2*g*(lf[t]*(1 - Cos[Phi1[t]]) + l2*(1 - Phi2[t]))

L = Simplify[T - U]

dLdPhi1 = D[L, Phi1[t]] // Simplify;
dLdPhiDot1 = D[L, Phi1'[t]] // Simplify;
ELEQ1 = dLdPhi1 == D[dLdPhiDot1, t];
ELEQ1 = Simplify[ELEQ1, {l1 != 0, m1 != 0, m2 != 0, l2 != 0, k1 != 0}]

dLdPhi2 = D[L, Phi2[t]] // Simplify;
dLdPhiDot2 = D[L, Phi2'[t]] // Simplify;
ELEQ2 = dLdPhi2 == D[dLdPhiDot2, t];
ELEQ2 = Simplify[ELEQ2, {l1 != 0, m1 != 0, m2 != 0, l2 != 0, k1 != 0}]

dLlf = D[L, lf[t]] // Simplify;
dLdlfDot = D[L, lf'[t]] // Simplify;
ELEQ3 = dLlf == D[dLdlfDot, t];
ELEQ3 = Simplify[ELEQ3, {l1 != 0, m1 != 0, m2 != 0, l2 != 0, k1 != 0}]

As one can see, we have three PDEs. Now I am trying to solve these equations, numerically, but I don't know how mathematica solves coupled PDEs numerically.

There are 4 degrees of freedom, so 4 equations of motions. In this I had 2 springs, not one like you. But you can easily change that if you want.

how does mathematica solve numerically coupled PDEs

Use NDSolve? Is this what you mean by the above question? This is how to do it for the above set up. NDSolve takes care of everything. Coupled or not coupled.

Assuming both springs have the same relaxed length of $$L$$. Starting by finding the Lagrangian $$\mathcal{L}=T-V$$. For $$m_{1}$$ \begin{align*} T_{1} & =\frac{1}{2}m_{1}\left( \dot{x}_{1}^{2}+\left( \left( L+x_{1}\right) \dot{\theta}_{1}\right) ^{2}\right) \\ V_{1} & =-m_{1}g\left( L+x_{1}\right) \cos\theta_{1}+\frac{1}{2}k_{1}% x_{1}^{2} \end{align*} And for $$m_{2}$$ \begin{align*} T_{2} & =\frac{1}{2}m_{2}\left( \left( \dot{x}_{2}+\dot{x}_{1}\cos\left( \theta_{1}-\theta_{2}\right) \right) ^{2}+\left( \dot{x}_{1}\sin\left( \theta_{1}-\theta_{2}\right) \right) ^{2}\right) \\ & +\frac{1}{2}m_{2}\left( \left( \left( L+x_{2}\right) \dot{\theta} _{2}+\left( L+x_{1}\right) \dot{\theta}_{1}\cos\left( \theta_{1}-\theta _{2}\right) \right) ^{2}+\left( \left( L+x_{1}\right) \dot{\theta} _{1}\sin\left( \theta_{1}-\theta_{2}\right) \right) ^{2}\right) \\ V_{2} & =-m_{2}g\left( \left( L+x_{1}\right) \cos\theta_{1}+\left( L+x_{2}\right) \cos\theta_{2}\right) +\frac{1}{2}k_{2}x_{2}^{2}% \end{align*} Hence $$\mathcal{L}=\left( T_{1}+T_{2}\right) -\left( V_{1}+V_{2}\right)$$ There are 4 generalized coordinates, $$x_{1},x_{2},\theta_{1},\theta_{2}$$. Now Mathematica is used to obtain the four equations of motion to help with the algebra. Once $$x_{1},x_{2},\theta_{1},\theta_{2}$$ are solved for, the position of each mass $$m_{1},m_{2}$$ is fully known at each time instance, and each mass motion can be animated. The four equations of motion are

\begin{align*} \frac{d}{dt}\left( \frac{\partial \mathcal{L}}{\partial\dot{x}_{1}}\right) -\frac{\partial \mathcal{L}}{\partial x_{1}} & =0\\ \frac{d}{dt}\left( \frac{\partial \mathcal{L}}{\partial\dot{x}_{2}}\right) -\frac{\partial \mathcal{L}}{\partial x_{2}} & =0\\ \frac{d}{dt}\left( \frac{\partial \mathcal{L}}{\partial\dot{\theta}_{1}}\right) -\frac{\partial \mathcal{L}}{\partial\theta_{1}} & =0\\ \frac{d}{dt}\left( \frac{\partial% \mathcal{L}}{\partial\dot{\theta}_{2}}\right) -\frac{\partial \mathcal{L}}{\partial\theta_{2}} & =0 \end{align*}

The rest is done using Mathematica to help with the algebra

ClearAll[x1, x2, θ1, θ2, t]
T1 = 1/2 m1 (x1'[t]^2 + ((L + x1[t]) θ1'[t])^2);
V1 = -m1 g (L + x1[t]) Cos[θ1[t]] + 1/2 k1 x1[t]^2;
T2 = 1/2 m2 ((x2'[t] +
x1'[t] Cos[θ1[t] - θ2[t]])^2 + (x1'[
t] Sin[θ1[t] - θ2[t]])^2) +
1/2 m2 (((L + x2[t]) θ2'[
t] + ((L + x1[t]) θ1'[
t] Cos[θ1[t] - θ2[t]]))^2 + ((L +
x1[t]) θ1'[t] Sin[θ1[t] - θ2[t]])^2);
V2 = -m2 g ((L + x1[t]) Cos[θ1[t]] + (L + x2[t]) Cos[θ2[
t]]) + 1/2 k2 x2[t]^2;

Set up Lagrangian

(lag = (T1 + T2) - (V1 + V2)) // Simplify

Equation for x1

(eq1 = D[D[lag, x1'[t]], t] - D[lag, x1[t]] == 0) // Simplify

Equation for x2

(eq2 = D[D[lag, x2'[t]], t] - D[lag, x2[t]] == 0) // Simplify

Equation for θ1

(eq3 = D[D[lag, θ1'[t]], t] - D[lag, θ1[t]] == 0) // Simplify

Equation for θ2

(eq4 = D[D[lag, θ2'[t]], t] - D[lag, θ2[t]] ==     0) // Simplify

Numerically solve the equations of motion

pars = {L -> 1, m1 -> 1, m2 -> 2, g -> 9.8, k1 -> 10, k2 -> 30};
ic = {θ1[0] == 5 Degree  , θ1'[0] == 0, θ2[0] ==
3 Degree, θ2'[0] == 0, x1[0] == 0, x1'[0] == 0, x2[0] == 0,
x2'[0] == 0};
eqs = Flatten[{eq1, eq2, eq3, eq4}] /. pars
numericalSolution=First@NDSolve[{eqs,ic},{x1,x2,θ1,θ2},{t,0,20}];

Plot solution (angles)

Plot[Evaluate[({θ1[t], θ2[t]} /. numericalSolution)*180/
Pi], {t, 0, 20}, PlotRange -> All,
AxesLabel -> {"time", "Angle (in degree)"}, ImageSize -> 400,
PlotLegends -> {"mass 1", "mass 2"}]

PLot solution for x1,x2

Plot[Evaluate[{x1[t], x2[t]} /. numericalSolution], {t, 0, 20},
PlotRange -> All,
AxesLabel -> {"time", "spring extensions in meters"},
ImageSize -> 400, PlotLegends -> {"mass 1", "mass 2"}]

You can do similar thing for your ode's.

I wanted to make a Manipulate of the above, but I did not have the time.

• But if I am given the coordinates for the mass m1 as: $\vec r_1(t)=l_f(t)[sin(\phi_1(t)),-cos(\phi_1(t))]$ and for m2: $\vec r_2(t)=[l_f(t)sin(\phi_1(t))+l_2sin(\phi_2(t)),-l_f(t)cos(\phi_1(t))-l_2cos(\phi_2(t))]$ then I do have only 3 eq. right? for phi1,phi2 and l_f(t). One more thing, relaxed length of a spring, is meant it's length with no weight on it right? Because if we have the weight, then we have the equlibrium length Commented Jan 15, 2023 at 20:22
• is meant it's length with no weight on it right? Yes. The relaxed length means the length of the spring when no extensive or compressive forces are acting on it. So the x will include any extension due to adding the mass and any pull from the second mass. Commented Jan 15, 2023 at 20:47
• then I do have only 3 eq. right? You should have as many equations as the number of degrees of freedom. Right? Clearly there are 4 DOF as you have no constrains anywhere. each mass can move along its x and also can rotate at same time. So 4 is total. This is in the setup I showed. I have 2 springs but you have only 1 spring. So, for your setup, yes, you will only have 3 DOF, so only 3 ODE's. You can use NDSolve to solve these. It is the same process, but you will have one less equation (no x2 change in your case). My setup is little more general. I used this setup since I had it before. Commented Jan 15, 2023 at 20:48
• btw, the whole point of this answer is just to show using NDSolve. Which is what you have asked about. But I did not just want to just say that and stop, so I used a setup I had similar to your to illustrate this to make it more complete. Commented Jan 15, 2023 at 21:01
• Thanks for the help, really appreciate it Commented Jan 15, 2023 at 21:02