# Why NDSolve's default method fails in solving this ODE system?

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:

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));


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?

• Does NDSolve give an error? If so, what is it? – Michael E2 Aug 22 '16 at 10:18
• Note that the code does not seem to have been copied correctly, the derivatives in particular. You may find this meta Q&A helpful. – Michael E2 Aug 22 '16 at 10:20
• 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. – Michael E2 Aug 22 '16 at 10:27
• 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 . – user6043040 Aug 22 '16 at 11:18
• If I use your "ExplicitRungeKutta" method options, it fails sooner, at t == 0.426: i.stack.imgur.com/GueyO.png – Michael E2 Aug 22 '16 at 12:25

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]


• 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. – Chris K Aug 22 '16 at 16:27
• @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. – Michael E2 Aug 22 '16 at 16:38
• 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. – user6043040 Aug 23 '16 at 2:41
• @user6043040 You're welcome! – Michael E2 Aug 23 '16 at 2:47