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 procedure-form solutions.

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 floating-point 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 hardware-float specific fast evaluator, is used. This can provide a significant efficiency boost, speeding up computation by a factor of as much as $100$ 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 $10$, on top of the $100$.

Note that when operating in hardware mode, the setting of Digits, as long as it is 15 or smaller, is irrelevant, so a computation run with $\mathrm{Digits}=10$ will produce the same result as one run with $\mathrm{Digits}=5$ or $\mathrm{Digits}=14$ for example.

Digits := 16:

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

$\mathrm{dsn16}:=\mathbf{proc}\left(\mathrm{x\_rkf45}\right)...\mathbf{end\; proc}$

tt := time():

dsn16(50*Pi);

$\left[t=157.0796326794896\,y\left(t\right)=\mathrm{-5.0356505}\times {10}^{\mathrm{-13}}\,\frac{\ⅆ}{\ⅆt}y\left(t\right)=1.000000000046060\right]$

t_sw := time()-tt;

$\mathrm{t\_sw}≔6.026$

Digits := 10:

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

$\mathrm{dsn10}:=\mathbf{proc}\left(\mathrm{x\_rkf45}\right)...\mathbf{end\; proc}$

tt := time():

dsn10(50*Pi);

$\left[t=157.079632679490\,y\left(t\right)=\mathrm{-4.83536686168504}\times {10}^{\mathrm{-14}}\,\frac{\ⅆ}{\ⅆt}y\left(t\right)=1.00000000004606\right]$

t_hw := time()-tt;

$\mathrm{t\_hw}≔0.035$

t_sw/t_hw;

$172.1714286$

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

$\mathrm{dsn10}:=\mathbf{proc}\left(\mathrm{x\_rkf45}\right)...\mathbf{end\; proc}$

tt := time():

dsn10(50*Pi);

$\left[t=157.079632679490\,y\left(t\right)=\mathrm{-4.83536686168504}\times {10}^{\mathrm{-14}}\,\frac{\ⅆ}{\ⅆt}y\left(t\right)=1.00000000004606\right]$

t_hc := time()-tt;

$\mathrm{t\_hc}≔0.016$

t_sw/t_hc, t_hw/t_hc;

$376.6250000,2.187500000$

Most procedures output by dsolve/numeric are compute-as-needed, 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 compute-as-needed strategy is that it is not viable to reverse the integration direction, as the stability of the IVP problem is not known.

For example, if the initial point for an IVP problem is $t=0$, a value is requested at $t=1000$, and then a subsequent request for a value at $t=999$ is made, the IVP system needs to be integrated all the way from $t=0$ for both requests. Conversely if the request for the value at $t=999$ is made first, and the value at $t=1000$ second, then for the second request, the integration only needs to be done from $t=999..1000$, as no reversal of the integration direction needs to be made for the second request.

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=1e-12, relerr=1e-12, maxfun=0, output=listprocedure);

$\mathrm{dsn}≔\left[t=\mathbf{proc}\left(t\right)...\mathbf{end\; proc}\,y\left(t\right)=\mathbf{proc}\left(t\right)...\mathbf{end\; proc}\,\frac{\ⅆ}{\ⅆt}y\left(t\right)=\mathbf{proc}\left(t\right)...\mathbf{end\; proc}\right]$

tt := time():

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

t_spec := time()-tt;

$\mathrm{t\_spec}≔0.057$

tt := time():

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

t_gen := time()-tt;

$\mathrm{t\_gen}≔5.234$

t_gen/t_spec;

$91.82456140$

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 expensive--especially for dense systems that are auto-generated, 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 round-off 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.

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

$\mathrm{eq1}≔\cos \left(t\right)\left(\frac{\ⅆ}{\ⅆt}x\left(t\right)\right)+\left(\frac{1}{2}-\sin \left(t\right)\right)\left(\frac{\ⅆ}{\ⅆt}y\left(t\right)\right)=\frac{1}{2}-y\left(t\right)$

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

