The dsolve command with options numeric and method=dverk78 finds a numerical solution using a seventh-eighth order continuous Runge-Kutta method.

The following options are available for the basic use form of the dverk78 method:

'output'

Specifies the desired output from dsolve, and the known option specifies user-defined known functions. For more information, see dsolve[numeric].

'number', 'procedure', 'start', 'initial', and 'procvars'

These options are used to specify the IVP using procedures. For more information, see dsolve[numeric,IVP].

'startinit','implicit', and 'optimize'

These options control the method and behavior of the computation. For more information on the first two, see dsolve[numeric,IVP], for the last, see dsolve[numeric].

'abserr', 'relerr', 'initstep', 'minstep', and 'maxstep'

Specify the desired accuracy of the solution, and allow for more detailed control of the step size. For more information, see dsolve[Error_Control]. The default values are $\mathrm{abserr}=1.\times {10}^{\mathrm{-8}}$, $\mathrm{relerr}=1.\times {10}^{\mathrm{-8}}$, $\mathrm{minstep}=0$, $\mathrm{maxstep}=2.0$, and $\mathrm{initstep}=\mathrm{maxstep}{\mathrm{relerr}}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$6$}\right.}$.

The dverk78 method is capable of working in arbitrary precision based on the setting of Digits, and can be used to obtain high accuracy solutions for ODE systems. As a note of caution, however, it is often necessary to work with a greater number of Digits than would be expected for the requested error tolerance, to prevent round-off error from giving an inaccurate error estimate. When the error tolerance is too strict for the current setting of Digits, it is detected by the algorithm, and an error is issued. Asymptotically dverk78 requires Digits set so that $\mathrm{tol}=\mathrm{O}\left({\mathrm{Float}\left(1\,-\mathrm{Digits}\right)}^{\raisebox{1ex}{$9$}\!\left/ \!\raisebox{-1ex}{$8$}\right.}\right)$.

Note: As for all other numeric dsolve methods, any setting of Digits that is less than or equal to hardware precision works in hardware precision, so for $\mathrm{Digits}<\mathrm{evalhf}\left(\mathrm{Digits}\right)$, the algorithm works in $\mathrm{trunc}\left(\mathrm{evalhf}\left(\mathrm{Digits}\right)\right)$ Digits (see evalhf).

'errorest'

Allows control over the method used to obtain error estimates, and can have the values $\mathrm{interpolant}$ or pair. The $\mathrm{interpolant}$ option is the default, and tells dverk78 to estimate the error using the interpolant of the method. The pair option tells dverk to estimate the error using the difference between the computed solution and a lower order estimate of the solution (hence 'pair'). Typically, the $\mathrm{interpolant}$ option gives a better error estimate than the pair option, but requires more work to compute. This may also be a better option for expensive systems, as dverk78 with the $\mathrm{interpolant}$ option performs 20 function evaluations per step, while the pair option only requires 13.

'maxfun'

An integer that specifies the maximum on the number of evaluations of the ODE or ODE system in one call to the returned procedure. This limits the amount of work done on any individual call (See dsolve[maxfun]). This option is relevant for the procedure-style output only. By default, this value is set to zero (disabled).

Results can be plotted using the function odeplot in the plots package.

$\mathrm{dsol1}≔\mathrm{dsolve}\left(\left\{y\left(0\right)=1\&comma;\mathrm{D}\left(y\right)\left(x\right)=y\left(x\right)\right\}\&comma;\mathrm{numeric}\&comma;\mathrm{method}=\mathrm{dverk78}\&comma;\mathrm{abserr}=1.\times {10}^{\mathrm{-8}}\&comma;\mathrm{relerr}=1.\times {10}^{\mathrm{-8}}\right)$

$\mathrm{dsol1}≔\mathbf{proc}\left(\mathrm{x\_dverk78}\right)...\mathbf{end\; proc}$

$\mathrm{dsol1}\left(0\right)$

$\left[x=0.\&comma;y\left(x\right)=1.\right]$

$\mathrm{dsol1}\left(0.4\right)$

$\left[x=0.400000000000000\&comma;y\left(x\right)=1.49182469764117\right]$

