NDSolve coupled Non Linear differential equation, Obtain solution

Question

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

Answer 1

Increased WorkingPrecision can be achieved by solving Eqb for the second derivatives,

eqs = Flatten@Solve[Eqb, {D[Z[t], t, t], D[θ[t], t, t], D[ϕ[t], t, t]}] /. Rule -> Equal;
(* `eqs` is too long to be reproduced here *)
solb = NDSolve[{eqs, incb}, varb, {t, 0, 40}, WorkingPrecision -> 30, MaxSteps -> 50000];
Plot[Evaluate[{(θ[t] /. solb)}], {t, 0, 10}, PlotLegends -> {"θ"}, PlotRange -> All]

which produces a figure indistinguishable from that in the question. I conclude, therefore, that the computation is producing the proper result. Possibly, the question contains an incorrect set of ODEs or boundary conditions. Incidentally, Z[t] and ϕ[t] are approximately symmetric about zero.

Plot[Evaluate[{Z[t], θ[t], ϕ[t]} /. solb], {t, 0, 10}, PlotRange -> All]