How to solve these ODEs using NDSolve?

Question

I have six odes and I cannot use DSolve. So I tried NDSolve. But it says there may be some errors.The code is such like this:

I1 = 2; I2 = 3; I3 = 4;
NDSolve[{I1*ω1'[t] + (I3 - I2)*ω2[t]*ω3[t] == 0, 
I2*ω2'[t] + (I1 - I3)*ω1[t]*ω3[t] == 0, 
I3*ω3'[t] + (I2 - I1)*ω2[t]*ω1[t] == 0, 
ω1[t] == φ'[t]*Sin[θ[t]]*Sin[ψ[t]] + θ'[t]*Cos[ψ[t]], 
ω2[t] == φ'[t]*Sin[θ[t]]*Cos[ψ[t]] - θ'[t]*Sin[ψ[t]], 
ω3[t] == φ'[t]*Cos[θ[t]] + ψ'[t], 
ω1[0] == 2, ω2[0] == 3, ω3[0] == 4, ψ[0] == 0, φ[0] == 0, θ[0] == Pi/6}, 
{ω1, ω2, ω3, ψ, φ, θ}, {t, 0, 120}]

I want to know how to avoid this error.

Answer 1

The computation fails at t = 0 because Sin[ψ[t]] vanishes there. However, this problem can be circumvented by eliminating Sin[ψ[t]] as follows. Define for convenience,

eq1 = ω1[t] - (φ'[t]*Sin[θ[t]]*Sin[ψ[t]] + θ'[t]*Cos[ψ[t]]);
eq2 = ω2[t] - (φ'[t]*Sin[θ[t]]*Cos[ψ[t]] - θ'[t]*Sin[ψ[t]]};

and construct the linear combinations,

eq1n = Simplify[eq1 Sin[ψ[t]] + eq2 Cos[ψ[t]]]
eq2n = Simplify[eq1 Cos[ψ[t]] - eq2 Sin[ψ[t]]]
(* Sin[ψ[t]] ω1[t] + Cos[ψ[t]] ω2[t] - Sin[θ[t]] φ'[t] *)
(* Cos[ψ[t]] ω1[t] - Sin[ψ[t]] ω2[t] - φ'[t] *)

(If this were not possible, the equations could not be solved, even in principle.)

Now replace eq1, eq2 by eq1n, eq2n.

I1 = 2; I2 = 3; I3 = 4;
s = NDSolveValue[{I1*ω1'[t] + (I3 - I2)*ω2[t]*ω3[t] == 0, 
    I2*ω2'[t] + (I1 - I3)*ω1[t]*ω3[t] == 0, 
    I3*ω3'[t] + (I2 - I1)*ω2[t]*ω1[t] == 0, 
    eq1n == 0, eq2n == 0, 
    ω3[t] == φ'[t]*Cos[θ[t]] + ψ'[t], 
    ω1[0] == 2, ω2[0] == 3, ω3[0] == 4, ψ[0] == 0, φ[0] == 0, θ[0] == Pi/6}, 
    {ω1[t], ω2[t], ω3[t], ψ[t], φ[t], θ[t]}, {t, 0, 120}];
Plot[Evaluate@s[[1 ;; 3]], {t, 0, 120}, ImageSize -> Large]
Plot[Evaluate@s[[4 ;; 6]], {t, 0, 120}, ImageSize -> Large]

Incidentally, the original equations also can be solved by modifying slightly the initial condition ψ[0] == 0 to ψ[0] == 10^-6.

And still another approach is to use the option,

Method -> {"EquationSimplification" -> "Residual"}

All give the same answer.

Answer 2

I suspect, the error message indicates that we are have 1/0 and I think this is because of Sin[] so I adopted the idea from here to change it in to Sinc.

Infinity::indet: Indeterminate expression ComplexInfinity+ComplexInfinity encountered.

I1 = 2; I2 = 3; I3 = 4;

Eq1 = I1*ω1'[t] + (I3 - I2)*ω2[t]*ω3[t] == 0;

Eq2 = I2*ω2'[t] + (I1 - I3)*ω1[t]*ω3[t] == 0;

Eq3 = I3*ω3'[t] + (I2 - I1)*ω2[t]*ω1[t] == 0;

Eq4 = ω1[t] == (φ'[t]*Sin[θ[t]]*Sin[ψ[t]] + θ'[t]*Cos[ψ[t]] ) /. {Sin[z_] :> 
    z*Sinc[z], Csc[z_] :> 1/(z*Sinc[z])};

Eq5 = ω2[t] == (φ'[t]*Sin[θ[t]]*Cos[ψ[t]] - θ'[t]*Sin[ψ[t]]) /. {Sin[z_] :> 
    z*Sinc[z], Csc[z_] :> 1/(z*Sinc[z])};

Eq6 = ω3[t] == φ'[t]*Cos[θ[t]] + ψ'[t];

sol = NDSolve[{Eq1, Eq2, Eq3, Eq4, Eq5, Eq6, ω1[0] == 2, ω2[0] == 3, ω3[0] == 
    4, ψ[0] == 0, φ[0] == 0, θ[0] == Pi/6}, {ω1, ω2, ω3, ψ, φ, θ}, {t, 0, 120}];

Plot[Evaluate[{ω1[t], ω2[t], ω3[t], ψ[t], φ[t], θ[t]} /. sol], {t, 0, 120}, PlotRange -> All]