MathGroup Archive 2012

[Date Index] [Thread Index] [Author Index]

Search the Archive

NDSolve very very slow

I'm using NDSolve to solve the ODE:


\[Lambda][l_] = l (l + 1);
rinf = 15000;
$MinPrecision = 40;
wp = $MinPrecision;
ac = $MinPrecision - 8;
pg = wp/2;

   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;

(*The solution:*)

\[CapitalPhi]out[\[Omega]_, l_] := \[CapitalPhi]out[\[Omega], l] = {\[CapitalPhi], \[CapitalPhi]'} /. 
   Block[{$MaxExtraPrecision = 100}, 
     NDSolve[{eq[\[Omega], l], \[CapitalPhi][rinf] == 
        init, \[CapitalPhi]'[rinf] == 
        dinit}, {\[CapitalPhi], \[CapitalPhi]'}, {r, 29,
        39}, WorkingPrecision -> wp, 
      AccuracyGoal -> ac, PrecisionGoal -> pg, 
      MaxSteps -> \[Infinity]]][[1]];

(*Run as:*)

\[CapitalPhi]out[2, 1]//AbsoluteTiming

for example, using the appropriate init, dinit from below.

The initial conditions, init and dinit, depend on the parameters \[Omega] and l,
and I compute them in a complicated way, but to save on posting all that code I will
give them here for just a sample of the troubling values where NDSolve is super slow (namely 
the large \[Omega] gets the slower it gets because of the solution being very oscillatory):


\[Omega]=2, l=1

init=-0.00003158476708483920624537825468234683185264 + 0.00005870985383714427760467962946161878449483 I; 
dinit=-0.0001174332599414352606295233010134701480989 - 0.00006318187181746517318131907026619247606411 I;

\[Omega]=15, l=5

init=-0.00006498016608463250053612003781430969895230 + 0.00001490041816106928145235191804618494074335 I;
dinit=-0.0002235317451490380993971487099054187817855 - 0.0009748334620044207867560932198118869943621 I;

\[Omega]=30, l=7

init=0.00006000392711252709381442292707073690333932 - 0.00002905121643121747797219318220828295775626 I;
dinit=0.0008716487129159803540979973990019187444996 +0.001800359797589481474882016847763171211252 I;

The time that NDSolve takes to solve this system goes from 5minutes in the \[Omega]=2 case, to
3 or 4hrs in the \[Omega]=30 case; I really need to speed this up drastically.

I can't help feel that C++ or another system may do it much quicker, but before I take that leap, I would
like to ensure Mathematica can't be speeded up for this calculation.

1) Could 'Compile' somehow be used to speed this up? I know supposedley Mathematica compiles 
automatically to byte code in NDSolve, and there is the Compile->True option (which I have tried
without any speed increase), but what about explicitely compiling the functions being fed into
NDSolve (I'm thinking of something like what was done here, could anyone show me how to adapt that to my code?)

2) Specifying somehow that my functions are only numeric and not symbolic with attributes?

3) Different 'Method' in NDSolve?

I should mention that for large \[Omega] I do need the result to be reliable to something like 14dp, and
that why I am working with high WorkingPrecision etc.

Thanks for any help and suggestions.

  • Prev by Date: Re: Eigenvalues works very slow
  • Next by Date: Re: Eigenvalues works very slow
  • Previous by thread: Re: CDF Security
  • Next by thread: Re: NDSolve very very slow