# How to solve these ODEs using NDSolve?

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 == 2, ω2 == 3, ω3 == 4, ψ == 0, φ == 0, θ == Pi/6},
{ω1, ω2, ω3, ψ, φ, θ}, {t, 0, 120}]


I want to know how to avoid this error.

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 == 2, ω2 == 3, ω3 == 4, ψ == 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 to ψ == 10^-6.

And still another approach is to use the option,

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


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 == 2, ω2 == 3, ω3 ==
4, ψ == 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]