Re: NDSolve useless?
- To: mathgroup at smc.vnet.net
- Subject: [mg63917] Re: NDSolve useless?
- From: Alberto Verga <Alberto.Verga at laposte.net>
- Date: Sun, 22 Jan 2006 00:52:31 -0500 (EST)
- References: <dqd9oi$mc5$1@smc.vnet.net> <dqsphh$ch9$1@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
Hi,
Thank you very much for this precise and relevant answer. Trying different
solutions I used vector valued functions, but I did not use a full compiled
rhs to introduce into NDSolve. Your example is for me very instructive (I
have little knowledge about numerical computations using mathematica.)
(To complete the program, add:
ListPlot[Transpose[xy[tf] /. sol], PlotJoined -> True];
you'll get a nice spiral.)
Regards,
Alberto.
<rknapp at wolfram.com> a écrit dans le message de news:
dqsphh$ch9$1 at smc.vnet.net...
> The problem here is that you are comparing apples and oranges.
>
> ode45 requires that you set up a function that returns the right hand
> side of the ODEs in first order normal form, which I assume is what you
> used to
> make comparisons.
>
> You have the option to do this in Mathematica (though it is not required),
> and as I show below, there is very little speed difference between NDSolve
> and ode45 with this specification.
>
> The reason Mathematica is so much slower for the problem specification you
> have set up is that you have created a large symbolic expression that
> Mathematica is required to process and optimize. We have made
> improvements so that this processing will be faster in the next release
> of Mathematica, but even so, it is still much faster, particularly in
> a problem of this type, to just specify the right hand sides of the
> equations:
>
> Define:
>
> NV = 128;
> tf = 1.0;
> delta = 0.1;
> a = 0.01;
> x0=Table[n/NV+a Sin[2p n/NV],{n,1,NV}];
> y0=Table[-a Sin[2p n/NV],{n,1,NV}];
>
> (* This gives the right hand side of the equations *)
> cf = Compile[{{x, _Real,
> 1}, {y, _Real, 1}}, Module[{difx, dify, deno,NV = Length[x]},
> Transpose[Table[
> difx = 2.*Pi*(x[[n]] - x);
> dify = 2.*Pi*(y[[n]] -y);
> deno = 2.*NV*(Cosh[dify] - Cos[difx] +delta^2);
> {Total[-Sinh[dify]/deno], Total[Sin[difx]/deno]},
> {n, 1,NV}]]
> ]];
>
>
> (* Using this prevents evaluation except during the solution process *)
> g[{x_?VectorQ, y_?VectorQ}] := cf[x, y];
>
> then, I get
>
> In[10]:=
> Timing[sol=First[NDSolve[{xy'[t] == g[xy[t]],xy[0] == {x0,
> y0}},xy,{t,0,tf}]]]
>
> Out[10]=
> {2.984 Second,{xy->InterpolatingFunction[{{0.,1.}},<>]}}
>
> the reason that this is slower than your time for ode45 is that the
> default tolerances
> for Mathematica are set more stringently. Essentially the some
> tolerance settings
> as for ode45 can be achieved with
>
> In[11]:=
> Timing[sol=First[NDSolve[{xy'[t] == g[xy[t]],xy[0] == {x0, y0}},xy,
> {t,0,tf}, PrecisionGoal->3, AccuracyGoal->6]]]
>
> Out[11]=
> {1.39 Second,{xy->InterpolatingFunction[{{0.,1.}},<>]}}
>
> which is quite comparable in time.
>