The dsolve command with the numeric or type=numeric option and a real-valued differential-algebraic initial value problem (DAE IVP) finds a numerical solution for the DAE IVP. If the optional equation method=numericmethod is provided (where numericmethod is one of rkf45_dae, ck45_dae, rosenbrock_dae, or mebdfi), dsolve uses that method to obtain the numerical solution.

In most cases dsolve is able to detect if the problem is a DAE system, as opposed to an ODE system, namely the cases in which pure algebraic equations in the dependent variables are present. If the input is a DAE system containing no purely algebraic equations, the method must be included to specify that the system is a DAE system.

Constrained mechanical systems often give rise to DAE problems. (See the pendulum example below.)

The return value of dsolve and the following high level or common options are discussed in dsolve[numeric] and dsolve[Error_Control].

The exception is that for the DAE solvers the absolute error tolerance can be specified as a per-component list. For expected behavior, the variables of the problem must also be specified as a list, and the entries of 'abserr' must correspond 1-1 to the variables in the list, or to the variables in the system converted to first order using the order of the variables in the list. For more information, see the ?dsolve[numeric,IVP] help page.

The default DAE IVP method is a modified Runge-Kutta Fehlberg method, which uses a base order 4-5 method, but has been modified to find solutions for DAE problems. The default stiff method is a Rosenbrock method, which uses a base order 3-4 method. For a description of the modifications done to these methods in extending them to DAE solution, see the ?DAE_extension help page. The other method available for DAE IVP is the dsolve[mebdfi] method, which is short for Modified Extended Backward-Differentiation Formula Implicit method.

In general, the DAE IVP solvers are very similar to the standard differential IVP solvers, so this page is primarily concerned with outlining the differences between them.

The DAE solvers are currently restricted to finding solutions for real-valued problems.

For use of any of the methods, the specified initial conditions must satisfy all hidden constraints of the system (that is, they must be a consistent set of initial conditions with respect to the DAE). In the event that they do not, an error results, and information is provided on the unsatisfied condition.

In some cases, it may be necessary to use fsolve to compute consistent initial conditions for the problem.

If the ?DAE_extension methods are in use, the differential option is set to true, and the system is sufficiently linear in the algebraic variables (i.e., variables which have no derivatives appearing in the input system), it is possible to skip initial conditions for those variables. If the initial conditions are skipped when they are required, an error will be produced.

The following options are also available for some or all of the DAE methods:

'minstep', 'maxstep', and 'initstep'

These options are discussed in dsolve[Error_Control], and are only available for the dsolve[mebdfi] method.

'startinit'= boolean

This option controls the behavior of the numerical integration with respect to consecutive calculations. This option is described in dsolve[numeric,IVP].

For the methods, rkf45_dae, ck45_dae, and rosenbrock_dae, when a 'range' has been specified, the solution is recomputed only when the new value of the independent variable is not between the initial point and any other previously requested solution point. This has the effect of never reversing the direction of integration, and making evaluation of the solution for an already computed interval quite inexpensive (the solution values are obtained by interpolation). The storage of the solution can also be enabled by using the storage argument. The startinit parameter also forces these methods to recompute the solution whenever a solution value is desired.

'events'= list

'event_pre'= keyword

'event_maxiter'= integer

These options are available only for the rkf45_dae, ck45_dae, and rosenbrock_dae methods, and are used to specify and control the behavior of events. These are the same as for standard IVP problems. For a description, see the ?dsolve[numeric,IVP] and ?dsolve[Events] help pages.

'projection', 'differential', and 'implicit'

The 'projection', 'differential' and 'implicit' options are specific to the extension methods, so they are discussed there.

In addition to the computation of solution values for the given input problem, the procedure returns (that is, output=procedurelist, listprocedure or operator) provide additional interactive features, including the ability to specify parameters. Information on these features is provided on the dsolve[numeric,interactive] page. The exception is that the mebdfi method cannot work with parametrized problems.

$\mathrm{VectorCalculus}\left[\mathrm{SetCoordinates}\right]\left(\mathrm{cartesian}\right)\:$

$\mathrm{VectorCalculus}\left[\mathrm{BasisFormat}\right]\left(\mathrm{false}\right)\:$

$r≔\mathrm{VectorCalculus}\left[\mathrm{Vector}\right]\left(\left[x\left(t\right)\,y\left(t\right)\right]\right)$

