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) + 
  • 1
    $\begingroup$ 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. $\endgroup$
    – bbgodfrey
    Commented Sep 5, 2015 at 0:49
  • 1
    $\begingroup$ Is this what you see? !Mathematica graphics $\endgroup$ Commented Sep 5, 2015 at 0:49
  • 1
    $\begingroup$ 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. $\endgroup$
    – Bill
    Commented Sep 5, 2015 at 1:04
  • $\begingroup$ Yes belisarius, for those initial conditions. $\endgroup$ Commented Sep 5, 2015 at 12:04
  • $\begingroup$ From the accepted answer it seems that the question was based on an erroneous equation - that's why I propose to close it. $\endgroup$
    – Jens
    Commented Sep 5, 2015 at 20:43

2 Answers 2


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) + 

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}]

Mathematica graphics

  • $\begingroup$ 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. $\endgroup$ Commented Sep 5, 2015 at 15:51
  • $\begingroup$ @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. $\endgroup$
    – bbgodfrey
    Commented Sep 5, 2015 at 15:56
  • $\begingroup$ Ah, no, it does work :) The line breaks caused Mathematica to read the input correctly. $\endgroup$ Commented Sep 5, 2015 at 16:17
  • $\begingroup$ @turbodiesel4598 This is a common problem when copying code from StackExchange. See Meta1493 for a solution. $\endgroup$
    – bbgodfrey
    Commented Sep 5, 2015 at 17:05
  • 2
    $\begingroup$ @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 $\endgroup$ Commented Sep 5, 2015 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. >>

enter image description here

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.

enter image description here

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.

  • $\begingroup$ 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. $\endgroup$ Commented Sep 5, 2015 at 12:51
  • 1
    $\begingroup$ @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. $\endgroup$
    – bbgodfrey
    Commented Sep 5, 2015 at 13:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.