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

Online Help

All Products    Maple    MapleSim


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

Efficiency in and with dsolve/numeric solution procedures

 

Description

Digits

Plotting

The implicit option for systems

Solutions depending on other solutions

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 Digits=10 will produce the same result as one run with Digits=5 or 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);

dsn16:=procx_rkf45...end proc

(1)

tt := time():

dsn16(50*Pi);

t=157.0796326794896,yt=5.035650510-13,ⅆⅆtyt=1.000000000046060

(2)

t_sw := time()-tt;

t_sw:=4.459

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

dsn10:=procx_rkf45...end proc

(4)

tt := time():

dsn10(50*Pi);

t=157.079632679490,yt=4.8353668616850410-14,ⅆⅆtyt=1.00000000004606

(5)

t_hw := time()-tt;

t_hw:=0.023

(6)

Compare:

t_sw/t_hw;

193.8695652

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

dsn10:=procx_rkf45...end proc

(8)

tt := time():

dsn10(50*Pi);

t=157.079632679490,yt=4.8353668616850410-14,ⅆⅆtyt=1.00000000004606

(9)

t_hc := time()-tt;

t_hc:=0.008

(10)

Compare:

t_sw/t_hc, t_hw/t_hc;

557.3750000,2.875000000

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

dsn:=t=proct...end proc,yt=proct...end proc,ⅆⅆtyt=proct...end proc

(12)

plot using specialized odeplot

tt := time():

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

t_spec := time()-tt;

t_spec:=0.062

(13)

plot using generic plot (adaptive)

tt := time():

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

t_gen := time()-tt;

t_gen:=3.182

(14)

Compare:

t_gen/t_spec;

51.32258065

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

eq1:=costⅆⅆtxt+12sintⅆⅆtyt=12yt

(16)

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

eq2:=sintⅆⅆtxt+costⅆⅆtyt=xt

(17)

Direct symbolic solution might use the first equation to solve for ⅆⅆtxt and the second for ⅆⅆtyt, but this is a problem when cost 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 ⅆⅆtxt and ⅆⅆtyt 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);

dsn:=procx_rkf45...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 yt 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};

sys1:=ⅆ2ⅆt2xt+xt=0,x0=0,Dx0=1

(19)

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

sol1:=t=proct...end proc,xt=proct...end proc,ⅆⅆtxt=proct...end proc

(20)

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

xf:=proct...end proc

(21)

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

sys2:=ⅆ2ⅆt2yt=xft,y0=0,Dy0=1

(22)

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

sol2:=procx_rkf45...end proc

(23)

tt := time():

sol2(100);

t=100.,yt=200.506608804987,ⅆⅆtyt=1.13767792205232

(24)

t_separate := time()-tt;

t_separate:=0.758

(25)

Now combining:

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

sysc:=ⅆ2ⅆt2xt+xt=0,ⅆ2ⅆt2yt=xt,x0=0,y0=0,Dx0=1,Dy0=1

(26)

solc := dsolve(sysc, numeric);

solc:=procx_rkf45...end proc

(27)

tt := time():

solc(100);

t=100.,xt=0.506372106885150,ⅆⅆtxt=0.862327180629681,yt=200.506372106886,ⅆⅆtyt=1.13767281937032

(28)

t_combined := time()-tt;

t_combined:=0.004

(29)

Compare:

t_separate/t_combined;

189.5000000

(30)
  

When the time scales are different (for example, the second system depends upon x1t, where xt 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);

true,1407

(31)

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

sol2p := dsolve(sys2p, numeric);

sol2p:=procx_rkf45...end proc

(32)

tt := time():

sol2p(100);

t=100.,yt=100.,ⅆⅆtyt=1.

(33)

t_piecewise := time()-tt;

t_piecewise:=0.005

(34)

Compare:

t_separate/t_piecewise;

151.6000000

(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.

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]

 


Download Help Document

Was this information helpful?



Please add your Comment (Optional)
E-mail Address (Optional)
What is ? This question helps us to combat spam