$r≔\left[\begin{array}{c}x\left(t\right)\\ y\left(t\right)\end{array}\right]$

$v≔\mathrm{VectorCalculus}\left[\mathrm{diff}\right]\left(r\,t\right)$

$v≔\left[\begin{array}{c}\frac{\ⅆ}{\ⅆt}x\left(t\right)\\ \frac{\ⅆ}{\ⅆt}y\left(t\right)\end{array}\right]$

$\mathrm{con}≔\mathrm{VectorCalculus}\left[\mathrm{DotProduct}\right]\left(r\,r\right)-l^2$

$\mathrm{con}≔{x\left(t\right)}^2+{y\left(t\right)}^2-l^2$

$T≔\frac{m}{2}\mathrm{VectorCalculus}\left[\mathrm{DotProduct}\right]\left(v\,v\right)$

$T≔\frac{m\left({\left(\frac{\ⅆ}{\ⅆt}x\left(t\right)\right)}^2+{\left(\frac{\ⅆ}{\ⅆt}y\left(t\right)\right)}^2\right)}{2}$

$V≔mgy\left(t\right)$

$V≔mgy\left(t\right)$

$L≔T-V$

$L≔\frac{m\left({\left(\frac{\ⅆ}{\ⅆt}x\left(t\right)\right)}^2+{\left(\frac{\ⅆ}{\ⅆt}y\left(t\right)\right)}^2\right)}{2}-mgy\left(t\right)$

$\mathrm{ML}≔L-\mathrm{\lambda }\left(t\right)\mathrm{con}$

$\mathrm{ML}≔\frac{m\left({\left(\frac{\ⅆ}{\ⅆt}x\left(t\right)\right)}^2+{\left(\frac{\ⅆ}{\ⅆt}y\left(t\right)\right)}^2\right)}{2}-mgy\left(t\right)-\mathrm{\lambda }\left(t\right)\left({x\left(t\right)}^2+{y\left(t\right)}^2-l^2\right)$

$Q≔\mathrm{remove}\left(\mathrm{has}\,\mathrm{VariationalCalculus}:-\mathrm{EulerLagrange}\left(\mathrm{ML}\,t\,\left[x\left(t\right)\,y\left(t\right)\,\mathrm{\lambda }\left(t\right)\right]\right)\,K\right)$

$Q≔\left\{-2x\left(t\right)\mathrm{\lambda }\left(t\right)-m\left(\frac{{\ⅆ}^2}{\ⅆt^2}x\left(t\right)\right)\,l^2-{x\left(t\right)}^2-{y\left(t\right)}^2\,-mg-2y\left(t\right)\mathrm{\lambda }\left(t\right)-m\left(\frac{{\ⅆ}^2}{\ⅆt^2}y\left(t\right)\right)\right\}$

$\mathrm{EQx}≔\mathrm{isolate}\left(\mathrm{select}\left(\mathrm{has}\,Q\,\mathrm{diff}\left(x\left(t\right)\,t\,t\right)\right)\left[1\right]\,\mathrm{diff}\left(x\left(t\right)\,t\,t\right)\right)$

$\mathrm{EQx}≔\frac{{\ⅆ}^2}{\ⅆt^2}x\left(t\right)=-\frac{2x\left(t\right)\mathrm{\lambda }\left(t\right)}{m}$

$\mathrm{EQy}≔\mathrm{isolate}\left(\mathrm{select}\left(\mathrm{has}\,Q\,\mathrm{diff}\left(y\left(t\right)\,t\,t\right)\right)\left[1\right]\,\mathrm{diff}\left(y\left(t\right)\,t\,t\right)\right)$

$\mathrm{EQy}≔\frac{{\ⅆ}^2}{\ⅆt^2}y\left(t\right)=-\frac{mg+2y\left(t\right)\mathrm{\lambda }\left(t\right)}{m}$

$\mathrm{EQc}≔\mathrm{remove}\left(\mathrm{has}\,Q\,\mathrm{diff}\right)\left[1\right]$

$\mathrm{EQc}≔l^2-{x\left(t\right)}^2-{y\left(t\right)}^2$

$\mathrm{Dcon}≔\mathrm{diff}\left(\mathrm{EQc}\,t\right)$

$\mathrm{Dcon}≔-2x\left(t\right)\left(\frac{\ⅆ}{\ⅆt}x\left(t\right)\right)-2y\left(t\right)\left(\frac{\ⅆ}{\ⅆt}y\left(t\right)\right)$

