3
$\begingroup$

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.

$\endgroup$
3
$\begingroup$

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]

enter image description here

enter image description here

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.

$\endgroup$
2
$\begingroup$

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]
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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