# NDSolve coupled Non Linear differential equation, Obtain solution

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?

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]

• Please edit the post and explain your question in clear English sentences. Don't just post some code. Posts on this site are required to be clear answerable questions. The clearer your phrasing, the more likely you will get a useful answer. – Szabolcs Mar 16 '16 at 18:36
• @Szabolcs just explained question with full explanation. Thank you for nice suggestion. – Bikramjit Mar 16 '16 at 22:00

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]


• Thank you for providing idea to solve equations with high precision, now I will try to change boundry conditions. – Bikramjit Mar 18 '16 at 3:43
• Dear @bbgodfrey can you please explain Flatten@Solve command, I shall be very thankful to you. – Bikramjit Mar 18 '16 at 17:30
• @Bikramjit Flatten@Solve{...] has the same effect as Flatten[Solve[...]]. Flatten removes an extra layer of curly brackets, { }. – bbgodfrey Mar 18 '16 at 23:14