$\mathrm{DDcon}≔\mathrm{eval}\left(\mathrm{diff}\left(\mathrm{Dcon}\,t\right)\,\left\{\mathrm{EQx}\,\mathrm{EQy}\right\}\right)$

$\mathrm{DDcon}≔-2{\left(\frac{\ⅆ}{\ⅆt}x\left(t\right)\right)}^2+\frac{4{x\left(t\right)}^2\mathrm{\lambda }\left(t\right)}{m}-2{\left(\frac{\ⅆ}{\ⅆt}y\left(t\right)\right)}^2+\frac{2y\left(t\right)\left(mg+2y\left(t\right)\mathrm{\lambda }\left(t\right)\right)}{m}$

$\mathrm{sys}≔\left\{y\left(0\right)=-l\,\mathrm{D}\left(x\right)\left(0\right)=\mathrm{vx}\right\}\mathbf{union}\mathrm{eval}\left(\mathrm{convert}\left(\left\{\mathrm{DDcon}\,\mathrm{Dcon}\,\mathrm{EQc}\right\}\,\mathrm{D}\right)\,t=0\right)$

$\mathrm{sys}≔\left\{-2x\left(0\right)\mathrm{D}\left(x\right)\left(0\right)-2y\left(0\right)\mathrm{D}\left(y\right)\left(0\right)\,l^2-{x\left(0\right)}^2-{y\left(0\right)}^2\,-2{\mathrm{D}\left(x\right)\left(0\right)}^2+\frac{4{x\left(0\right)}^2\mathrm{\lambda }\left(0\right)}{m}-2{\mathrm{D}\left(y\right)\left(0\right)}^2+\frac{2y\left(0\right)\left(mg+2y\left(0\right)\mathrm{\lambda }\left(0\right)\right)}{m}\,y\left(0\right)=-l\,\mathrm{D}\left(x\right)\left(0\right)=\mathrm{vx}\right\}$

$\mathrm{ini}≔\left[\mathrm{solve}\left(\mathrm{sys}\,\left\{\mathrm{\lambda }\left(0\right)\,x\left(0\right)\,y\left(0\right)\,\mathrm{D}\left(x\right)\left(0\right)\,\mathrm{D}\left(y\right)\left(0\right)\right\}\right)\right]\left[1\right]$

$\mathrm{ini}≔\left\{\mathrm{\lambda }\left(0\right)=\frac{m\left(gl+{\mathrm{vx}}^2\right)}{2l^2}\,x\left(0\right)=0\,y\left(0\right)=-l\,\mathrm{D}\left(x\right)\left(0\right)=\mathrm{vx}\,\mathrm{D}\left(y\right)\left(0\right)=0\right\}$

$\mathrm{dsys}≔\mathrm{eval}\left(\left\{\mathrm{EQc}\,\mathrm{EQx}\,\mathrm{EQy}\right\}\,\left\{g=9.8\,l=1\,m=1\,\mathrm{vx}=\frac{1}{10}\right\}\right)$

$\mathrm{dsys}≔\left\{1-{x\left(t\right)}^2-{y\left(t\right)}^2\,\frac{{\ⅆ}^2}{\ⅆt^2}x\left(t\right)=-2x\left(t\right)\mathrm{\lambda }\left(t\right)\,\frac{{\ⅆ}^2}{\ⅆt^2}y\left(t\right)=-9.8-2y\left(t\right)\mathrm{\lambda }\left(t\right)\right\}$

$\mathrm{dini}≔\mathrm{eval}\left(\mathrm{ini}\,\left\{g=9.8\,l=1\,m=1\,\mathrm{vx}=\frac{1}{10}\right\}\right)$

$\mathrm{dini}≔\left\{\mathrm{\lambda }\left(0\right)=4.905000000\,x\left(0\right)=0\,y\left(0\right)=\mathrm{-1}\,\mathrm{D}\left(x\right)\left(0\right)=\frac{1}{10}\,\mathrm{D}\left(y\right)\left(0\right)=0\right\}$

$\mathrm{dsol1}≔\mathrm{dsolve}\left(\mathrm{dsys}\mathbf{union}\mathrm{dini}\,\mathrm{numeric}\right)$

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

$\mathrm{dsol1}\left(\frac{1}{2}\right)$

