Efficiency in and with dsolve/numeric solution procedures

Description


•

There are many considerations with respect to efficiency when dealing with numerical ODE solutions in Maple. This page serves as a central repository for information on the efficient use of dsolve/numeric for computation of solution values, plots, and general consideration for all dsolve/numeric procedureform solutions.



Digits


•

The setting of the environment variable Digits is critical for the efficiency of dsolve/numeric. When Digits is set to the value of evalhf(Digits) (currently 15 on all supported platforms) or smaller, the computation of the solution proceeds using hardware floatingpoint operations, using 53 binary Digits of precision.


When Digits is 15 or smaller, we say that the computation is running in hardware mode, and when possible, evalhf, a hardwarefloat specific fast evaluator, is used. This can provide a significant efficiency boost, speeding up computation by a factor of as much as in some cases.


Furthermore, for the methods that support it, the compile option can be used when operating in hardware mode, which can speed up computation by a further factor of , on top of the .

Run with software precision (16 Digit base10)
>

dsn16 := dsolve({diff(y(t),t,t)+y(t),y(0)=0,D(y)(0)=1}, numeric, abserr=1e12, relerr=1e12, maxfun=0);

 (1) 
 (2) 
 (3) 
Run with hardware precision
>

dsn10 := dsolve({diff(y(t),t,t)+y(t),y(0)=0,D(y)(0)=1}, numeric, abserr=1e12, relerr=1e12, maxfun=0);

 (4) 
 (5) 
 (6) 
Compare:
 (7) 
Run with hardware precision and compile
>

dsn10 := dsolve({diff(y(t),t,t)+y(t),y(0)=0,D(y)(0)=1}, numeric, abserr=1e12, relerr=1e12, maxfun=0, compile=true);

 (8) 
 (9) 
 (10) 
Compare:
 (11) 


Plotting


•

Most procedures output by dsolve/numeric are computeasneeded, meaning that the system is integrated either from the initial point or from the last computed solution point as solution values are requested.


The word most is used here because certain problems or options force or imply storage of the solution. These include solutions for BVP problems, IVP solutions using the range option for precomputing the solution, or solutions using the piecewise output form.


One serious issue with the computeasneeded strategy is that it is not viable to reverse the integration direction, as the stability of the IVP problem is not known.


For this reason (among others), there is a special purpose IVP solution plotting routine, plots[odeplot], which produces efficient plots for ODE solutions. The use of the standard plot routines, like plot, can be inefficient as they have no way of knowing about the restriction to the order of the points being computed. In addition, a special purpose routine allows for specialized output modes for the ODE integrators, supporting a tighter (and even more efficient) bond between the solution process and the plotting.

>

dsn := dsolve({diff(y(t),t,t)+y(t),y(0)=0,D(y)(0)=1}, numeric, abserr=1e12, relerr=1e12, maxfun=0, output=listprocedure);

 (12) 
plot using specialized odeplot
>

plots[odeplot](dsn,[t,y(t)],0..50*Pi,numpoints=500);

 (13) 
plot using generic plot (adaptive)
>

plot(op([2,2],dsn),0..50*Pi);

 (14) 
Compare:
 (15) 


The implicit option for systems


•

For IVP problems the input system needs to be converted to solved first order form, as described in dsolve/numeric/IVP. For some systems this symbolic solution step can be prohibitively expensiveespecially for dense systems that are autogenerated, or for systems that arise from a homotopy problem.


Fortunately, for the cases where the system is linear in the derivatives (after first order conversion) the implicit option can be used. The implicit option delays the symbolic backsolve (performed once at the time when the solution procedure is being formulated) into many purely numerical backsolves (performed at each evaluation of the derivatives during the integration).


In addition to avoiding symbolic expression swell, and the associated roundoff error resulting from large expressions, this approach also has the advantage that it also chooses a more stable choice of solving variables in the elimination process, due to pivoting on the local numerical values of the system coefficients instead of on the unknown symbolic coefficients.

To illustrate the idea, consider the first order ODE system:
>

eq1 := cos(t)*diff(x(t),t)+(1/2sin(t))*diff(y(t),t)=1/2y(t);

 (16) 
>

eq2 := sin(t)*diff(x(t),t)+cos(t)*diff(y(t),t)=x(t);

 (17) 
Direct symbolic solution might use the first equation to solve for and the second for , but this is a problem when passes through zero. Note that dsolve/numeric actually performs a solve with backsolve, and some limited simplifications here though, so this system does not actually pose a problem.
The implicit option will, at each point where the values of and are required, evaluate the coefficients then solve, avoiding this type of issue.
>

dsn := dsolve({eq1,eq2,x(0)=0,y(0)=0}, numeric, implicit=true);

 (18) 
>

plots[odeplot](dsn,[x(t),y(t)],0..6*Pi);



Solutions depending on other solutions


•

Sometimes it is necessary to compute a dsolve/numeric solution of an ODE problem that depends upon another dsolve/numeric solution. For example, suppose you obtain a solution for from a dsolve/numeric problem that is then used as part of another problem as, say, a forcing function.


When the time scales are equivalent, a solution can be obtained much more efficiently by simply combining the two systems. This is because a system that depends on a dsolve/numeric solution procedure cannot be integrated under evalhf (see the first section), so it will integrate much more slowly.

For example, solving separately:
>

sys1 := {diff(x(t),t,t)+x(t)=0, x(0)=0, D(x)(0)=1};

 (19) 
>

sol1 := dsolve(sys1, numeric, output=listprocedure);

 (20) 
 (21) 
>

sys2 := {diff(y(t),t,t)=xf(t), y(0)=0, D(y)(0)=1};

 (22) 
>

sol2 := dsolve(sys2, numeric, known={xf});

 (23) 
 (24) 
>

t_separate := time()tt;

 (25) 
Now combining:
>

sysc := sys1 union eval(sys2,xf(t)=x(t));

 (26) 
>

solc := dsolve(sysc, numeric);

 (27) 
 (28) 
>

t_combined := time()tt;

 (29) 
Compare:
 (30) 

When the time scales are different (for example, the second system depends upon , where is the solution of another ODE problem), then for the default stiff and nonstiff methods, piecewise form output is available for a prespecified range. Note that it may be necessary to compute the piecewise solution well past the desired integration limit for the second problem, as the solver may step past the current region during the integration.

From the same example:
>

sol1p := dsolve(sys1, numeric, output=piecewise, range=0..110):

>

type(xp,specfunc(anything,piecewise)), nops(xp);

 (31) 
>

sys2p := {diff(y(t),t,t)=xf, y(0)=0, D(y)(0)=1}:

>

sol2p := dsolve(sys2p, numeric);

 (32) 
 (33) 
>

t_piecewise := time()tt;

 (34) 
Compare:
>

t_separate/t_piecewise;

 (35) 

Finally, a straightup range solution (no piecewise) can be used if desired. This will at least avoid recomputation of the solution, but will not run in evalhf mode, so it will be the least efficient of the alternatives.



See Also


Digits, dsolve/ck45, dsolve/classical, dsolve/dverk78, dsolve/gear, dsolve/lsode, dsolve/numeric, dsolve/numeric/IVP, dsolve/rkf45, dsolve/rosenbrock, evalhf, piecewise, plots[odeplot]

