The dsolve DAE extension methods are standard numerical methods designed for ODE IVP that have been extended to the solution of DAE problems.

Specifically, the solvers have been extended to the solution of regular ODE and limited index 1 DAE (that can be isolated for a required dependent variable). The system is reduced to this specific index 1 system via index reduction.

Note: The index can be thought of as how far the DAE system is from being an ODE system.

The DAE extension methods then use the modified solver at each step. At the end of each step, the solution is projected back so that it satisfies any additional constraints of the problem. These additional constraints arise in the process of index reduction, and will only be present for systems for which index reduction is required.

dsolve has three extension methods, rkf45_dae, which is an extension of the rkf45 method, ck45_dae, which is an extension of the ck45 method, and rosenbrock_dae, which is an extension of the rosenbrock method.

As extensions, these methods retain many of the features of the original solvers, for example, events and range solutions. This page is concerned only with highlighting the differences between these methods, and any extension specific options.

One difference is that, by default, optimization is active for the extension methods, as the system is usually quite complex after relevant postprocessing has occurred.

For setup of the problem, it is required that the initial conditions are consistent, that is, that they satisfy any hidden constraints of the problem. In any cases where the constraints are not sufficiently satisfied, an error results, and the violated condition(s) are printed, allowing adjustment of the initial data to satisfy them.

In addition, if differential is false, it is possible to automatically compute values for any of the index 1 variables of the problem. Initial conditions for these variables need not be provided.

There are two major differences between the capabilities of the extension solvers and their original counterparts: the extension solvers can only be used with real valued problems and the initial values for the problem cannot be changed interactively.

Three options are specific to the DAE extension solvers:

'projection'= boolean

This option specifies whether, after each step, the solution is projected back onto the set of constraints that are not in direct use by the solver. This is true by default. If the DAE IVP is highly stable, and speed is an issue, setting this to false can provide a reasonable solution in less time.

'differential'= boolean

This option specifies that index reduction be performed until the system is fully differential (that is, does not take advantage of the index 1 capabilities of the modified solvers). This is false by default. The use of the index 1 capabilities of the solvers often provides a more stable and less stiff system than the full index reduction, and is usually faster.

'implicit'= boolean

This option specifies that parts of the extended system can be left in unsolved (implicit) form until they are evaluated. The solver determines which parts of the system must be put into solved form, in order to render the DAE problem numerically solvable, and which parts can be left in implicit form. This is essentially a smart form of the implicit option used for IVP problems (see dsolve[numeric,IVP]). This is a helpful option. It reduces the complication of the system derived quantities that must be computed to perform the numerical integration. The cost, though, is a linear solve every time the expressions must be evaluated. For problems in which the system is nearly in solved form with respect to the leading derivatives, this is somewhat slower than the direct method. For problems that are dense in the leading derivatives, this can be significantly faster.

Other issues are also discussed in the dsolve[numeric,DAE] page.

$\mathrm{dsys}≔\left\{{x\left(t\right)}^2+{y\left(t\right)}^2=1\,\mathrm{diff}\left(x\left(t\right)\,t\,t\right)=-2\mathrm{\lambda }\left(t\right)x\left(t\right)\,\mathrm{diff}\left(y\left(t\right)\,t\,t\right)=-2\mathrm{\lambda }\left(t\right)y\left(t\right)-{\mathrm{\pi }}^2\,\mathrm{\lambda }\left(0\right)=\frac{1}{200}+\frac{{\mathrm{\pi }}^2}{2}\,x\left(0\right)=0\,y\left(0\right)=-1\,\mathrm{D}\left(x\right)\left(0\right)=\frac{1}{10}\,\mathrm{D}\left(y\right)\left(0\right)=0\right\}$

$\mathrm{dsys}≔\left\{{x\left(t\right)}^2+{y\left(t\right)}^2=1\,\frac{{\ⅆ}^2}{\ⅆt^2}x\left(t\right)=-2\mathrm{\lambda }\left(t\right)x\left(t\right)\,\frac{{\ⅆ}^2}{\ⅆt^2}y\left(t\right)=-2\mathrm{\lambda }\left(t\right)y\left(t\right)-{\mathrm{\pi }}^2\,\mathrm{\lambda }\left(0\right)=\frac{1}{200}+\frac{{\mathrm{\pi }}^2}{2}\,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}\,\mathrm{numeric}\,\mathrm{method}=\mathrm{rkf45\_dae}\,\mathrm{abserr}=1.\times {10}^{\mathrm{-7}}\,\mathrm{relerr}=1.\times {10}^{\mathrm{-7}}\,\mathrm{maxfun}=0\right)$

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