$\mathrm{dsol1}\left(1.2\right)$

$\left[x=1.20000000000000\&comma;y\left(x\right)=3.32011692267802\right]$

$\mathrm{dsol1}\left(1.5\right)$

$\left[x=1.50000000000000\&comma;y\left(x\right)=4.48168906989179\right]$

$\mathrm{dsys2}≔\left\{\mathrm{diff}\left(x\left(t\right)\&comma;t\right)=-y\left(t\right)\&comma;\mathrm{diff}\left(y\left(t\right)\&comma;t\right)=x\left(t\right)\right\}$

$\mathrm{dsys2}≔\left\{\frac{\&DifferentialD;}{\&DifferentialD;t}x\left(t\right)=-y\left(t\right)\&comma;\frac{\&DifferentialD;}{\&DifferentialD;t}y\left(t\right)=x\left(t\right)\right\}$

$\mathrm{init2}≔\left\{x\left(0\right)=1\&comma;y\left(0\right)=0\right\}$

$\mathrm{init2}≔\left\{x\left(0\right)=1\&comma;y\left(0\right)=0\right\}$

$\mathrm{dsol2}≔\mathrm{dsolve}\left(\mathrm{dsys2}\mathbf{union}\mathrm{init2}\&comma;\mathrm{numeric}\&comma;\mathrm{method}=\mathrm{dverk78}\&comma;\mathrm{output}=\mathrm{Array}\left(\left[0\&comma;0.7\&comma;1.1\&comma;1.5\&comma;2.2\&comma;2.5\right]\right)\right)$

$\mathrm{dsol2}≔\left[\begin{array}{c}\left[\begin{array}{ccc}t& x\left(t\right)& y\left(t\right)\end{array}\right]\\ \left[\begin{array}{ccc}0.& 1.& 0.\\ 0.7& 0.764842187284791& 0.644217687237135\\ 1.1& 0.453596121533758& 0.891207359932299\\ 1.5& 0.0707372017908811& 0.997494986546686\\ 2.2& \mathrm{-0.588501116104254}& 0.808496403994437\\ 2.5& \mathrm{-0.801143614267719}& 0.598472144707831\end{array}\right]\end{array}\right]$

$\left[\mathrm{dsol2}\left[1,1\right]\left[1\right]=\mathrm{dsol2}\left[2,1\right]\left[6,1\right]\&comma;\mathrm{dsol2}\left[1,1\right]\left[2\right]=\mathrm{dsol2}\left[2,1\right]\left[6,2\right]\right]$

$\left[t=2.5\&comma;x\left(t\right)=\mathrm{-0.801143614267719}\right]$

$\mathrm{dsys3}≔\left\{\mathrm{diff}\left(x\left(t\right)\&comma;\mathrm{`\$`}\left(t\&comma;3\right)\right)=x\left(t\right)+y\left(t\right)\&comma;\mathrm{diff}\left(y\left(t\right)\&comma;t\right)=-x\left(t\right)+y\left(t\right)\right\}$

$\mathrm{dsys3}≔\left\{\frac{{\&DifferentialD;}^3}{\&DifferentialD;t^3}x\left(t\right)=x\left(t\right)+y\left(t\right)\&comma;\frac{\&DifferentialD;}{\&DifferentialD;t}y\left(t\right)=-x\left(t\right)+y\left(t\right)\right\}$

$\mathrm{init3}≔\left\{x\left(0\right)=1\&comma;y\left(0\right)=1.3\&comma;\mathrm{D}\left(x\right)\left(0\right)=-1\&comma;{\mathrm{D}}^{\left(2\right)}\left(x\right)\left(0\right)=0.5\right\}$

$\mathrm{init3}≔\left\{x\left(0\right)=1\&comma;y\left(0\right)=1.3\&comma;\mathrm{D}\left(x\right)\left(0\right)=\mathrm{-1}\&comma;{\mathrm{D}}^{\left(2\right)}\left(x\right)\left(0\right)=0.5\right\}$

