2
$\begingroup$

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.

enter image description here

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]
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Szabolcs Mar 16 '16 at 18:36
  • $\begingroup$ @Szabolcs just explained question with full explanation. Thank you for nice suggestion. $\endgroup$ – Bikramjit Mar 16 '16 at 22:00
7
$\begingroup$

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]

enter image description here

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

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.