Re: NDSolve profiling
- To: mathgroup at smc.vnet.net
- Subject: [mg127017] Re: NDSolve profiling
- From: Ben <benjamin.r.lewis at gmail.com>
- Date: Sun, 24 Jun 2012 04:25:27 -0400 (EDT)
- Delivered-to: l-mathgroup@mail-archive0.wolfram.com
- References: <jrrvgu$4sd$1@smc.vnet.net> <jrup92$fnq$1@smc.vnet.net>
Narasimham, let me give you one specific example, hopefully someone
can help me to pinpoint my problem.
Here is a two-dimensional ODE defined by some complicated functions f1
and f2 (reproduced below). The outcome of NDSolve can be utterly
different depending on whether this problem is expressed in component
form or vector form:
laws={u''[t]==f1[u[t],v[t],u'[t],v'[t]],
v''[t]==f2[u[t],v[t],u'[t],v'[t]]};
Timing[NDSolve[Join[{u[0]==v[0]==0},
{u'[0]==10000,v'[0]==10000.0001},laws],{u,v},{t,-10,10},StoppingTest-
>v[t]^2<0.9` +u[t]^2]]
>> {8.003,{{u->InterpolatingFunction[{{-0.804964,0.804964}...
f3[x__List]:=({f1@@#1,f2@@#1}&)[Join[x]]; st[{u_,v_}]:=v^2<0.9` +u^2;
Timing[NDSolve[{y[0]=={0,0},y'[0]=={10000,
10000.0001},y''[t]==f3[y[t],y'[t]]},{y},{t,-10,10},StoppingTest-
>st[y[t]]]]
>> NDSolve::mxst: Maximum number of 10000 steps reached at the point t == 0.328...
>> {291.753,{{y->InterpolatingFunction[{{-0.827141,0.328212}...
Notice that the same functions are used in both forms. One form fails
early (at +0.3 instead of +0.8), and obviously uses far too many steps
even to get that far. It doesn't always have to be this way, for
example a different choice of initial conditions (like, y'=(0,1))
appears to produce the expected symmetry of behaviour between both
approaches. Any idea what's going on??
By the way, another difference is that in the component form, I
suspect f1 & f2 are only explicitly evaluated in the equation
processing stage, unlike the vector form. Since I wanted to exercise
finer control over how the guts of f1&f2 are evaluated during the
iterate stage, it is troublesome for me that NDSolve becomes
unreliable in that particular form.
Note in f1 & f2 below, that ProductLog appears repeatedly with the
same argument. In my real equations this occurs more than forty times,
so I want to make sure that it isn't being recalculated each time
(since NDSolve of this kind of equation is the main bottleneck in what
I'm doing).
-Ben
f2 = Function[{u, v, pu,
pv}, -pu ((
pv (-E^(1 + ProductLog[(u^2 - v^2)/E]) - u^2 + v^2) (2 u + (
2 u)/(1 + ProductLog[(u^2 - v^2)/E])))/(
2 (E^(1 + ProductLog[(u^2 - v^2)/E]) + u^2 - v^2)^2) + (
pu (-E^(1 + ProductLog[(u^2 - v^2)/E]) - u^2 + v^2) (-2 v - (
2 v)/(1 + ProductLog[(u^2 - v^2)/E])))/(
2 (E^(1 + ProductLog[(u^2 - v^2)/E]) + u^2 - v^2)^2)) -
pv ((pu (-E^(1 + ProductLog[(u^2 - v^2)/E]) - u^2 + v^2) (2 u + (
2 u)/(1 + ProductLog[(u^2 - v^2)/E])))/(
2 (E^(1 + ProductLog[(u^2 - v^2)/E]) + u^2 - v^2)^2) + (
pv (-E^(1 + ProductLog[(u^2 - v^2)/E]) - u^2 + v^2) (-2 v - (
2 v)/(1 + ProductLog[(u^2 - v^2)/E])))/(
2 (E^(1 + ProductLog[(u^2 - v^2)/E]) + u^2 - v^2)^2))]; f1 =
Function[{u, v, pu,
pv}, -pv (-((pv (2 u + (2 u)/(1 + ProductLog[(u^2 - v^2)/E])))/(
2 (E^(1 + ProductLog[(u^2 - v^2)/E]) + u^2 - v^2))) - (
pu (-2 v - (2 v)/(1 + ProductLog[(u^2 - v^2)/E])))/(
2 (E^(1 + ProductLog[(u^2 - v^2)/E]) + u^2 - v^2))) -
pu (-((pu (2 u + (2 u)/(1 + ProductLog[(u^2 - v^2)/E])))/(
2 (E^(1 + ProductLog[(u^2 - v^2)/E]) + u^2 - v^2))) - (
pv (-2 v - (2 v)/(1 + ProductLog[(u^2 - v^2)/E])))/(
2 (E^(1 + ProductLog[(u^2 - v^2)/E]) + u^2 - v^2)))];