$\mathrm{t1}≔\mathrm{time}\left(\right)\:$

$\mathrm{r1}≔\mathrm{dsol1}\left(2000\right)$

$\mathrm{r1}≔\left[t=2000.\,\mathrm{\lambda }\left(t\right)=4.93868361756221\,x\left(t\right)=\mathrm{-0.0123807482061333}\,\frac{\ⅆ}{\ⅆt}x\left(t\right)=0.0922928793706959\,y\left(t\right)=\mathrm{-0.999923355301742}\,\frac{\ⅆ}{\ⅆt}y\left(t\right)=\mathrm{-0.00114274175827234}\right]$

$\mathrm{t1}≔\mathrm{time}\left(\right)-\mathrm{t1}$

$\mathrm{t1}≔0.783$

$\mathrm{dsol2}≔\mathrm{dsolve}\left(\mathrm{dsys}\,\mathrm{numeric}\,\mathrm{method}=\mathrm{rkf45\_dae}\,\mathrm{projection}=\mathrm{false}\,\mathrm{abserr}=1.\times {10}^{\mathrm{-7}}\,\mathrm{relerr}=1.\times {10}^{\mathrm{-7}}\,\mathrm{maxfun}=0\right)$

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

$\mathrm{t2}≔\mathrm{time}\left(\right)\:$

$\mathrm{r2}≔\mathrm{dsol2}\left(2000\right)$

$\mathrm{r2}≔\left[t=2000.\,\mathrm{\lambda }\left(t\right)=4.94243519773441\,x\left(t\right)=0.0129846594419248\,\frac{\ⅆ}{\ⅆt}x\left(t\right)=0.0914835724074168\,y\left(t\right)=\mathrm{-0.999134421601209}\,\frac{\ⅆ}{\ⅆt}y\left(t\right)=0.00118969565389328\right]$

$\mathrm{t2}≔\mathrm{time}\left(\right)-\mathrm{t2}\:$

$\mathrm{map}\left(\mathrm{rhs}\,\mathrm{r1}\right)-\mathrm{map}\left(\mathrm{rhs}\,\mathrm{r2}\right)$

$\left[0.\,\mathrm{-0.00375158017219857}\,\mathrm{-0.0253654076480581}\,0.000809306963279091\,\mathrm{-0.000788933700533234}\,\mathrm{-0.00233243741216563}\right]$

$\mathrm{t1},\mathrm{t2}$

$0.783,0.590$

$\mathrm{dsol3}≔\mathrm{dsolve}\left(\mathrm{dsys}\,\mathrm{numeric}\,\mathrm{method}=\mathrm{rkf45\_dae}\,\mathrm{differential}=\mathrm{true}\,\mathrm{abserr}=1.\times {10}^{\mathrm{-7}}\,\mathrm{relerr}=1.\times {10}^{\mathrm{-7}}\,\mathrm{maxfun}=0\right)$

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

$\mathrm{t3}≔\mathrm{time}\left(\right)\:$

$\mathrm{r3}≔\mathrm{dsol3}\left(2000\right)$

$\mathrm{r3}≔\left[t=2000.\,\mathrm{\lambda }\left(t\right)=4.93868637305339\,x\left(t\right)=\mathrm{-0.0123954832429250}\,\frac{\ⅆ}{\ⅆt}x\left(t\right)=0.0923324756218421\,y\left(t\right)=\mathrm{-0.999923173037470}\,\frac{\ⅆ}{\ⅆt}y\left(t\right)=\mathrm{-0.00114459377988479}\right]$

$\mathrm{t3}≔\mathrm{time}\left(\right)-\mathrm{t3}\:$

$\mathrm{map}\left(\mathrm{rhs}\,\mathrm{r1}\right)-\mathrm{map}\left(\mathrm{rhs}\,\mathrm{r3}\right)$

$\left[0.\,\mathrm{-2.75549117834117}\times {10}^{\mathrm{-6}}\,0.0000147350367917121\,\mathrm{-0.0000395962511461900}\,\mathrm{-1.82264271897381}\times {10}^{\mathrm{-7}}\,1.85202161244891\times {10}^{\mathrm{-6}}\right]$

$\mathrm{t1},\mathrm{t3}$

$0.783,0.731$

