Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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:

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

The ODE is:

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

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

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

Solve as

\[CapitalPhi]inExt[\[Omega]\[Omega]_ , 
 l_] := \[CapitalPhi]inExt[\[Omega]\[Omega], l] = \[CapitalPhi] /. 
 NDSolve[{eq[\[Omega]\[Omega], l], \[CapitalPhi][rH] == 
    N[HorizonICs[l, \[Omega]\[Omega]][[1]], wp], \[CapitalPhi]'[
     rH] == N[HorizonICs[l, \[Omega]\[Omega]][[2]], 
     wp]}, \[CapitalPhi], {r, rH, rinf}, WorkingPrecision -> wp, 
  AccuracyGoal -> ac, MaxSteps -> \[Infinity], 
  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;
\[Lambda][l_] = l (l + 1);
rinf = 10000;
ninfphase = 
50; 

Set init conditions:

Infinitycs = Module[{n = ninfphase, c}, 
Clear[c];
veqexp = 
CoefficientList[
 Series[(-2 - l r - l^2 r + 2 (r + I r^3 \[Omega]) 
\!\(\*SuperscriptBox["v", "\[Prime]",
MultilineFunction->None]\)[r] + (-2 + r) r^2 
\!\(\*SuperscriptBox["v", "\[Prime]",
MultilineFunction->None]\)[r]^2 + (-2 + r) r^2 
\!\(\*SuperscriptBox["v", "\[Prime]\[Prime]",
MultilineFunction->None]\)[r])/
   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, \[Infinity], 
   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, \[Omega]\[Omega]_?NumericQ] := Module[{c2},
Do[c2[i] = 
Infinitycs[[i]] /. {l -> ll, \[Omega] -> \[Omega]\[Omega]}, {i, 1,
 ninfphase}];
vtrunc = Sum[c2[i]/r^i, {i, 1, ninfphase}];
init = 1/r Exp[I \[Omega]\[Omega] rstar[r] + vtrunc] /. r -> rinf;
dinit = 
D[1/r Exp[I \[Omega]\[Omega] rstar[r] + vtrunc], r] /. r -> rinf;
Clear[c2];
{init, dinit}]

Solve it

\[CapitalPhi]out[\[Omega]\[Omega]_, 
l_] := \[CapitalPhi]out[\[Omega]\[Omega], l] = \[CapitalPhi] /. 
Block[{$MaxExtraPrecision = 100}, 
  NDSolve[{eq[\[Omega]\[Omega], l], \[CapitalPhi][rinf] == 
     N[InfinityICs[l, \[Omega]\[Omega]][[1]], wp], \[CapitalPhi]'[
      rinf] == 
     N[InfinityICs[l, \[Omega]\[Omega]][[2]], 
      wp]}, \[CapitalPhi], {r, rinf, r0}, WorkingPrecision -> wp, 
   AccuracyGoal -> ac, MaxSteps -> \[Infinity]]][[1]];
share|improve this question
1  
You could partition the interval yourself ... –  belisarius Oct 12 '12 at 20:39
    
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) –  drN Oct 12 '12 at 21:26
    
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. –  fpghost Oct 12 '12 at 22:35
1  
Do you mind including the boundary conditions perhaps. And the values of M and omega etc.. –  drN Oct 12 '12 at 22:39
    
To repeat my question earlier: why can't you do NDSolve[(* stuff *), {r, rinf, rinf}, (* more stuff *)]? –  J. M. Oct 12 '12 at 22:59
show 3 more comments

1 Answer

up vote 2 down vote accepted

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}]
   0.07731217497157500942

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}];
ByteCount[yi]
   1296

yn = y /. First@NDSolve[{y'[x] == y[x] Cos[x + y[x]], y[0] == 1}, y, {x, 0, 5}];
ByteCount[yn]
   3288
share|improve this answer
    
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. –  fpghost Oct 12 '12 at 22:40
    
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. –  J. M. Oct 12 '12 at 22:44
    
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. –  fpghost Oct 12 '12 at 22:55
    
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? –  J. M. Oct 12 '12 at 23:11
    
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. –  fpghost Oct 12 '12 at 23:14
show 10 more comments

Your Answer

 
discard

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.