Efficiency in and with dsolve/numeric solution procedures - Maple Programming Help

Home : Support : Online Help : Mathematics : Differential Equations : dsolve : numeric : dsolve/numeric/efficiency

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

Run with software precision (16 Digit base-10)

 > 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)\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{...}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathbf{end proc}}$ (1)
 > tt := time():
 > dsn16(50*Pi);
 $\left[{t}{=}{157.0796326794896}{,}{y}{}\left({t}\right){=}{-5.0356505}{}{{10}}^{{-13}}{,}\frac{{ⅆ}}{{ⅆ}{t}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{y}{}\left({t}\right){=}{1.000000000046060}\right]$ (2)
 > t_sw := time()-tt;
 ${\mathrm{t_sw}}{≔}{8.393}$ (3)

Run with hardware precision

 > 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)\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{...}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathbf{end proc}}$ (4)
 > tt := time():
 > dsn10(50*Pi);
 $\left[{t}{=}{157.079632679490}{,}{y}{}\left({t}\right){=}{-4.83536686168504}{}{{10}}^{{-14}}{,}\frac{{ⅆ}}{{ⅆ}{t}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{y}{}\left({t}\right){=}{1.00000000004606}\right]$ (5)
 > t_hw := time()-tt;
 ${\mathrm{t_hw}}{≔}{0.045}$ (6)

Compare:

 > t_sw/t_hw;
 ${186.5111111}$ (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=1e-12, relerr=1e-12, maxfun=0, compile=true);
 ${\mathrm{dsn10}}{:=}{\mathbf{proc}}\left({\mathrm{x_rkf45}}\right)\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{...}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathbf{end proc}}$ (8)
 > tt := time():
 > dsn10(50*Pi);
 $\left[{t}{=}{157.079632679490}{,}{y}{}\left({t}\right){=}{-4.83536686168504}{}{{10}}^{{-14}}{,}\frac{{ⅆ}}{{ⅆ}{t}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{y}{}\left({t}\right){=}{1.00000000004606}\right]$ (9)
 > t_hc := time()-tt;
 ${\mathrm{t_hc}}{≔}{0.015}$ (10)

Compare:

 > t_sw/t_hc, t_hw/t_hc;
 ${559.5333333}{,}{3.000000000}$ (11)

Plotting

 • 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)\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{...}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathbf{end proc}}{,}{y}{}\left({t}\right){=}{\mathbf{proc}}\left({t}\right)\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{...}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathbf{end proc}}{,}\frac{{ⅆ}}{{ⅆ}{t}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{y}{}\left({t}\right){=}{\mathbf{proc}}\left({t}\right)\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{...}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathbf{end proc}}\right]$ (12)

plot using specialized odeplot

 > tt := time():
 > plots[odeplot](dsn,[t,y(t)],0..50*Pi,numpoints=500); > t_spec := time()-tt;
 ${\mathrm{t_spec}}{≔}{0.052}$ (13)

plot using generic plot (adaptive)

 > tt := time():
 > plot(op([2,2],dsn),0..50*Pi); > t_gen := time()-tt;
 ${\mathrm{t_gen}}{≔}{3.438}$ (14)

