Trying to solve a coupled differential equations of motion. When there is no Working Precision mention it runs and give following plot for θ[t] but this is wrong, θ must be symmetric about 0.
later WorkingPrecision is set to 30
and now it gives me following error.
For the method NDSolve IDA, only machine real code is available. Unable to continue with complex values or beyond floating-point exceptions.
Are there any other options to generate solution?
Thanks in Advance
Code is as following.
(* External Loading *)
w[t_] := (20/1000)*Sin[Ω*t] + (20/1000)*Cos[Ω*t];
β[t_] := ((40/1000)*Sin[Ω*t] - (40/1000)*Cos[Ω*t])/(2*a);
Ω = 22;
(* Physical Parameters *)
m = 2350; ms = 7920; Js = 4495; l = 1.1; h = 1.78; hs = 19/10;
a = 24/20; c = 20000.;
k = 2200000; g = 9.81; di = 2.4; θc = ArcCos[l/di];
n = 3; d = 2258.78; b = 5324000; p = 3;
(* Set of Equations *)
Eqz = 2 k (Z[t]-w[t])+l m Cos[θ[t]+ϕ[t]] θ'[t]^2+2 l m Cos[θ[t]+ϕ[t]] θ'[t] ϕ'[t]+(-h m Cos[ϕ[t]]+l m Cos[θ[t]+ϕ[t]]+Cos[ϕ[t]] hs ms) ϕ'[t]^2+2 c (Z'[t]- w'[t])+(m+ms)Z''[t]+l m Sin[θ[t]+ϕ[t]] θ''[t]+(-h m Sin[ϕ[t]]+l m Sin[θ[t]+ϕ[t]] +Sin[ϕ[t]] hs ms) ϕ''[t];
Eqϕ = -g h m Sin[ϕ[t]] + g hs ms Sin[ϕ[t]]+g l m Sin[θ[t]+ϕ[t]]+2 a^2 k (ϕ[t] - β[t])+h l m Sin[θ[t]] θ'[t]^2+2 a^2 c (ϕ'[t]-β'[t])+2 h l m Sin[θ[t]] θ'[t] ϕ'[t]-h m Sin[ϕ[t]] Z''[t]+hs ms Sin[ϕ[t]] Z''[t]+l m Sin[θ[t]+ϕ[t]] Z''[t]+l^2 m θ''[t]-h l m Cos[θ[t]] θ''[t]+(Js+(h^2+l^2) m-2 h l m Cos[θ[t]]) ϕ''[t];
Eqθ = b (θ[t]/θc)^(2 n - 1) + d (θ[t]/θc)^(2 p) θ'[t]+l m (g Sin[θ[t]+ϕ[t]]-h Sin[θ[t]] ϕ'[t]^2+Sin[θ[t]+ϕ[t]] Z''[t]+l θ''[t]+(l - h Cos[θ[t]]) ϕ''[t]);
(* Formulation to Obtain Solution *)
Eqb = {Eqz == 0, Eqϕ == 0, Eqθ == 0};
Eqb = Chop[Eqb];
Eqb = Eqb /. {n1_Real -> Round[n1, 10^-11]};
incb = {Z[0] == 0, Z'[0] == 0, ϕ[0] == 0, ϕ'[0] == 0, θ[0] == 0,θ'[0] == 0};
varb = {Z, ϕ, θ};
solb = NDSolve[{ Eqb, incb}, varb, {t, 0, 40},Method -> {"EquationSimplification" -> "Residual"},WorkingPrecision -> 30, MaxSteps -> ∞]
Plot[Evaluate[{(θ[t] /. solb)}], {t, 0, 10},PlotLegends -> {"θ"}, PlotRange -> All]