Mathematica Stack Exchange is a question and answer site for users of Mathematica. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

As mentioned in the documentation for NDSolve it's often convenient to reduce a 2nd order ODE to a system of first order equations. When I do this however I seem to see a significant speed reduction in finding the solution.

Is there any reason this would be the case?

An Example:

Some definitions:

rstar[r_] := r + 2 M Log[r/(2 M) - 1];
rH = 200001/100000;
r0 = 10;
\[Lambda][l_] = l (l + 1);

Take the equation

 eq[\[Omega]_,l_] := \[CapitalPhi]''[r] + (2 (r - M))/(
 r (r - 2 M)) \[CapitalPhi]'[
  r] + ((\[Omega]^2 r^2)/(r - 2 M)^2 - \[Lambda][l]/(
   r (r - 2 M))) \[CapitalPhi][r] == 0;

Solve it with certain ICs:

init=-0.0000894423075560122420468703835499 + 
0.0000447222944185058822813688948339 I;
dinit=-4.464175354293244250869336196691640386266791`30.*^-6 - 
 8.950483248390306670770345406047835993931665`30.*^-6 I;

 sol := \[CapitalPhi] /. 
 Block[{$MaxExtraPrecision = 100}, 
 NDSolve[{eq[1/10, 1], \[CapitalPhi][rinf] == 
    init, \[CapitalPhi]'[rinf] == dinit}, \[CapitalPhi], {r, r0, 
   rinf}, WorkingPrecision -> wp, AccuracyGoal -> ac, 
  MaxSteps -> \[Infinity]]][[1]];

Now as set of first order equations

It turns out this system can be written as a first order set in terms of a related dependent variable $r^*$, and an effective potential $V$. Some more definitions:

init2=-0.8944230755601224204687038354990773373534 + 
0.4472229441850588228136889483392836606307 I;
dinit2=-0.04472224961131835705979008430399621833410 - 
0.08944221816744666391325700074861130268693 I;

r[rs_] := 2 (M + M ProductLog[E^(-1 + rs/(2 M))]);
V[rs_, \[Omega]_,l_] := \[Omega]^2 - (1 - (2 M)/r[rs]) (\[Lambda] [l]/(r[rs])^2 + (2 M)/(r[rs])^3);
rsH = N[rstar[rH], wp];
rsinf = N[rstar[rinf], wp];
rs0 = N[rstar[r0], wp];

Solve the first order system

  sol2 :=  {R, Rp} /. 
  NDSolve[{Rp[rs] == R'[rs], Rp'[rs] == -V[rs, 1/10, 1] R[rs], 
  R[rsinf] == init2, Rp[rsinf] == dinit2}, {R, Rp}, {rs, rsinf, 
  rs0}, WorkingPrecision -> wp, AccuracyGoal -> ac, 
  MaxSteps -> \[Infinity]][[1]];

Now run:

 sol // AbsoluteTiming
 sol2 // AbsoluteTiming

I find that sol2 takes roughly three times as long, and for even large $\omega,\ell$ this time difference gets more pronounced

share|improve this question
an example would be helpful – acl Oct 15 '12 at 11:18
I'm interested in this, too. I thought Mathematica automatically decomposed the 2nd order ODE into a set of 1st order ODEs, and then those are solved. Maybe reducing the order by hand, taking into account certain simplifications that only you know, gives something faster than what Mathematica can deduce. On second read: Or maybe Mathematica finds a reduction that is better than what you have accomplished. – Eric Brown Oct 15 '12 at 15:13
I've added my working example (let me know if anything is missing; times I get for second order ~7s and ~21s for the first order reduction. I should note that my first order set is found by introducing a new dependent var $r^*$ related to the old $r$, and the $R:=r \Phi$ is the new independent var. These ICs are only good for $(\omega,\ell)=(1/10,1)$. The problem gets much much worse when $(\omega,\ell)$ get larger. – fpghost Oct 15 '12 at 18:35

As acl points out, it would be nice to have an example illustrating when this happens. Here's an example (a basic pendulum equation) illustrating that it need not happen.

x2[t_] = x[t] /. First[
     NDSolve[{x''[t] == 2 Sin[x[t]] - 2 x[t], x[0] == 1, x'[0] == 0},
      x[t], {t, 0, 20}]]; // AbsoluteTiming
{x1[t_], y1[t_]} = {x[t], y[t]} /. First[
     NDSolve[{x'[t] == y[t],
       y'[t] == 2 Sin[x[t]] - 2 x[t], x[0] == 1, y[0] == 0},
      {x[t], y[t]}, {t, 0, 20}]]; // AbsoluteTiming
ParametricPlot[{{x1[t], y1[t]}, {x2[t], x2'[t]}},
 {t, 0, 20}, PlotStyle -> {
   Directive[{Thickness[0.02], Blue}],
   Directive[{Thickness[0.005], Yellow}]}]

enter image description here

share|improve this answer
Nice example that! However, if we use the method of lines option of NDSolve, doesn't it "break" the PDE into a set of ODEs? – drN Oct 15 '12 at 14:57
Hi drN, this is an ODE. PDE would have something like: eqns,{var1,var2,var3}, {x,0,L}, {t,0,T} :-) – Eric Brown Oct 15 '12 at 15:05
@EricBrown Gotcha! But isn't it the nature of MethodOfLines to generate a set of equations from a single ODE or PDE? – drN Oct 15 '12 at 20:26
Yipper. Somewhere in the Mathematica documentation (I think Advanced NDSolve) there is a nice graphic for this. Method of Lines is for reducing PDEs -> set of ODES. But the resulting ODEs (and the ODE in the Mark's comment above) may be greater than 1st order, and then those are reduced to first order by techniques like substitution which result in as many 1st order equations as the order of the original nth order ODE. – Eric Brown Oct 16 '12 at 3:37

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

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