3
$\begingroup$

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];
M=1;
rinf=10000;
rH = 200001/100000;
r0 = 10;
wp=30;
ac=wp-8;
\[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

Improve this question
$\endgroup$
3
  • 3
    $\begingroup$ an example would be helpful $\endgroup$
    – acl
    Commented Oct 15, 2012 at 11:18
  • $\begingroup$ 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. $\endgroup$
    – Eric Brown
    Commented Oct 15, 2012 at 15:13
  • $\begingroup$ 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. $\endgroup$
    – fpghost
    Commented Oct 15, 2012 at 18:35

1 Answer 1

Reset to default
2
$\begingroup$

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

Improve this answer
$\endgroup$
4
  • $\begingroup$ 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? $\endgroup$
    – dearN
    Commented Oct 15, 2012 at 14:57
  • $\begingroup$ Hi drN, this is an ODE. PDE would have something like: eqns,{var1,var2,var3}, {x,0,L}, {t,0,T} :-) $\endgroup$
    – Eric Brown
    Commented Oct 15, 2012 at 15:05
  • $\begingroup$ @EricBrown Gotcha! But isn't it the nature of MethodOfLines to generate a set of equations from a single ODE or PDE? $\endgroup$
    – dearN
    Commented Oct 15, 2012 at 20:26
  • $\begingroup$ 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. $\endgroup$
    – Eric Brown
    Commented Oct 16, 2012 at 3:37

Your Answer

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

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