$\mathrm{eq2}≔\sin \left(t\right)\left(\frac{\ⅆ}{\ⅆt}x\left(t\right)\right)+\cos \left(t\right)\left(\frac{\ⅆ}{\ⅆt}y\left(t\right)\right)=x\left(t\right)$

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

$\mathrm{dsn}:=\mathbf{proc}\left(\mathrm{x\_rkf45}\right)...\mathbf{end\; proc}$

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

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 $y\left(t\right)$ 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.

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

$\mathrm{sys1}≔\left\{\frac{{\ⅆ}^2}{\ⅆt^2}x\left(t\right)+x\left(t\right)=0\,x\left(0\right)=0\,\mathrm{D}\left(x\right)\left(0\right)=1\right\}$

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

$\mathrm{sol1}≔\left[t=\mathbf{proc}\left(t\right)...\mathbf{end\; proc}\,x\left(t\right)=\mathbf{proc}\left(t\right)...\mathbf{end\; proc}\,\frac{\ⅆ}{\ⅆt}x\left(t\right)=\mathbf{proc}\left(t\right)...\mathbf{end\; proc}\right]$

xf := op([2,2],sol1);

$\mathrm{xf}:=\mathbf{proc}\left(t\right)...\mathbf{end\; proc}$

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

$\mathrm{sys2}≔\left\{\frac{{\ⅆ}^2}{\ⅆt^2}y\left(t\right)=\mathrm{xf}\left(t\right)\,y\left(0\right)=0\,\mathrm{D}\left(y\right)\left(0\right)=1\right\}$

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

$\mathrm{sol2}:=\mathbf{proc}\left(\mathrm{x\_rkf45}\right)...\mathbf{end\; proc}$

tt := time():

sol2(100);

$\left[t=100.\,y\left(t\right)=200.506608804987\,\frac{\ⅆ}{\ⅆt}y\left(t\right)=1.13767792205232\right]$

t_separate := time()-tt;

$\mathrm{t\_separate}≔0.883$

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

$\mathrm{sysc}≔\left\{\frac{{\ⅆ}^2}{\ⅆt^2}x\left(t\right)+x\left(t\right)=0\,\frac{{\ⅆ}^2}{\ⅆt^2}y\left(t\right)=x\left(t\right)\,x\left(0\right)=0\,y\left(0\right)=0\,\mathrm{D}\left(x\right)\left(0\right)=1\,\mathrm{D}\left(y\right)\left(0\right)=1\right\}$

solc := dsolve(sysc, numeric);

$\mathrm{solc}:=\mathbf{proc}\left(\mathrm{x\_rkf45}\right)...\mathbf{end\; proc}$

tt := time():

solc(100);

$\left[t=100.\,x\left(t\right)=\mathrm{-0.506372106885150}\,\frac{\ⅆ}{\ⅆt}x\left(t\right)=0.862327180629681\,y\left(t\right)=200.506372106886\,\frac{\ⅆ}{\ⅆt}y\left(t\right)=1.13767281937032\right]$

t_combined := time()-tt;

$\mathrm{t\_combined}≔0.015$

t_separate/t_combined;

$58.86666667$

When the time scales are different (for example, the second system depends upon $x\left(1-t\right)$, where $x\left(t\right)$ is the solution of another ODE problem), then for the default stiff and non-stiff methods, piecewise form output is available for a pre-specified 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.

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

xp := op([2,2],sol1p):

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

$\mathrm{true},1407$

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

sol2p := dsolve(sys2p, numeric);

$\mathrm{sol2p}:=\mathbf{proc}\left(\mathrm{x\_rkf45}\right)...\mathbf{end\; proc}$

tt := time():

sol2p(100);

$\left[t=100.\,y\left(t\right)=100.\,\frac{\ⅆ}{\ⅆt}y\left(t\right)=1.\right]$

t_piecewise := time()-tt;

$\mathrm{t\_piecewise}≔0.007$

t_separate/t_piecewise;

$126.1428571$

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