# Coupled differential equations from Lagrangian — singularity being reached very quickly

I was wondering if there was a way to help the solutions to these differential equations not hit a singularity as quickly as they currently do. The equations are very complex (derived from a Lagrangian), so it could be that Mathematica is struggling because of that. Depending on the initial conditions, I'm only able to solve the equations for values of time (t) up to half a second or a few seconds if lucky. This is my input:

s = NDSolve[{θ''[t] + φ'[t]*Sin[φ[t]]*(Cos[θ[t]]*Cos[t] +
Sin[θ[t]]*Sin[t]) + Cos[φ[t]]*(θ'[t]*Sin[θ[t]]*Cos[t] +
Cos[θ[t]]*Sin[t] - θ'[t]*Cos[θ[t]]*Sin[t] - Sin[θ[t]]*Cos[t]) ==
(φ'[t]*Sin[φ[t]]*(Cos[θ[t]]*Cos[t] + Sin[θ[t]]*Sin[t]) - θ'[t]*
Cos[φ[t]]*(Cos[θ[t]]*Sin[t] - Sin[θ[t]]*Cos[t])) + (φ'[t]^2)*
Cos[θ[t]]*Sin[θ[t]] - 9.81*Sin[θ[t]], φ'[t]*Cos[φ[t]]*(Sin[θ[t]]*Cos[t] -
Cos[θ[t]]*Sin[t]) + Sin[φ[t]]*(θ'[t]*Cos[θ[t]]*Cos[t] -
Sin[θ[t]]*Sin[t] + θ'[t]*Sin[θ[t]]*Sin[t] - Cos[θ[t]]*Cos[t]) +
φ''[t]*((Sin[θ[t]])^2) + 2*θ'[t]*φ'[t]*Sin[θ[t]]*Cos[θ[t]] ==
((φ'[t]^2)*Cos[φ[t]]*(Sin[θ[t]]*Cos[t] - Cos[θ[t]]*Sin[t]) + θ'[t]*φ'[t]*
Sin[φ[t]]*(Cos[θ[t]]*Cos[t] + Sin[θ[t]]*Sin[t])),
θ[0] == Pi/8, φ[0] == Pi/8, θ'[0] == 0, φ'[0] == 0}, {θ, φ}, {t, 0.5}]

Plot[Evaluate[{θ[t], φ[t]} /. s], {t, 0, 0.5}, PlotStyle -> Automatic]


Sorry if this is a little messy, but I'm simply wanting to know whether Mathematica can do better than come up with a better numerical solution than one which is only reasonable for one second or so.

Whatever I do, it is possible that my equations could be wrong - they're derived from this Lagrangian:

l = φ'[t] Sin[φ[t]] (Sin[θ[t]] Cos[t] - Cos[θ[t]] Sin[t]) - θ'[ t] Cos[φ[t]]
(Cos[θ[t]] Cos[t] + Sin[θ[t]] Sin[t]) + 0.5*(θ'[t]^2 + φ'[t]^2 Sin[θ[t]]^2 + 1) +
9.81*Cos[θ[t]]

• Welcome to Mathematica.SE! I hope you will become a regular contributor. To get started, 1) take the introductory Tour now, 2) when you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge, 3) remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign, and 4) give help too, by answering questions in your areas of expertise. – bbgodfrey Sep 5 '15 at 0:49
• Is this what you see? !Mathematica graphics – Dr. belisarius Sep 5 '15 at 0:49
• If you do Simplify[{your system},0<t<Pi] you can make your system about 4x simpler. That can sometimes help with singularity problems or with understanding the system. But in this case it doesn't fix the problem at t==1.37. Even increasing WorkingPrecision->32 doesn't help. – Bill Sep 5 '15 at 1:04
• Yes belisarius, for those initial conditions. – turbodiesel4598 Sep 5 '15 at 12:04
• From the accepted answer it seems that the question was based on an erroneous equation - that's why I propose to close it. – Jens Sep 5 '15 at 20:43

Starting from your Lagrangian, the lazy way:

l = φ'[t] Sin[φ[t]] (Sin[θ[t]] Cos[t] - Cos[θ[t]] Sin[t]) - θ'[ t] Cos[φ[t]]
(Cos[θ[t]] Cos[t] + Sin[θ[t]] Sin[t]) + 0.5*(θ'[t]^2 + φ'[t]^2 Sin[θ[t]]^2 + 1) +
9.81*Cos[θ[t]]

Needs["VariationalMethods"]
eqs = EulerEquations[l, {θ[t], φ[t]}, t]
s = NDSolve[Join[eqs, {θ[0] == Pi/8, φ[0] ==  Pi/8, θ'[0] == 0,  φ'[0] == 0}],
{θ, φ}, {t, 0, 5}]

Plot[Evaluate[{θ[t], φ[t]} /. s], {t, 0, 5}]


• Wow! Now that is both awesome and looks totally reasonable - thanks! However, I am getting a bunch of errors when I try to use this myself. – turbodiesel4598 Sep 5 '15 at 15:51
• @turbodiesel4598 I have run belisarius' code and obtained the same answer he did. I also compared the equations from EulerEquations[l, {θ[t], φ[t]}, t] with those in your question and found that the second equations differed. – bbgodfrey Sep 5 '15 at 15:56
• Ah, no, it does work :) The line breaks caused Mathematica to read the input correctly. – turbodiesel4598 Sep 5 '15 at 16:17
• @turbodiesel4598 This is a common problem when copying code from StackExchange. See Meta1493 for a solution. – bbgodfrey Sep 5 '15 at 17:05
• @turbodiesel4598 This is what is happening mathematica.stackexchange.com/q/5563/193. BTW you need to include @Username to your comments for the other user to be alerted – Dr. belisarius Sep 5 '15 at 17:45

As noted in a comment by Bill, the calculation described in the question becomes stiff at

NDSolve::ndsz: At t == 1.3677154513713385, step size is effectively zero; singularity or stiff system suspected. >>


The natural course of action is to use Method -> "StiffnessSwitching" or Method -> "ImplicitRungeKutta", but doing so causes the calculation to fail instead at about t == 0.56.

This latter problem can be traced to the term φ''[t]*((Sin[θ[t]])^2), which vanishes at θ == 0, as it does in the second plot. It may be possible to handle this singularity by expanding the equations about this point. However, the OP first should check whether the equation and initial conditions are correct.

Why the first calculation, using the default Method, integrated past this point is uncertain. Perhaps, it was taking sufficiently large time steps there that it jumped over the singularity.

• I might try solving about that point then. Whatever I do, it is possible that my equations could be wrong - they're derived from this Lagrangian: L = φ'[t]sin[φ[t]](sin[θ[t]]cos[t] - cos[θ[t]]sin[t]) - θ'[t]cos[φ[t]](cos[θ[t]]cos[t] + sin[θ[t]]sin[t]) + 0.5*(θ'[t]² + φ'[t]²sin²[θ[t]] + 1) + 9.81*cos[θ[t]] And I may have made mistakes when finding the equations of motion. In fact, as of today, I have a good reason to think that there may be a mistake somewhere, as setting φ[0] = 0 causes φ to be totally stationary, although it does allow the equation to be solved without any singularities. – turbodiesel4598 Sep 5 '15 at 12:51
• @turbodiesel4598 If you derived the equations of motion from L using Mathematica, you might wish to put that derivation in your question. If not, try deriving them that way. Also, this problem looks like a pendulum of some sort. If so, it may have a constant of the motion, which would simplify the integration. The character of the numerical solution also makes me think that there may be a constant of the motion. – bbgodfrey Sep 5 '15 at 13:30