3
$\begingroup$

Recently, someone showed me the following ODE system for a try in Mathematica's NDSolve,

eqns={z'[t]==-2 Sin[θ[t]] (1-(3 ξ1[t]^2)/1000),
θ'[t]==z[t] (12+(2 Cos[θ[t]] (1-(3 ξ1[t]^2)/1000))/Sqrt[1-z[t]^2]),
ξ1'[t]==3000 Cos[ϕ1[t]]-50 ξ1[t]-50 Sin[ϕ1[t]-ϕ2[t]] ξ2[t],ϕ1'[t]==300-1000 (1/250+(3 Cos[θ[t]] Sqrt[1-z[t]^2])/1000)-(3000 Sin[ϕ1[t]])/ξ1[t]-(50 Cos[ϕ1[t]-ϕ2[t]] ξ2[t])/ξ1[t],ξ2'[t]==50 Sin[ϕ1[t]-ϕ2[t]] ξ1[t]+50 ξ2[t],ϕ2'[t]==300-(50 Cos[ϕ1[t]-ϕ2[t]] ξ1[t])/ξ2[t]}


ibcs={z[0] == 7/10, θ[0] == 0, ξ1[0] == 107/10, ξ2[0] == 
  107/10, ϕ1[0] == π/2, ϕ2[0] == π/2}

But NDsolve's default settings:

sol = NDSolve[{eqns,ibcs}, {z, θ, ξ1, ξ2, ϕ1, ϕ2}, {t, Pi}]

does not solve this system.

The default NDSolve setting as above does not work properly(This is what causes my courisity): it is very demanding in CPU time and memory, and obtains poor/wrong numerical results in a very limited interval of $t$, e.g., $t\in[0,0.01]$ with warning message like: maximum steps reached at 0.01.

I have to try to specify the Method option of NDSolve in order to solve it easily. For example When using the following option, it works!:

Method->{"ExplicitRungeKutta","Coefficients" ->"EmbeddedExplicitRungeKuttaCoefficients","DifferenceOrder" -> 5,"StiffnessTest" -> False};

which gives: enter image description here

I also tried Matlab's ode45, and found the system could be easily solved by its default options settings:

function tryODE45
options = odeset('RelTol',1e-4,'AbsTol',1e-4*ones(1,6));
t0=[(7/10),0,(107/10),(107/10),(1/2)*pi,(1/2)*pi];
[T,Y] = ode45(@odeEqns,[0 .35],t0,options);
figure,
subplot(231),plot(T,Y(:,1),'r'),title('z');
subplot(232),plot(T,Y(:,2),'b'),title('\theta');
subplot(233),plot(T,Y(:,3),'m'),title('\xi_1');
subplot(234),plot(T,Y(:,4),'b'),title('\xi_2');
subplot(235),plot(T,Y(:,5),'g'),title('\phi_1');
subplot(236),plot(T,Y(:,6),'b'),title('\phi_2');

function dy=odeEqns(~,y)
dy = zeros(6,1);    % a column vector
dy(1)=-2*sin(y(2))*(1+(-3/1000)*y(3)^2);
dy(2)=y(1)*(12+2*cos(y(2))*(1+(-1)*y(1)^2)^(-1/2)*(1+(-3/1000)*y(3)^2));
dy(3)=(3000*cos(y(5))+(-50)*y(3)+(-50)*sin(y(5)+(-1)*y(6))*y(4));
dy(4)=(50*sin(y(5)+(-1)*y(6))*y(3)+50*y(4));
dy(5)=(300+(-1000)*((1/250)+(3/1000)*cos(y(2))*(1+(-1)*y(1)^2)^(1/2))+(-3000)*sin(y(5))*y(3)^(-1)+(-50)*cos(y(5)+(-1)*y(6))*y(3)^(-1)*y(4));
dy(6)=(300+(-50)*cos(y(5)+(-1)*y(6))*y(3)*y(4)^(-1));

enter image description here

My question is: why NDSolve has problem when handling such an ODE system? Is there any hint or tricks in order to properly use NDSolve in future cases?