$\mathrm{dsol3}≔\mathrm{dsolve}\left(\mathrm{dsys3}\mathbf{union}\mathrm{init3}\&comma;\mathrm{numeric}\&comma;\mathrm{method}=\mathrm{dverk78}\&comma;\mathrm{output}=\mathrm{listprocedure}\&comma;\mathrm{maxstep}=0.01\right)\&colon;$

$\mathrm{dsol3x}≔\mathrm{subs}\left(\mathrm{dsol3}\&comma;x\left(t\right)\right)\&colon;$$\mathrm{dsol3y}≔\mathrm{subs}\left(\mathrm{dsol3}\&comma;y\left(t\right)\right)\&colon;$

$\left[\mathrm{dsol3x}\left(0\right)\&comma;\mathrm{dsol3y}\left(0\right)\right]$

$\left[1.\&comma;1.3\right]$

$\left[\mathrm{dsol3x}\left(1\right)\&comma;\mathrm{dsol3y}\left(1\right)\right]$

$\left[0.623059643663728\&comma;2.30932410466714\right]$

$\left[\mathrm{dsol3x}\left(1.5\right)\&comma;\mathrm{dsol3y}\left(1.5\right)\right]$

$\left[1.36464098811480\&comma;3.23223976483577\right]$

$\mathrm{dsol4a}≔\mathrm{dsolve}\left(\left\{\mathrm{diff}\left(y\left(x\right)\&comma;x\&comma;x\right)=-y\left(x\right)\&comma;y\left(0\right)=0\&comma;\mathrm{D}\left(y\right)\left(0\right)=1\right\}\&comma;\mathrm{numeric}\&comma;\mathrm{method}=\mathrm{dverk78}\&comma;\mathrm{tolerance}=1.\times {10}^{\mathrm{-9}}\&comma;\mathrm{errorest}=\mathrm{interpolant}\right)$

$\mathrm{dsol4a}≔\mathbf{proc}\left(\mathrm{x\_dverk78}\right)...\mathbf{end\; proc}$

$\mathrm{tt}≔\mathrm{time}\left(\right)\&colon;$

$\mathrm{dsol4a}\left(1000\right)$

$\left[x=1000.\&comma;y\left(x\right)=0.826879438181235\&comma;\frac{\&DifferentialD;}{\&DifferentialD;x}y\left(x\right)=0.562379146833518\right]$

$\mathrm{time}\left(\right)-\mathrm{tt}$

$0.220$

$\mathrm{dsol4b}≔\mathrm{dsolve}\left(\left\{\mathrm{diff}\left(y\left(x\right)\&comma;x\&comma;x\right)=-y\left(x\right)\&comma;y\left(0\right)=0\&comma;\mathrm{D}\left(y\right)\left(0\right)=1\right\}\&comma;\mathrm{numeric}\&comma;\mathrm{method}=\mathrm{dverk78}\&comma;\mathrm{tolerance}=1.\times {10}^{\mathrm{-9}}\&comma;\mathrm{errorest}=\mathrm{pair}\right)$

$\mathrm{dsol4b}≔\mathbf{proc}\left(\mathrm{x\_dverk78}\right)...\mathbf{end\; proc}$

$\mathrm{tt}≔\mathrm{time}\left(\right)\&colon;$

$\mathrm{dsol4b}\left(1000\right)$

$\left[x=1000.\&comma;y\left(x\right)=0.826879446231662\&comma;\frac{\&DifferentialD;}{\&DifferentialD;x}y\left(x\right)=0.562379141799196\right]$

$\mathrm{time}\left(\right)-\mathrm{tt}$

$0.118$

Enright, W.H. "The Relative Efficiency of Alternative Defect Control Schemes for High Order Continuous Runge-Kutta Formulas." Technical Report 252/91, Dept. of Computer Science, University of Toronto, June, 1991.

Verner, J.H. "Explicit Runge-Kutta Methods with Estimates of the Local Truncation Error." SIAM Journal of Numerical Analysis, Aug. 1978.

Much appreciation and thanks go to Dr. Jim Verner who provided us with the coefficients for the dverk78 method, and the related interpolant. Without these coefficients the method would be restricted to hardware precision.