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

1 Answer 1

8
$\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$
3
  • $\begingroup$ Thank you for providing idea to solve equations with high precision, now I will try to change boundry conditions. $\endgroup$
    – Bikramjit
    Mar 18, 2016 at 3:43
  • $\begingroup$ Dear @bbgodfrey can you please explain Flatten@Solve command, I shall be very thankful to you. $\endgroup$
    – Bikramjit
    Mar 18, 2016 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, 2016 at 23:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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