NDSolve very very slow
- To: mathgroup at smc.vnet.net
- Subject: [mg128496] NDSolve very very slow
- From: popov.ghost at gmail.com
- Date: Thu, 25 Oct 2012 01:42:14 -0400 (EDT)
- Delivered-to: l-mathgroup@mail-archive0.wolfram.com
- Delivered-to: l-mathgroup@wolfram.com
- Delivered-to: mathgroup-newout@smc.vnet.net
- Delivered-to: mathgroup-newsend@smc.vnet.net
I'm using NDSolve to solve the ODE:
(*Definitions:*)
M=1;
\[Lambda][l_] = l (l + 1);
rinf = 15000;
$MinPrecision = 40;
wp = $MinPrecision;
ac = $MinPrecision - 8;
pg = wp/2;
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;
(*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):
e.g.
\[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
http://forums.wolfram.com/mathgroup/archive/2006/Jan/msg00373.html, 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.