I need to solve a second order ODE numerically. The ODE depends on two parameters (a,b). Things work fine when 'a' is small, but for large 'a' the solutions are oscillating rapidly and Mathematica takes a long time to solve or eventually just runs out of memory.

I need to integrate over quite a large range (10,000) and that is part of the problem, but I actually only need the value of the InterpolatingFunction produced at the end point. Is there a way to tell Mathematica I just want this last point? And not store the rest (very large) InterpolatingFunction in memory? i.e. just integrate so far, take that point use as ICs for next leg, then integrate to next pt, take that as IC and integrate onwards,etc...

Or just some other strategy for using NDSolve with such highly oscillatory solutions.

Some definitions:

rstar[r_] := r + 2 M Log[r/(2 M) - 1];
$MinPrecision = 45;
wp = $MinPrecision;
ac = $MinPrecision - 8;
λ[l_] = l (l + 1);
rinf = 10000;
rH = 200001/100000;
nH = 200;

The ODE is:

eq[ω_, l_] := Φ''[r] + (2 (r - M))/(
 r (r - 2 M)) Φ'[
 r] + ((ω^2 r^2)/(r - 2 M)^2 - λ[l]/(
 r (r - 2 M))) Φ[r] == 0

Without going into detail about why I have these initial conditions, they are :

HorizonICs[l_?IntegerQ, ω_?NumericQ] := 
 Module[{ΦinrH, dΦinrH},
   b[0] = 1; b[-1] = 0; b[-2] = 0;
   b[n_] := 
  b[n] = Simplify[
 1/(2 n (n - 4 I ω)) (-2 I ω b[-2 + n] + 
    2 I n ω b[-2 + n] + l b[-1 + n] + l^2 b[-1 + n] + 
    n b[-1 + n] - n^2 b[-1 + n] - 4 I ω b[-1 + n] + 
    8 I n ω b[-1 + n])];
  uintrunc[r_, n_] := Sum[b[i] (r - 2 M)^i, {i, 0, n}];
  ΦinrH = (Exp[-I ω rstar[rH]] uintrunc[rH, nH])/(
  2 M);
 dΦinrH = 
  D[(Exp[-I ω rstar[r]] uintrunc[r, nH])/(2 M), r] /. 
   r -> rH;

Solve as

ΦinExt[ωω_ , 
 l_] := ΦinExt[ωω, l] = Φ /. 
 NDSolve[{eq[ωω, l], Φ[rH] == 
    N[HorizonICs[l, ωω][[1]], wp], Φ'[
     rH] == N[HorizonICs[l, ωω][[2]], 
     wp]}, Φ, {r, rH, rinf}, WorkingPrecision -> wp, 
  AccuracyGoal -> ac, MaxSteps -> ∞, 
  Method -> "StiffnessSwitching"][[1]];

You can also look for another solution that has the property of being simple near infinity and we set the ICs there. This is the particular solution that seems to be really really slow and causes memory crash.

Some definitions:

M = 1;
r0 = 5/2;
rstar[r_] := r + 2 M Log[r/(2 M) - 1];
$MinPrecision = 45;
wp = $MinPrecision;
ac = $MinPrecision - 8;
λ[l_] = l (l + 1);
rinf = 10000;
ninfphase = 

Set init conditions:

Infinitycs = Module[{n = ninfphase, c}, 
veqexp = 
 Series[(-2 - l r - l^2 r + 2 (r + I r^3 ω) 
\!\(\*SuperscriptBox["v", "′",
MultilineFunction->None]\)[r] + (-2 + r) r^2 
\!\(\*SuperscriptBox["v", "′",
MultilineFunction->None]\)[r]^2 + (-2 + r) r^2 
\!\(\*SuperscriptBox["v", "′′",
   r /. {v'[r_] :> Sum[-i c[i]/r^(i + 1), {i, 1, n}], 
    v''[r_] :> 
     Sum[i (i + 1) c[i]/r^(i + 2), {i, 1, n}]}, {r, ∞, 
   n - 1}], r^-1];
Do[c[i] = c[i] /. Simplify[Solve[veqexp[[i]] == 0, c[i]][[1]]]; 
Print, {i, 1, n}] ;
Table[c[i], {i, 1, n}]];

InfinityICs[ll_?IntegerQ, ωω_?NumericQ] := Module[{c2},
Do[c2[i] = 
Infinitycs[[i]] /. {l -> ll, ω -> ωω}, {i, 1,
vtrunc = Sum[c2[i]/r^i, {i, 1, ninfphase}];
init = 1/r Exp[I ωω rstar[r] + vtrunc] /. r -> rinf;
dinit = 
D[1/r Exp[I ωω rstar[r] + vtrunc], r] /. r -> rinf;
{init, dinit}]

Solve it

l_] := Φout[ωω, l] = Φ /. 
Block[{$MaxExtraPrecision = 100}, 
  NDSolve[{eq[ωω, l], Φ[rinf] == 
     N[InfinityICs[l, ωω][[1]], wp], Φ'[
      rinf] == 
     N[InfinityICs[l, ωω][[2]], 
      wp]}, Φ, {r, rinf, r0}, WorkingPrecision -> wp, 
   AccuracyGoal -> ac, MaxSteps -> ∞]][[1]];
  • 1
    $\begingroup$ You could partition the interval yourself ... $\endgroup$ Oct 12, 2012 at 20:39
  • $\begingroup$ Sounds a bit like the Burgers equation (although that is a PDE). What does your equation look like? What does your mathematica code look like? I wonder if it is because of your history length in case you are using the % for the previous step. (Try setting $HistoryLength=0) $\endgroup$
    – dearN
    Oct 12, 2012 at 21:26
  • $\begingroup$ I added the ODE, it's the Klein-Gordon equation on Schwarzschild. I don't use % at all in code. It's just that the solutions to these equations get mega oscillatory as omega gets bigger, and I integrate over such a large range of 10000 in r, that I think the interpolatingfunction being generated starts to hog lots and lots of memory, as mathematica samples tonnes of points perhaps. $\endgroup$
    – fpghost
    Oct 12, 2012 at 22:35
  • 1
    $\begingroup$ Do you mind including the boundary conditions perhaps. And the values of M and omega etc.. $\endgroup$
    – dearN
    Oct 12, 2012 at 22:39
  • 1
    $\begingroup$ To repeat my question earlier: why can't you do NDSolve[(* stuff *), {r, rinf, rinf}, (* more stuff *)]? $\endgroup$ Oct 12, 2012 at 22:59

1 Answer 1


I need to integrate over quite a large range (10,000) and that is part of the problem, but I actually only need the value of the InterpolatingFunction[] produced at the end point. Is there a way to tell Mathematica I just want this last point?

One way to go about it is to have the start and end of the integration interval be identical. Consider the following:

y[5] /. First @ NDSolve[{y'[x] == y[x] Cos[x + y[x]], y[0] == 1}, y, {x, 5, 5}]

To see that the approach saves space, here's a comparison:

yi = y /. First@NDSolve[{y'[x] == y[x] Cos[x + y[x]], y[0] == 1}, y, {x, 5, 5}];

yn = y /. First@NDSolve[{y'[x] == y[x] Cos[x + y[x]], y[0] == 1}, y, {x, 0, 5}];
  • $\begingroup$ thanks but that won't work for me. I should have mentioned that I need to set up some initial conditions at a certain radius (the dependent variable) and then integrate outwards to a large distance, and it's the point at large r where I need the solutions. But I don't need the rest of the 10,000, so Mathematica could dump it, as soon as it gets enough data to leap frog to the next point if you see what I mean. $\endgroup$
    – fpghost
    Oct 12, 2012 at 22:40
  • $\begingroup$ Right; that's why I suggested that you have the start and end points be the same. If you'll edit your question to show the NDSolve[] snippet that you say has taxed your memory, I can edit my answer to match. $\endgroup$ Oct 12, 2012 at 22:44
  • $\begingroup$ that's not possible, the IC for the solution I want have to start at a certain radius (the event horizon) and I integrate it out to the "infinity". The solution with ICs set at the infinity and integrated inwards are different, and there is no known way to set the ICs anywhere else..Sorry if that's not a clear explanation, but that much at least is a given in my problem: ICs must be at r=2, point I want must be at r=10000. Edited my post now that should be a working example. Try large omega values to observe the lag. $\endgroup$
    – fpghost
    Oct 12, 2012 at 22:55
  • $\begingroup$ NDSolve[] will still start at the initial condition you gave (in your new example, Φ[rH] == (* stuff *)); as mentioned in the docs, the integration interval supplied need not contain the value in the initial conditions, and NDSolve[] will go forward or backward as appropriate. That's why I don't understand why you say this won't work. Have you tried it? $\endgroup$ Oct 12, 2012 at 23:11
  • $\begingroup$ Oh, actually, sorry I misunderstood you, maybe that will work then if NDSolve behaves like that. I didn't realise it could skip forward/backward. Unfortunately I have to dash now, but I will try this ASAP. thanks. $\endgroup$
    – fpghost
    Oct 12, 2012 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.