$\left[t=0.500000000000000\,\mathrm{\lambda }\left(t\right)=4.89750023467855\,x\left(t\right)=0.0319392699415791\,\frac{\ⅆ}{\ⅆt}x\left(t\right)=0.000564539406239432\,y\left(t\right)=\mathrm{-0.999489811490763}\,\frac{\ⅆ}{\ⅆt}y\left(t\right)=0.0000180384378593043\right]$

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

$\left[t=1.\,\mathrm{\lambda }\left(t\right)=4.90499905723651\,x\left(t\right)=0.000360891517794992\,\frac{\ⅆ}{\ⅆt}x\left(t\right)=\mathrm{-0.0999937504894365}\,y\left(t\right)=\mathrm{-0.999999934755133}\,\frac{\ⅆ}{\ⅆt}y\left(t\right)=\mathrm{-0.0000360794101890703}\right]$

$\mathrm{dsol2}≔\mathrm{dsolve}\left(\mathrm{dsys}\mathbf{union}\mathrm{dini}\,\mathrm{numeric}\,\mathrm{method}=\mathrm{rosenbrock\_dae}\right)$

$\mathrm{dsol2}≔\mathbf{proc}\left(\mathrm{x\_rosenbrock\_dae}\right)...\mathbf{end\; proc}$

$\mathrm{dsol2}\left(\frac{1}{2}\right)$

$\left[t=0.500000000000000\,\mathrm{\lambda }\left(t\right)=4.89750019118586\,x\left(t\right)=0.0319392364763694\,\frac{\ⅆ}{\ⅆt}x\left(t\right)=0.000564626004686083\,y\left(t\right)=\mathrm{-0.999489813543334}\,\frac{\ⅆ}{\ⅆt}y\left(t\right)=0.0000180482728717227\right]$

$\mathrm{dsol2}\left(1\right)$

$\left[t=1.\,\mathrm{\lambda }\left(t\right)=4.90499906207811\,x\left(t\right)=0.000360968463263150\,\frac{\ⅆ}{\ⅆt}x\left(t\right)=\mathrm{-0.0999936779452163}\,y\left(t\right)=\mathrm{-0.999999936691318}\,\frac{\ⅆ}{\ⅆt}y\left(t\right)=\mathrm{-0.0000360897908606990}\right]$

$\mathrm{dsol3}≔\mathrm{dsolve}\left(\mathrm{dsys}\mathbf{union}\mathrm{dini}\,\mathrm{numeric}\,\mathrm{method}=\mathrm{mebdfi}\right)$

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

$\mathrm{dsol3}\left(\frac{1}{2}\right)$

$\left[t=0.500000000000000\,\mathrm{\lambda }\left(t\right)=4.89749908804820855\,x\left(t\right)=0.0319389748677520596\,\frac{\ⅆ}{\ⅆt}x\left(t\right)=0.000565199414553933343\,y\left(t\right)=\mathrm{-0.999489820801930273}\,\frac{\ⅆ}{\ⅆt}y\left(t\right)=0.0000180581765179766362\right]$

$\mathrm{dsol3}\left(1\right)$

$\left[t=1.\,\mathrm{\lambda }\left(t\right)=4.90499581413915475\,x\left(t\right)=0.000361181747066772996\,\frac{\ⅆ}{\ⅆt}x\left(t\right)=\mathrm{-0.0999926673290293111}\,y\left(t\right)=\mathrm{-0.99999993019115307}\,\frac{\ⅆ}{\ⅆt}y\left(t\right)=\mathrm{-0.0000360957473722881540}\right]$