$\endgroup$
  • $\begingroup$ Does NDSolve give an error? If so, what is it? $\endgroup$ – Michael E2 Aug 22 '16 at 10:18
  • $\begingroup$ Note that the code does not seem to have been copied correctly, the derivatives in particular. You may find this meta Q&A helpful. $\endgroup$ – Michael E2 Aug 22 '16 at 10:20
  • 2
    $\begingroup$ Does it work without "StiffnessTest" -> False? That feels like a hint to me. -- Why does the Mathematica code ask for a solution out to t == Pi, but the MATLAB code only goes up to 0.35? The default for Mathematica works fine if I set the time limit to 0.35. $\endgroup$ – Michael E2 Aug 22 '16 at 10:27
  • $\begingroup$ Thank you @MichaelE2. Because there are many Greek symbols, I use steampiano.net/msc to convert them. This might be the reason to the issues you mentioned. -- This is a very simple ODE problem if using Matlab or NDSolve with the method options I proposed: delivering results almost instantly. While the default NDSolve does not work and gives warning message like: maximum steps reached at 0.01 . $\endgroup$ – user6043040 Aug 22 '16 at 11:18
  • 2
    $\begingroup$ If I use your "ExplicitRungeKutta" method options, it fails sooner, at t == 0.426: i.stack.imgur.com/GueyO.png $\endgroup$ – Michael E2 Aug 22 '16 at 12:25
5
$\begingroup$

To clarify my comments, I thought I'd post the complete code I'm using. As of now, I don't see a problem. The main difference between the OP's results is that the OP's first code integrates over the interval {t, 0, Pi}, but in the plots and MATLAB code, the interval of integration only goes up to t == 0.35. That's a huge difference. The step size at t = 0.35 is already down to 2.90964*10^-6 due to the high frequency of oscillations of z[t]. I also find that neither the default setting not the Runge-Kutta option can integrate out to t == Pi; in fact, the RK method gives up sooner than the default (LSODA).

eqns = {(z')[t] == -2 Sin[θ[t]] (1 - (3 ξ1[t]^2)/1000),
   (θ')[t] == z[t] (12 + (2 Cos[θ[t]] (1 - (3 ξ1[t]^2)/1000))/Sqrt[1 - z[t]^2]),
   (ξ1')[t] == 3000 Cos[ϕ1[t]] - 50 ξ1[t] - 50 Sin[ϕ1[t] - ϕ2[t]] ξ2[t],
   (ϕ1')[t] == 300 - 1000 (1/250 + (3 Cos[θ[t]] Sqrt[1 - z[t]^2])/1000) -
     (3000 Sin[ϕ1[t]])/ξ1[t] - (50 Cos[ϕ1[t] - ϕ2[t]] ξ2[t])/ξ1[t],
   (ξ2')[t] == 50 Sin[ϕ1[t] - ϕ2[t]] ξ1[t] + 50 ξ2[t],
   (ϕ2')[t] == 300 - (50 Cos[ϕ1[t] - ϕ2[t]] ξ1[t])/ξ2[t]};

ibcs = {z[0] == 7/10, θ[0] == 0, ξ1[0] == 107/10, ξ2[0] == 107/10,
   ϕ1[0] == π/2, ϕ2[0] == π/2};

{sol} = NDSolve[{eqns, ibcs}, {z, θ, ξ1, ξ2, ϕ1, ϕ2},
    {t,(*Pi*)0.35}]; // AbsoluteTiming
ByteCount@sol
(*
  {0.041426, Null}
  1256632
*)

GraphicsGrid@Partition[
  Plot[#, {t, 0, 0.35}, PlotStyle -> Red] & /@ 
   Through[{z, θ, ξ1, ξ2, ϕ1, ϕ2}[t] /. sol],
  3]

Mathematica graphics

$\endgroup$
  • 2
    $\begingroup$ Ich spreche kein Matlab, but is OP's options = odeset('RelTol',1e-4,'AbsTol',1e-4*ones(1,6)); equivalent to AccuracyGoal->4, PrecisionGoal->4? If so, that speeds up the Mathematica code about three times, while still producing similar graphs. $\endgroup$ – Chris K Aug 22 '16 at 16:27
  • $\begingroup$ @ChrisK I don't know much Matlab, but that seems to be correct, although it seems to use them in a different way. In Matlab, the condition is error < max(RelTol*abs(y(i)), AbsTol(i)) and in Mathematica it's equivalent to the sum instead of max. $\endgroup$ – Michael E2 Aug 22 '16 at 16:38
  • $\begingroup$ Thank you! It seems to be this reason. According to the documents of NDSolve, there are also similar implementation as matlab's ode45 in NDSolve. $\endgroup$ – user6043040 Aug 23 '16 at 2:41
  • $\begingroup$ @user6043040 You're welcome! $\endgroup$ – Michael E2 Aug 23 '16 at 2:47

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.