$\mathrm{rates}≔\left\{K\left[s\right]=\frac{\mathrm{FLB\_ZHU}\left(t\right)}{\mathrm{FLB}\left(t\right)\mathrm{ZHU}\left(t\right)}\,r\left[1\right]\left(t\right)=k\left[1\right]{\mathrm{FLB}\left(t\right)}^4\mathrm{sqrt}\left(\mathrm{CO2}\left(t\right)\right)\,r\left[2\right]\left(t\right)=k\left[2\right]\mathrm{FLBT}\left(t\right)\mathrm{ZHU}\left(t\right)\,r\left[3\right]\left(t\right)=\frac{k\left[2\right]}{K}\mathrm{FLB}\left(t\right)\mathrm{ZLA}\left(t\right)\,r\left[4\right]\left(t\right)=k\left[3\right]\mathrm{FLB}\left(t\right){\mathrm{ZHU}\left(t\right)}^2\,r\left[5\right]\left(t\right)=k\left[4\right]{\mathrm{FLB\_ZHU}\left(t\right)}^2\mathrm{sqrt}\left(\mathrm{CO2}\left(t\right)\right)\right\}\:$

$\mathrm{inflow}≔F\left[\mathrm{In}\right]\left(t\right)=\mathrm{klA}\left(\frac{\mathrm{pCO2}}{H}-\mathrm{CO2}\left(t\right)\right)\:$

$\mathrm{conc}≔\left\{\mathrm{diff}\left(\mathrm{CO2}\left(t\right)\,t\right)=-0.5r\left[1\right]\left(t\right)-r\left[4\right]\left(t\right)-0.5r\left[5\right]\left(t\right)+F\left[\mathrm{In}\right]\left(t\right)\,\mathrm{diff}\left(\mathrm{FLB}\left(t\right)\,t\right)=-2r\left[1\right]\left(t\right)+r\left[2\right]\left(t\right)-r\left[3\right]\left(t\right)-r\left[4\right]\left(t\right)\,\mathrm{diff}\left(\mathrm{FLBT}\left(t\right)\,t\right)=r\left[1\right]\left(t\right)-r\left[2\right]\left(t\right)+r\left[3\right]\left(t\right)\,\mathrm{diff}\left(\mathrm{ZHU}\left(t\right)\,t\right)=-r\left[2\right]\left(t\right)+r\left[3\right]\left(t\right)-2r\left[4\right]\left(t\right)\,\mathrm{diff}\left(\mathrm{ZLA}\left(t\right)\,t\right)=r\left[2\right]\left(t\right)-r\left[3\right]\left(t\right)+r\left[5\right]\left(t\right)\right\}\:$

$\mathrm{ics}≔\left\{\mathrm{CO2}\left(0\right)=0.00123\,\mathrm{FLB}\left(0\right)=0.444\,\mathrm{FLBT}\left(0\right)=0\,\mathrm{ZHU}\left(0\right)=0.007\,\mathrm{ZLA}\left(0\right)=0\right\}\:$

$\mathrm{pars}≔\left\{H=737\,K=34.4\,\mathrm{klA}=3.3\,\mathrm{pCO2}=0.9\,K\left[s\right]=115.83\,k\left[1\right]=18.7\,k\left[2\right]=0.58\,k\left[3\right]=0.09\,k\left[4\right]=0.42\right\}\:$

$\mathrm{rates}≔\mathrm{eval}\left(\mathrm{rates}\,\mathrm{pars}\right)\:$

$\mathrm{inflow}≔\mathrm{eval}\left(\mathrm{inflow}\,\mathrm{pars}\right)\:$

$\mathrm{dsol4}≔\mathrm{dsolve}\left(\left\{\mathrm{inflow}\,\mathrm{conc}\left[\right]\,\mathrm{ics}\left[\right]\,\mathrm{rates}\left[\right]\right\}\,\mathrm{numeric}\,\mathrm{method}=\mathrm{rosenbrock\_dae}\,\mathrm{abserr}=1.\times {10}^{\mathrm{-10}}\,\mathrm{relerr}=1.\times {10}^{\mathrm{-10}}\right)$

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

$\mathrm{plots}\left[\mathrm{odeplot}\right]\left(\mathrm{dsol4}\,\left[t\,\mathrm{CO2}\left(t\right)\right]\,0..0.5\right)$

$\mathrm{eval}\left(\left[\mathrm{CO2}\left(t\right)\,\mathrm{ZLA}\left(t\right)\right]\,\mathrm{dsol4}\left(180\right)\right)$

$\left[0.00120383147146856\,0.0170801088535847\right]$

$\mathrm{plots}\left[\mathrm{odeplot}\right]\left(\mathrm{dsol4}\,\left[\left[t\,10\mathrm{CO2}\left(t\right)\right]\,\left[t\,\mathrm{ZLA}\left(t\right)\right]\right]\,0..180\right)$