dsolve/numeric
find numerical solution of ordinary differential equations
Calling Sequence
Parameters
Description
Examples
References
dsolve(odesys, numeric, vars, options)
dsolve(numeric, procopts, options)
odesys
-
set or list; ordinary differential or differential-algebraic equation(s) and initial/boundary conditions
numeric
name; instruct dsolve to find a numerical solution
vars
(optional) any indeterminate function of one variable, or a set or list of them, representing the unknowns of the ODE problem
procopts
(required if odesys is not present) options that specify a procedure defined system (procedure, initial, start, number, and procvars). For more information, see dsolve[numeric,IVP].
options
(optional) equations of the form keyword = value
The dsolve command with the numeric or type=numeric option finds a numerical solution for the ODE or ODE system. If the optional equation method=numericmethod is provided (where numericmethod is one of rkf45, ck45, rosenbrock, bvp, rkf45_dae, ck45_dae, rosenbrock_dae, dverk78, lsode, gear, taylorseries, mebdfi, or classical), dsolve uses that method to obtain the numerical solution.
Both initial and boundary value problems can be numerically solved, as well as initial differential algebraic problems. For a system input, the type of problem is automatically detected, but for some DAE problems it may be necessary to specify the method.
By default, a procedure is returned that can be used to obtain solution values if given the value of the independent variable.
The default for initial value problems (IVP) is a Runge-Kutta Fehlberg method that produces a fifth order accurate solution (For more information, see rkf45 or numeric,IVP). For boundary value problems (BVP), a finite difference technique with Richardson extrapolation is used (For more information, see numeric,BVP). For differential-algebraic IVP problems (DAE), a modification of the rkf45 method is used (For more information, see DAE_extension).
All IVP methods can be used for complex-valued IVPs with a real-valued independent variable (though for the default stiff and nonstiff IVP methods, it may be necessary to specify that the problem is complex via the complex option). None of the BVP or DAE methods can currently be used for complex-valued problems, requiring the system be converted to a real system before calling dsolve.
The IVP and DAE methods have additional capabilities for the returned solution procedure. Specifically, it is possible to query the last computed solution value (useful for problems with singularities), query the initial data, and change the initial data for some solution procedures (DAE methods support a subset of this). The IVP and DAE IVP integrators can also handle parametrized problems. For more information, see numeric,interactive as well as numeric,IVP and numeric,DAE.
When the set of dependent variables vars is specified as a list, then this specifies the order of the dependent variables as they appear in the output. The order is otherwise alphabetical in the dependent variable name.
Options
The high-level options, most common to IVPs, BVPs, and DAEs are as follows.
'output' =
keyword or array
'stiff' =
boolean
'events' =
list
'event_pre' =
keyword
'event_maxiter' =
integer
'range' =
numeric..numeric
'abserr' =
numeric or list
'relerr' =
'known' =
name or list of names
'optimize' =
'output'= keyword or array
Keyword that can take the values procedurelist, listprocedure, operator, or piecewise, or an array or Array that gives the values of the independent variable at which solution values are desired.
The default keyword is procedurelist, which gives the output from dsolve as a procedure. This procedure accepts the value of the independent variable as an argument, and it returns a list of the solution values of the form variable=value, where the left-hand sides are the names of the independent variable, the dependent variable(s) and their derivatives (for higher order equations), and the right-hand sides are the corresponding computed solution values.
The listprocedure keyword gives the output as a list of equations of the form variable=procedure, where the left-hand sides are the names of the independent variable, the dependent variable(s) and derivatives, and the right-hand sides are procedures that can be used to compute the corresponding solution components. This output form is most useful when the returned procedure is to be used later with, for example, fsolve.
The operator keyword gives the output as a list of equations of the form operator=procedure, where the left-hand sides are operators that can be applied to an independent variable value to give the function or derivative evaluated at a point, and the right-hand sides are procedures that can be used to compute the corresponding solution components. This output form is most useful for shortcut evaluation forms where the returned list is to be evaluated at a point (for example, result⁡1).
The piecewise keyword is available only for non-stiff and stiff default IVP and DAE methods (rkf45, ck45, rosenbrock, rkf45_dae, ck45_dae, and rosenbrock_dae) and the taylorseries method. It provides output as a list of equations of the form variable=pwfunc, where the left-hand sides are the names of the independent variable, the dependent variable(s) and derivatives, and the right-hand sides are piecewise polynomial functions describing the corresponding solution components. These piecewise functions are obtained from the method interpolants for each step of the computation. In addition, the form of the piecewise polynomials can be specified using an index on the piecewise output request. Specification of output=piecewise[horner] provides output in terms of horner-form polynomials (the default for taylorseries), while specification of output=piecewise[polynomial] provides output in terms of standard polynomial form (the default for all other methods). Use of this output form requires that the 'range' argument be used to specify to dsolve/numeric the desired range of the piecewise function.
If an array is used, it must be a vector of floats that defines the values of the independent variable at which solution values are desired. When this option is given, a 2 x 1 matrix is returned (instead of a solution procedure). The 1,1 entry is an array containing the name of the independent variable and the names of the dependent variable(s) and derivatives. The 2,1 entry is a matrix. The first column of this matrix is a copy of the output vector, that is, the values of the independent variable. The other columns are values of the dependent variable(s) and derivatives corresponding to the vector of names in 1,1. Row i of this matrix is the vector containing both the value of the independent variable and the values of the dependent variables evaluated at element i of the output array.
Note: If an Array is used, then the output is the same as for an array except the newer Array and Matrix datatypes are used in the output. For a large number of output points and/or solution components, the output data may not be directly visible because there is a default size at which the representation of an Array or Matrix is displayed instead of an inline display of the data itself. For more information, see interface(rtablesize).
'stiff'= boolean
Boolean that is used only for IVP and DAE. Setting stiff=true without selecting a method specifies that the default stiff methods (rosenbrock or rosenbrock_dae) must be used instead of the default non-stiff methods (rkf45 or rkf45_dae). When the method is also specified, a consistency check is performed to verify that the method matches with the 'stiff' value. For more information on the phenomena of stiffness, see dsolve[Stiffness].
'events'= list
'event_pre'= keyword
'event_maxiter'= integer
Arguments that are available only for real-valued IVPs or DAEs using the non-stiff and stiff methods, rkf45, ck45, rosenbrock, rkf45_dae, ck45_dae, or rosenbrock_dae, and are used to specify and control the behavior of events. These options are discussed in detail in dsolve[Events].
'range'= numeric...numeric
Values that specify the range of the independent variable over which solution values are desired.
For IVPs and DAEs this option is used only by the non-stiff and stiff methods (rkf45, ck45, rosenbrock, rkf45_dae, ck45_dae, rosenbrock_dae) and the taylorseries method. It has two purposes for procedure-type output. If 'range' is used, then the call to dsolve computes the solution over the desired range before returning, storing that solution for later calls to the returned procedure, which then compute the return values through interpolation.
Notes:
- When computing a numerical solution for a problem that has large regions where the solution is changing slowly, and small regions where the solution is changing rapidly, use of 'range' combined with the refine option of odeplot (rkf45, ck45, and rosenbrock only) allows better plotting of the details of the solution. If 'range' is not used, then the call to dsolve returns a procedure that does not store the computed solution, but rather computes the solution whenever a point is requested.
- For long-time integration problems it is suggested that 'range' not be used, as the storage of the entire solution can consume a fair bit of memory.
For the BVP method, this option is only needed when the BVP method is to be used to solve an IVP with a global error bound, as for BVP the information can be inferred from the boundary conditions in deqns. For more information, see numeric,BVP.
'abserr'= numeric or list
Numeric value that gives a limit on the absolute error tolerance for a successful step in the IVP and DAE cases, and a measure of the maximum error between the computed solution and the exact solution in the BVP case (in all but exceptional systems). It is supported by all methods except the classical methods (as all classical methods are implemented with a fixed step size and no error control). The list form of abserr is available for the rkf45, ck45, rosenbrock, and mebdfi IVP and DAE methods, and allows specification of absolute error tolerances that are different for each dependent variable, or for each solution component. More detail on this can be found in the dsolve[numeric,IVP] page and the dsolve[Error_Control] page.
'relerr'= numeric
Numeric value that gives a limit on the relative error tolerance for a successful step for an IVP or DAE. This option works in conjunction with 'abserr'. It is supported by all methods except classical, taylorseries, and bvp.
Defaults for 'abserr' and 'relerr' are specific to each method, but it should be noted that the default error tolerance for the non-stiff methods rkf45 and ck45 is now fixed (rather than controlled by the value of Digits).
For a detailed discussion about error control, see dsolve[Error_Control].
'known'= name or list of names
Provides a list of functions that should be considered as known when examining the system. This allows for specification of a system containing user-defined functions. Each user defined function must always return numerical values when given numeric input, and must return unevaluated when given name input. For an illustration of the use of this option, see the Examples section below.
Note: In most cases, use of a system with a user-defined function prevents the use of evalhf within dsolve, so obtaining a solution using this facility will run noticeably slower than specification of the system in terms of known functions (when it is possible to express the system in this way).
All methods except taylorseries, rosenbrock, and rosenbrock_dae support this option directly, but for these two methods additional information is required. The key problem is that these two methods require information on the derivatives of the input ODE system, which is not directly available when functions are procedure defined. It is still possible to use these methods with 'known', as long as the derivatives of the function are defined via use of one or more `diff/` rules. The rosenbrock method requires only a single derivative, but the taylorseries method requires derivatives to arbitrary order, so with taylorseries it is useful only if some finite derivative of the 'known' functions can be represented as a function instead of a procedure, in which case, all derivatives up to that function derivative must have a diff rule defined. The rosenbrock_dae method may require more than a single derivative, as the DAE preprocessing may need to differentiate the input system.
For taylorseries, rosenbrock, and rosenbrock_dae use of multiple argument functions, or functions that also depend on the dependent variables or derivatives of the problem are tricky, and have some limitations. The details of this are deferred to the respective help pages.
'optimize'= boolean
Boolean that tells dsolve to optimize the input system for computation. This adds an up-front cost, with the possible reward that the core computation itself will be faster. This is false by default for IVP and BVP problems, and true by default for DAE problems (which typically derive the most benefit from optimization). Systems containing conditional evaluation functions (such as piecewise) should not use optimize=true, as this will sometimes result in a slower computation. This option is available only for system defined problems, and for all methods except taylorseries (which uses its own form of optimization).
All other options, including options for specification of procedures for the evaluation of the ODE system, and options specific to IVP or BVP are discussed in dsolve[numeric,IVP] and dsolve[numeric,BVP], or the help pages specific to each method.
Notes
When Digits<=evalhf⁡Digits, numerical solutions are computed in double precision (hardware floats). This is the default, as on initialization Maple sets Digits to 10. Explicit setting of Digits>evalhf⁡Digits causes dsolve to compute numerical solutions using Maple floats instead of hardware floats. The precision of a computation is fixed in the call to dsolve, so for procedure outputs, once the procedure is created, further changes to the setting of Digits should have no effect on the computed solution.
Numerical solutions can be used in combination with fsolve, but the Digits setting on the call to fsolve must be lower (usually by 5 or more) than the precision for the dsolve procedure (because fsolve requires a higher accuracy on the residual than the current Digits setting). The rkf45, ck45, rosenbrock, and taylorseries methods, when computing solution values with an interpolant, work around this problem by computing the interpolant with the current setting of Digits when called from fsolve. An example of the problem is presented in the last example below.
Results can be plotted by using the odeplot function in the plots package.
ODE systems can also be interactively solved and plotted with the dsolve[interactive] Maplet interface.
Solve a system of ODEs
This system describes simple harmonic motion with no friction.
dsys≔diff⁡x⁡t,t=y⁡t,diff⁡y⁡t,t=−x⁡t,x⁡0=1,y⁡0=0:
dsn1≔dsolve⁡dsys,numeric:
dsn1⁡1
t=1.,x⁡t=0.540302331778567,y⁡t=−0.841471101155307
If you specify the set of dependent variables as a list, those variables appear in the same order in the output.
dsn2≔dsolve⁡dsys,numeric,y⁡t,x⁡t:
dsn2⁡1
t=1.,y⁡t=−0.841471101155307,x⁡t=0.540302331778567
Straightforward second order IVP:
deq1≔t+12⁢diff⁡y⁡t,t,t+t+1⁢diff⁡y⁡t,t+t+12−0.25⁢y⁡t=0
deq1≔t+12⁢ⅆ2ⅆt2y⁡t+t+1⁢ⅆⅆty⁡t+t+12−0.25⁢y⁡t=0
ic1≔y⁡0=0.6713967071418030,D⁡y⁡0=0.09540051444747446:
dsol1≔dsolve⁡deq1,ic1,numeric,range=0..1
dsol1 ≔ procx_rkf45...end proc
dsol1⁡0
t=0.,y⁡t=0.671396707141803,ⅆⅆty⁡t=0.0954005144474745
dsol1⁡0.5
t=0.5,y⁡t=0.649838013736792,ⅆⅆty⁡t=−0.170529562812705
dsol1⁡1
t=1.,y⁡t=0.513016090105034,ⅆⅆty⁡t=−0.363039820324515
Same second order IVP with operator output:
dsol1b≔dsolve⁡deq1,ic1,numeric,output=operator
dsol1b≔t=proct...end proc,y=proct...end proc,D⁡y=proct...end proc
dsol1b⁡0
t=0.,y⁡0=0.671396707141803,D⁡y⁡0=0.0954005144474745
dsol1b⁡0.5
t=0.5,y⁡0.5=0.649838013736792,D⁡y⁡0.5=−0.170529562812705
dsol1b⁡1
t=1.,y⁡1=0.513016090105034,D⁡y⁡1=−0.363039820324514
First order variable coefficient IVP:
dsol2≔dsolve⁡diff⁡y⁡x,x=y⁡x⁢cos⁡x,y⁡0=1,numeric,output=listprocedure
dsol2≔x=procx...end proc,y⁡x=procx...end proc
fy2≔eval⁡y⁡x,dsol2
fy2 ≔ procx...end proc
fy2⁡0,fy2⁡π4,fy2⁡π2,fy2⁡3⁢π2,fy2⁡π
1.,2.02811545714854,2.71828201727134,0.367879858344611,0.999999785849056
First order IVP with piecewise output:
dsol2p≔dsolve⁡diff⁡y⁡x,x=y⁡x⁢cos⁡x,y⁡0=1,numeric,output=piecewise,range=0..1
dsol2p≔x=x,y⁡x=undefinedx<0.0.999614858592639+1.02775126621593⁢x+0.499406940429015⁢x−0.02776212315712852−0.0145000342870509⁢x−0.02776212315712853−0.133835263967526⁢x−0.02776212315712854x≤0.05552424631425700860.989101552131176+1.14719850211587⁢x+0.481091117076617⁢x−0.1490132628421082−0.0904918748856021⁢x−0.1490132628421083−0.172313265589485⁢x−0.1490132628421084x≤0.2425022793699588840.948397466098112+1.31597042278255⁢x+0.387428683654558⁢x−0.3400573880490192−0.244571952705452⁢x−0.3400573880490193−0.220713724578106⁢x−0.3400573880490194x≤0.4376124967280781710.897425216810995+1.43475629979104⁢x+0.183288214102761⁢x−0.5413675335024242−0.432634746442101⁢x−0.5413675335024243−0.225886957923084⁢x−0.5413675335024244x≤0.6451225702767702240.890827299362533+1.44839243200584⁢x−0.133255972479318⁢x−0.7436285408746302−0.598745668466927⁢x−0.7436285408746303−0.155195537009859⁢x−0.7436285408746304x≤0.8421345114724896330.981744432894418+1.34133782273687⁢x−0.477242160985440⁢x−0.9210672557362472−0.673533875002395⁢x−0.9210672557362473−0.0140365582324127⁢x−0.9210672557362474x≤1.undefinedotherwise
Stiff first-order nonlinear problem:
dsys3≔diff⁡x⁡t,t=x⁡t2−t6,x⁡0=−0.1
dsys3≔ⅆⅆtx⁡t=x⁡t2−t6,x⁡0=−0.1
dsol3≔dsolve⁡dsys3,type=numeric,stiff=true,output=Array⁡0,0.25,0.5,0.75,1
dsol3≔tx⁡t0.−0.1000000000000000.250000000000000−0.09756970475011670.500000000000000−0.09634073253273390.750000000000000−0.1117447901123991.−0.229279026435938
Second-order linear BVP:
dsol4≔dsolve⁡diff⁡y⁡x,x,x=3⁢y⁡x,y⁡0=1,y⁡2=1.2,numeric
dsol4 ≔ procx_bvp...end proc
dsol4⁡0
x=0.,y⁡x=1.,ⅆⅆxy⁡x=−1.60520425463385
dsol4⁡0.5
x=0.5,y⁡x=0.492273627127114,ⅆⅆxy⁡x=−0.551072761778113
dsol4⁡1
x=1.,y⁡x=0.377412857949490,ⅆⅆxy⁡x=0.0632677085684644
dsol4⁡2
x=2.,y⁡x=1.20000000000000,ⅆⅆxy⁡x=1.97400119243652
Second-order nonlinear BVP:
dsys5≔diff⁡y⁡x,x,x+abs⁡y⁡x=0,y⁡1=−1,D⁡y⁡0=1
dsys5≔ⅆ2ⅆx2y⁡x+y⁡x=0,y⁡1=−1,D⁡y⁡0=1
dsol5≔dsolve⁡dsys5,numeric,output=Array⁡0,0.25,0.5,0.75,1
dsol5≔xy⁡xⅆⅆxy⁡x0.−1.409648429587521.000000000000000.25−1.201317539681560.6753185441771470.5−1.068460865510630.3930647861123990.75−1.002731527130150.1355057948096101.−1.00000000000000−0.113539882361526
DAE system for the simple 2-D pendulum:
dsys6≔x⁡t2+y⁡t2=1,diff⁡x⁡t,t,t=−2⁢λ⁡t⁢x⁡t,diff⁡y⁡t,t,t=−2⁢λ⁡t⁢y⁡t−π2,x⁡0=0,y⁡0=−1,D⁡x⁡0=110,D⁡y⁡0=0
dsys6≔x⁡t2+y⁡t2=1,ⅆ2ⅆt2x⁡t=−2⁢λ⁡t⁢x⁡t,ⅆ2ⅆt2y⁡t=−2⁢λ⁡t⁢y⁡t−π2,x⁡0=0,y⁡0=−1,D⁡x⁡0=110,D⁡y⁡0=0
dsol6≔dsolve⁡dsys6,numeric
dsol6 ≔ procx_rkf45_dae...end proc
dsol6⁡1
t=1.,λ⁡t=4.93980221499583,x⁡t=6.33398801850387×10−6,ⅆⅆtx⁡t=−0.100000139089421,y⁡t=−0.999999999836514,ⅆⅆty⁡t=−6.24735803827156×10−7
Example of use with fsolve:
dsys7≔diff⁡y⁡x,x=2⁢y⁡x,y⁡0=1
dsys7≔ⅆⅆxy⁡x=2⁢y⁡x,y⁡0=1
Digits≔20
dsol7≔dsolve⁡dsys7,numeric,method=gear,abserr=1.×10−16,relerr=1.×10−16,output=listprocedure
dsol7≔x=procx...end proc,y⁡x=procx...end proc
yproc≔rhs⁡dsol72
yproc ≔ procx...end proc
Attempt with same Digits setting, returns unevaluated.
fsolve(yproc(x)=1e5, x=4..7);
fsolve⁡yproc⁡x=100000.,x,4..7
The same command with Digits set to 15 succeeds.
Digits≔15
5.75646273248511
Digits≔10
Use of known option:
f := proc(x) local t; if not type(evalf(x),'numeric') then 'procname'(x); else evalf(Int(exp(-t^2/10),t=0..x)); end if; end proc;
f ≔ procxlocalt;ifnottype⁡evalf⁡x,'numeric'then'procname'⁡xelseevalf⁡Int⁡exp⁡−1/10*t^2,t=0..xend ifend proc
Note that the procedure is set up for an unevaluated return if the input is not numeric.
dsys8≔diff⁡y⁡x,x=y⁡x+f⁡x,y⁡0=0
dsys8≔ⅆⅆxy⁡x=y⁡x+f⁡x,y⁡0=0
dsolve⁡dsys8,numeric
Error, (in `dsolve/numeric/DAE/make_proc`) number of unknown functions and equations must match, got 2 functions {f, y}, and 1 equations
This failed because 'known' was not specified. Here, it is specified.
dsol8≔dsolve⁡dsys8,numeric,known=f
dsol8 ≔ procx_rkf45...end proc
dsol8⁡1
x=1.,y⁡x=0.708149294799617
dsol8⁡2
x=2.,y⁡x=4.19167993111610
A method that requires derivatives,
dsolve⁡dsys8,numeric,method=rosenbrock,known=f
Error, (in `dsolve/numeric/checkknown`) use of 'rosenbrock' method with procedure defined function f requires a differentiation rule `diff/f`
which failed because a differentiation rule is required. We then clear the remember table of diff and define this as:
forget⁡diff
`diff/f` := proc() # Chain rule diff(args[1],args[2])*exp(-(args[1])^2/10); end proc:
which works as follows:
diff⁡f⁡sin⁡x,x
cos⁡x⁢ⅇ−sin⁡x210
diff⁡Int⁡exp⁡−t210,t=0..sin⁡x,x
and now:
dsol8b≔dsolve⁡dsys8,numeric,method=rosenbrock,known=f
dsol8b ≔ procx_rosenbrock...end proc
dsol8b⁡1
x=1.,y⁡x=0.708149269188472
dsol8b⁡2
x=2.,y⁡x=4.19168036703643
Ascher, U.; Mattheij, R.; and Russell, R. "Numerical Solution of Boundary Value Problems for Ordinary Differential Equations." SIAM Classics in Applied Mathematics. Vol. 13. (1995).
Ascher, U., and Petzold, L. "Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations." SIAM, Philadelphia. 1998.
Bailey, P.B.; Shampine, L.F.; and Waltman, P.E. Nonlinear Two Point Boundary Value Problems. New York: Academic Press, 1968.
Boyce, W.E., and DiPrima, R.C. Elementary Differential Equations and Boundary Value Problems. New York: John Wiley & Sons, 1997.
Cash, J.R. "The Integration of Stiff IVP in ODE Using Modified Extended BDF." Computers and Mathematics with Applications. Vol. 9. (1983): 645-657.
Gear, C.W. Numerical Initial Value Problems in Ordinary Differential Equations. Prentice-Hall, 1971.
Hindmarsh, Alan C.; Stepleman, R.S.; et al, eds. Odepack, a Systemized Collection of ODE Solvers. Amsterdam: North-Holland, 1983.
Hubbard J.H., and West, B.H. Differential Equations: A Dynamical Systems Approach. New York: Springer, 1990. Part I. One Dimensional Equations.
Hull, T.E.; Enright, W.H.; Fellen, B.M.; and Sedgwick, A.E. "Comparing Numerical Methods for Ordinary Differential Equations." SIAM J. Numer. Anal. Vol. 9. (1972): 603-637.
Shampine, L.F., and Corless, R.M. "Initial Value Problems for ODEs in Problem Solving Environments." J. Comp. Appl. Math. Vol. 125(1-2). (2000): 31-40.
See Also
Digits
dsolve/DAE_extension
dsolve/Error_Control
dsolve[ck45]
dsolve[classical]
dsolve[dverk78]
dsolve[Events]
dsolve[gear]
dsolve[interactive]
dsolve[lsode]
dsolve[maxfun]
dsolve[mebdfi]
dsolve[numeric,BVP]
dsolve[numeric,DAE]
dsolve[numeric,interactive]
dsolve[numeric,IVP]
dsolve[rkf45]
dsolve[rosenbrock]
dsolve[Stiffness]
dsolve[taylorseries]
fsolve
interface
plots[odeplot]
Download Help Document