$\mathrm{dsysd}≔\left\{\mathrm{diff}\left(\mathrm{y1}\left(t\right)\,t\,t\right)+9.8+2\mathrm{\λ1}\left(t\right)\mathrm{y1}\left(t\right)+2\mathrm{\λ2}\left(t\right)\left(\mathrm{y1}\left(t\right)-\mathrm{y2}\left(t\right)\right)\,\mathrm{diff}\left(\mathrm{x1}\left(t\right)\,t\,t\right)+2\mathrm{\λ1}\left(t\right)\mathrm{x1}\left(t\right)+2\mathrm{\λ2}\left(t\right)\left(\mathrm{x1}\left(t\right)-\mathrm{x2}\left(t\right)\right)\,{\left(\mathrm{x1}\left(t\right)-\mathrm{x2}\left(t\right)\right)}^2+{\left(\mathrm{y1}\left(t\right)-\mathrm{y2}\left(t\right)\right)}^2-\frac{1}{4}\,{\mathrm{x1}\left(t\right)}^2+{\mathrm{y1}\left(t\right)}^2-1\,\mathrm{diff}\left(\mathrm{y2}\left(t\right)\,t\,t\right)+9.8-2\mathrm{\λ2}\left(t\right)\left(\mathrm{y1}\left(t\right)-\mathrm{y2}\left(t\right)\right)\,\mathrm{diff}\left(\mathrm{x2}\left(t\right)\,t\,t\right)-2\mathrm{\λ2}\left(t\right)\left(\mathrm{x1}\left(t\right)-\mathrm{x2}\left(t\right)\right)\right\}$

$\mathrm{dsysd}≔\left\{\frac{{\ⅆ}^2}{\ⅆt^2}\mathrm{x2}\left(t\right)-2\mathrm{\λ2}\left(t\right)\left(\mathrm{x1}\left(t\right)-\mathrm{x2}\left(t\right)\right)\,{\left(\mathrm{x1}\left(t\right)-\mathrm{x2}\left(t\right)\right)}^2+{\left(\mathrm{y1}\left(t\right)-\mathrm{y2}\left(t\right)\right)}^2-\frac{1}{4}\,{\mathrm{x1}\left(t\right)}^2+{\mathrm{y1}\left(t\right)}^2-1\,\frac{{\ⅆ}^2}{\ⅆt^2}\mathrm{x1}\left(t\right)+2\mathrm{\λ1}\left(t\right)\mathrm{x1}\left(t\right)+2\mathrm{\λ2}\left(t\right)\left(\mathrm{x1}\left(t\right)-\mathrm{x2}\left(t\right)\right)\,\frac{{\ⅆ}^2}{\ⅆt^2}\mathrm{y2}\left(t\right)+9.8-2\mathrm{\λ2}\left(t\right)\left(\mathrm{y1}\left(t\right)-\mathrm{y2}\left(t\right)\right)\,\frac{{\ⅆ}^2}{\ⅆt^2}\mathrm{y1}\left(t\right)+9.8+2\mathrm{\λ1}\left(t\right)\mathrm{y1}\left(t\right)+2\mathrm{\λ2}\left(t\right)\left(\mathrm{y1}\left(t\right)-\mathrm{y2}\left(t\right)\right)\right\}$

$\mathrm{ics}≔\left\{\mathrm{x1}\left(0\right)=0\,\mathrm{x2}\left(0\right)=0\,\mathrm{y1}\left(0\right)=-1\,\mathrm{y2}\left(0\right)=-\frac{3}{2}\,\mathrm{D}\left(\mathrm{x1}\right)\left(0\right)=-3\,\mathrm{D}\left(\mathrm{x2}\right)\left(0\right)=4\,\mathrm{D}\left(\mathrm{y1}\right)\left(0\right)=0\,\mathrm{D}\left(\mathrm{y2}\right)\left(0\right)=0\right\}$

$\mathrm{ics}≔\left\{\mathrm{x1}\left(0\right)=0\,\mathrm{x2}\left(0\right)=0\,\mathrm{y1}\left(0\right)=\mathrm{-1}\,\mathrm{y2}\left(0\right)=-\frac{3}{2}\,\mathrm{D}\left(\mathrm{x1}\right)\left(0\right)=\mathrm{-3}\,\mathrm{D}\left(\mathrm{x2}\right)\left(0\right)=4\,\mathrm{D}\left(\mathrm{y1}\right)\left(0\right)=0\,\mathrm{D}\left(\mathrm{y2}\right)\left(0\right)=0\right\}$

$\mathrm{dsold}≔\mathrm{dsolve}\left(\mathrm{dsysd}\mathbf{union}\mathrm{ics}\,\mathrm{numeric}\right)$

$\mathrm{dsold}≔\mathbf{proc}\left(\mathrm{x\_rkf45\_dae}\right)...\mathbf{end\; proc}$

$\mathrm{plots}\left[\mathrm{odeplot}\right]\left(\mathrm{dsold}\,\left[\mathrm{x2}\left(t\right)\,\mathrm{y2}\left(t\right)\right]\,0..5\,\mathrm{numpoints}=1000\right)$