Coupled differential equations from Lagrangian — singularity being reached very quickly

Question

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

Answer 1

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

Answer 2

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.