# “Singularity or stiff system suspected” and “Input value lies outside of range” errors when plotting complex pendulum

This is my first time posting a question so apologies for any issues.

I'm a first year physics student who just started using mathematica last week. Part of my physics course is to create a computer model of a complex pendulum. My pendulum is based off the angry arms exhibit from Questacon https://www.questacon.edu.au/visiting/galleries/wonderworks/exhibits/angry-arms

Essentially you have a solid T-Bar which rotates around the intersection. This t-bar then has 3 pendulums off the edge of each arm.

To model this, I defined the legrangian and Euler-legrangian of the system. I then used ndsolve to find the solutions to the differential equation, and used that solution to plot it.

Euler-Legrange Calculation:

(*The height and velocity of the bob of a pendulum, given then angle of the pendulum and the angle of the tbar*)
h[\[Theta]_, \[Theta]tbar_,t_] = dt*(1 - Cos[\[Theta]tbar[t]]) + dp*(1 - Cos[\[Theta][t]]);
v[\[Theta]_, \[Theta]tbar_,t_] = \[Theta]tbar'[t]*((dp*Sin[\[Theta][t]])^2 + (dt - dp*Cos[\[Theta][t]])^2)^(1/2) + \[Theta][t]*dp;
(*The total Kinetic Energy (U) and Potential Energy (V) of the system*)
U[\[Theta]1_, \[Theta]2_, \[Theta]3_, \[Theta]4_, t_] = 1/2 m*(v[\[Theta]1, \[Theta]4, t]^2 + v[\[Theta]2, \[Theta]4, t]^2 + v[\[Theta]3, \[Theta]4, t]^2);
V[\[Theta]1_, \[Theta]2_, \[Theta]3_, \[Theta]4_, t_] = m*g*(h[\[Theta]1, \[Theta]4, t] + h[\[Theta]2, \[Theta]4, t] + h[\[Theta]3, \[Theta]4, t]);
(*The calculation of the Legrangian and Euler-Legrangian*)
L[\[Theta]1_, \[Theta]2_, \[Theta]3_, \[Theta]4_, t_] = U[\[Theta]1, \[Theta]2, \[Theta]3, \[Theta]4, t] - V[\[Theta]1, \[Theta]2, \[Theta]3, \[Theta]4, t];
el = EulerEquations[L[\[Theta]l, \[Theta]r, \[Theta]m, \[Theta]t,t], {\[Theta]l[t], \[Theta]r[t], \[Theta]m[t], \[Theta]t[t]}, t];


Solving Euler-Legrange

Clear[\[Theta]l,\[Theta]r,\[Theta]m,\[Theta]t,dt,dp,g,m,tMax,nFrames]
pm = dt {Sin[\[Theta]t[t]], -Cos[\[Theta]t[t]]};
pl = dt {Sin[\[Theta]t[t] - N[Pi/2]], -Cos[\[Theta]t[t] - N[Pi/2]]};
pr = dt {Sin[\[Theta]t[t] + N[Pi/2]], -Cos[\[Theta]t[t] + N[Pi/2]]};
ppm = pm + dp {Sin[\[Theta]m[t]], -Cos[\[Theta]m[t]]};
ppr = pr + dp {Sin[\[Theta]r[t]], -Cos[\[Theta]r[t]]};
ppl = pl + dp {Sin[\[Theta]l[t]], -Cos[\[Theta]l[t]]};
dt = 1;
dp = 1/2;
g = 9.8;
m = 1;
tMax = 20;
initial = {\[Theta]l[0] == Pi/3, \[Theta]r[0] == 0, \[Theta]m[0] == 0, \[Theta]t[0] == Pi/2};
sol = First[NDSolve[Join[el,initial], {\[Theta]l[t], \[Theta]r[t], \[Theta]m[t], \[Theta]t[t]}, {t, 0, tMax}], Method -> Shooting]


Modelling

ParametricPlot[Evaluate[{ppl, ppr, ppm} /. sol],{t, 0, tMax},PlotStyle -> {Blue, Red, Green}, ImageSize -> Medium, PlotLegends -> {"Trajectory of leftmost pendulum", "Trajectory of rightmost pendulum", "Trajectory of middle pendulum"}]

nFrames = tMax*20;
frames = Table[Graphics[{
Gray, Thick, Line[{{0, 0}, pl, ppl}], Line[{{0, 0}, pm, ppm}],
Line[{{0, 0}, pr, ppr}],
Darker[Black], Disk[{0, 0}, 0.1],
Darker[Red], Disk[pr, 0.1], Disk[ppr, 0.1],
Darker[Blue], Disk[pl, 0.1], Disk[ppl, 0.1],
Darker[Green], Disk[pm, 0.1], Disk[ppm, 0.1]
} /. sol, PlotRange -> {{-2.5, 2.5}, {-2.5, 2.5}},
Background -> Lighter[Cyan]
],{t, tMax/nFrames, tMax, tMax/nFrames}];
ListAnimate[frames]


For some initial values (usually having everything but [\theta]t 0) the parametric plot and simulation both work perfectly and the only error I get is:

NDSolve::ntdvdae: Cannot solve to find an explicit formula for the derivatives. NDSolve will try solving the system as differential-algebraic equations.


For other initial values (Having more non zero initial values) the plot still runs alright but the pendulums aren't where I expect them and I get the additional error

NDSolve::ivres: NDSolve has computed initial values that give a zero residual for the differential-algebraic system, but some components are different from those specified. If you need them to be satisfied, giving initial conditions for all dependent variables and their derivatives is recommended


Finally for some seemingly random initial values, the parametric plot becomes a crazy ball, and the simulation freaks out. I also get:

NDSolve::ndsz: At t == 0.00453348, step size is effectively zero; singularity or stiff system suspected.
InterpolatingFunction::dmval: Input value {1/20} lies outside the range of data in the interpolating function. Extrapolation will be used.
InterpolatingFunction::dprec: The precision of input value {161/20} and/or the interpolation grid is insufficient to compute the value.


Note, the t value and input value number change depending on the initial values.

I've looked around for a solution but haven't been able to find anything. If you have any idea's let me know.

• Needs["VariationalMethods`"] is missing from code. – bbgodfrey Aug 12 '16 at 5:02
• Apologies, I had that at the top of the code, just forgot to copy paste it into the question – Patrick Armstrong Aug 13 '16 at 2:18