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)))];