Compare:

 > t_gen/t_spec;
 ${66.11538462}$ (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 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.

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

 > eq1 := cos(t)*diff(x(t),t)+(1/2-sin(t))*diff(y(t),t)=1/2-y(t);
 ${\mathrm{eq1}}{≔}{\mathrm{cos}}{}\left({t}\right){}\left(\frac{{ⅆ}}{{ⅆ}{t}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{x}{}\left({t}\right)\right){+}\left(\frac{{1}}{{2}}{-}{\mathrm{sin}}{}\left({t}\right)\right){}\left(\frac{{ⅆ}}{{ⅆ}{t}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{y}{}\left({t}\right)\right){=}\frac{{1}}{{2}}{-}{y}{}\left({t}\right)$ (16)
 > eq2 := sin(t)*diff(x(t),t)+cos(t)*diff(y(t),t)=x(t);
 ${\mathrm{eq2}}{≔}{\mathrm{sin}}{}\left({t}\right){}\left(\frac{{ⅆ}}{{ⅆ}{t}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{x}{}\left({t}\right)\right){+}{\mathrm{cos}}{}\left({t}\right){}\left(\frac{{ⅆ}}{{ⅆ}{t}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{y}{}\left({t}\right)\right){=}{x}{}\left({t}\right)$ (17)

Direct symbolic solution might use the first equation to solve for $\frac{ⅆ}{ⅆt}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}x\left(t\right)$ and the second for $\frac{ⅆ}{ⅆt}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}y\left(t\right)$, but this is a problem when $\mathrm{cos}\left(t\right)$ 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 $\frac{ⅆ}{ⅆt}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}x\left(t\right)$ and $\frac{ⅆ}{ⅆt}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}y\left(t\right)$ 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);
 ${\mathrm{dsn}}{:=}{\mathbf{proc}}\left({\mathrm{x_rkf45}}\right)\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{...}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathbf{end proc}}$ (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 $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.

For example, solving separately:

 > sys1 := {diff(x(t),t,t)+x(t)=0, x(0)=0, D(x)(0)=1};
 ${\mathrm{sys1}}{≔}\left\{\frac{{{ⅆ}}^{{2}}}{{ⅆ}{{t}}^{{2}}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{x}{}\left({t}\right){+}{x}{}\left({t}\right){=}{0}{,}{x}{}\left({0}\right){=}{0}{,}{\mathrm{D}}{}\left({x}\right){}\left({0}\right){=}{1}\right\}$ (19)
 > sol1 := dsolve(sys1, numeric, output=listprocedure);
 ${\mathrm{sol1}}{≔}\left[{t}{=}{\mathbf{proc}}\left({t}\right)\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{...}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathbf{end proc}}{,}{x}{}\left({t}\right){=}{\mathbf{proc}}\left({t}\right)\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{...}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathbf{end proc}}{,}\frac{{ⅆ}}{{ⅆ}{t}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{x}{}\left({t}\right){=}{\mathbf{proc}}\left({t}\right)\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{...}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathbf{end proc}}\right]$ (20)
 > xf := op([2,2],sol1);
 ${\mathrm{xf}}{:=}{\mathbf{proc}}\left({t}\right)\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{...}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathbf{end proc}}$ (21)
 > sys2 := {diff(y(t),t,t)=xf(t), y(0)=0, D(y)(0)=1};
 ${\mathrm{sys2}}{≔}\left\{\frac{{{ⅆ}}^{{2}}}{{ⅆ}{{t}}^{{2}}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{y}{}\left({t}\right){=}{\mathrm{xf}}{}\left({t}\right){,}{y}{}\left({0}\right){=}{0}{,}{\mathrm{D}}{}\left({y}\right){}\left({0}\right){=}{1}\right\}$ (22)
 > sol2 := dsolve(sys2, numeric, known={xf});
 ${\mathrm{sol2}}{:=}{\mathbf{proc}}\left({\mathrm{x_rkf45}}\right)\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{...}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathbf{end proc}}$ (23)
 > tt := time():
 > sol2(100);
 $\left[{t}{=}{100.}{,}{y}{}\left({t}\right){=}{200.506608804987}{,}\frac{{ⅆ}}{{ⅆ}{t}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{y}{}\left({t}\right){=}{1.13767792205232}\right]$ (24)
 > t_separate := time()-tt;
 ${\mathrm{t_separate}}{≔}{0.532}$ (25)

Now combining:

 > sysc := sys1 union eval(sys2,xf(t)=x(t));
 ${\mathrm{sysc}}{≔}\left\{\frac{{{ⅆ}}^{{2}}}{{ⅆ}{{t}}^{{2}}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{x}{}\left({t}\right){+}{x}{}\left({t}\right){=}{0}{,}\frac{{{ⅆ}}^{{2}}}{{ⅆ}{{t}}^{{2}}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{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\}$ (26)
 > solc := dsolve(sysc, numeric);
 ${\mathrm{solc}}{:=}{\mathbf{proc}}\left({\mathrm{x_rkf45}}\right)\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{...}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathbf{end proc}}$ (27)
 > tt := time():
 > solc(100);
 $\left[{t}{=}{100.}{,}{x}{}\left({t}\right){=}{-0.506372106885150}{,}\frac{{ⅆ}}{{ⅆ}{t}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{x}{}\left({t}\right){=}{0.862327180629681}{,}{y}{}\left({t}\right){=}{200.506372106886}{,}\frac{{ⅆ}}{{ⅆ}{t}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{y}{}\left({t}\right){=}{1.13767281937032}\right]$ (28)
 > t_combined := time()-tt;
 ${\mathrm{t_combined}}{≔}{0.008}$ (29)

Compare:

 > t_separate/t_combined;
 ${66.50000000}$ (30)
 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.

From the same example:

 > sol1p := dsolve(sys1, numeric, output=piecewise, range=0..110):
 > xp := op([2,2],sol1p):
 > type(xp,specfunc(anything,piecewise)), nops(xp);
 ${\mathrm{true}}{,}{1407}$ (31)
 > 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)\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{...}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathbf{end proc}}$ (32)
 > tt := time():
 > sol2p(100);
 $\left[{t}{=}{100.}{,}{y}{}\left({t}\right){=}{100.}{,}\frac{{ⅆ}}{{ⅆ}{t}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{y}{}\left({t}\right){=}{1.}\right]$ (33)
 > t_piecewise := time()-tt;
 ${\mathrm{t_piecewise}}{≔}{0.008}$ (34)

Compare:

 > t_separate/t_piecewise;
 ${66.50000000}$ (35)
 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.