 DAE extension - Maple Help

dsolve/numeric/DAE_extension

find numerical solution of ordinary differential-algebraic initial value problems Calling Sequence dsolve(daesys, numeric, vars, method=meth, options) Parameters

 daesys - set or list; ordinary differential equation(s), algebraic equation(s) and initial conditions numeric - literal; instruct dsolve to find a numerical solution vars - (optional) dependent variable or a set or list of dependent variables for daesys method - literal; keyword for specification of the method meth - name; specification of the method, rkf45_dae, ck45_dae, or rosenbrock_dae options - (optional) equations of the form keyword = value Description

 • 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. Examples

The simple pendulum problem in natural coordinates:

 > $\mathrm{dsys}≔\left\{\frac{{ⅆ}^{2}}{ⅆ{t}^{2}}x\left(t\right)=-2\mathrm{λ}\left(t\right)x\left(t\right),\frac{{ⅆ}^{2}}{ⅆ{t}^{2}}y\left(t\right)=-2\mathrm{λ}\left(t\right)y\left(t\right)-{\mathrm{Pi}}^{2},{x\left(t\right)}^{2}+{y\left(t\right)}^{2}=1,x\left(0\right)=0,\mathrm{D}\left(x\right)\left(0\right)=\frac{1}{10},y\left(0\right)=-1,\mathrm{D}\left(y\right)\left(0\right)=0,\mathrm{λ}\left(0\right)=\frac{1}{200}+\frac{{\mathrm{Pi}}^{2}}{2}\right\}$
 ${\mathrm{dsys}}{≔}\left\{{{x}{}\left({t}\right)}^{{2}}{+}{{y}{}\left({t}\right)}^{{2}}{=}{1}{,}\frac{{{ⅆ}}^{{2}}}{{ⅆ}{{t}}^{{2}}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{x}{}\left({t}\right){=}{-}{2}{}{\mathrm{\lambda }}{}\left({t}\right){}{x}{}\left({t}\right){,}\frac{{{ⅆ}}^{{2}}}{{ⅆ}{{t}}^{{2}}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{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){=}{-1}{,}{\mathrm{D}}{}\left({x}\right){}\left({0}\right){=}\frac{{1}}{{10}}{,}{\mathrm{D}}{}\left({y}\right){}\left({0}\right){=}{0}\right\}$ (1)
 > $\mathrm{dsol1}≔\mathrm{dsolve}\left(\mathrm{dsys},\mathrm{numeric},\mathrm{method}=\mathrm{rkf45_dae},\mathrm{abserr}=1.{10}^{-7},\mathrm{relerr}=1.{10}^{-7},\mathrm{maxfun}=0\right)$
 ${\mathrm{dsol1}}{≔}{\mathbf{proc}}\left({\mathrm{x_rkf45_dae}}\right)\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{...}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathbf{end proc}}$ (2)
 > $\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){=}{-0.0123807482061333}{,}\frac{{ⅆ}}{{ⅆ}{t}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{x}{}\left({t}\right){=}{0.0922928793706959}{,}{y}{}\left({t}\right){=}{-0.999923355301742}{,}\frac{{ⅆ}}{{ⅆ}{t}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{y}{}\left({t}\right){=}{-0.00114274175827234}\right]$ (3)
 > $\mathrm{t1}≔\mathrm{time}\left(\right)-\mathrm{t1}$
 ${\mathrm{t1}}{≔}{1.007}$ (4)

Compare with projection=false:

 > $\mathrm{dsol2}≔\mathrm{dsolve}\left(\mathrm{dsys},\mathrm{numeric},\mathrm{method}=\mathrm{rkf45_dae},\mathrm{projection}=\mathrm{false},\mathrm{abserr}=1.{10}^{-7},\mathrm{relerr}=1.{10}^{-7},\mathrm{maxfun}=0\right)$
 ${\mathrm{dsol2}}{≔}{\mathbf{proc}}\left({\mathrm{x_rkf45_dae}}\right)\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{...}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathbf{end proc}}$ (5)
 > $\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}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{x}{}\left({t}\right){=}{0.0914835724074168}{,}{y}{}\left({t}\right){=}{-0.999134421601209}{,}\frac{{ⅆ}}{{ⅆ}{t}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{y}{}\left({t}\right){=}{0.00118969565389328}\right]$ (6)
 > $\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.}{,}{-0.00375158017219857}{,}{-0.0253654076480581}{,}{0.000809306963279091}{,}{-0.000788933700533234}{,}{-0.00233243741216563}\right]$ (7)
 > $\mathrm{t1},\mathrm{t2}$
 ${1.007}{,}{0.814}$ (8)

Compare time with differential=true:

 > $\mathrm{dsol3}≔\mathrm{dsolve}\left(\mathrm{dsys},\mathrm{numeric},\mathrm{method}=\mathrm{rkf45_dae},\mathrm{differential}=\mathrm{true},\mathrm{abserr}=1.{10}^{-7},\mathrm{relerr}=1.{10}^{-7},\mathrm{maxfun}=0\right)$
 ${\mathrm{dsol3}}{≔}{\mathbf{proc}}\left({\mathrm{x_rkf45_dae}}\right)\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{...}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathbf{end proc}}$ (9)
 > $\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){=}{-0.0123954832429250}{,}\frac{{ⅆ}}{{ⅆ}{t}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{x}{}\left({t}\right){=}{0.0923324756218421}{,}{y}{}\left({t}\right){=}{-0.999923173037470}{,}\frac{{ⅆ}}{{ⅆ}{t}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{y}{}\left({t}\right){=}{-0.00114459377988479}\right]$ (10)
 > $\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.}{,}{-2.75549117834117}{×}{{10}}^{{-6}}{,}{0.0000147350367917121}{,}{-0.0000395962511461900}{,}{-1.82264271897381}{×}{{10}}^{{-7}}{,}{1.85202161244891}{×}{{10}}^{{-6}}\right]$ (11)
 > $\mathrm{t1},\mathrm{t3}$
 ${1.007}{,}{1.076}$ (12)

The Chemical Akzo Nobel problem

This DAE system describes a chemical process in which FLB and ZHU are mixed while carbon dioxide is added over time. The ZLA concentration is of interest.

The rate equations for each reaction are given by:

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

The carbon dioxide inflow is governed by

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

The concentrations are described by the following differential equations

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

and the Initial values are given by:

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

The parameters are as follows:

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

To begin, we apply the parameters to the rates and inflow:

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

then pass all equations and initial conditions to dsolve with moderately tight tolerances choosing the rosenbrock_dae method, as the problem is known to be stiff:

 > $\mathrm{dsol4}≔\mathrm{dsolve}\left(\left\{{\mathrm{conc}}_{[]},{\mathrm{rates}}_{[]},\mathrm{inflow},{\mathrm{ics}}_{[]}\right\},\mathrm{numeric},\mathrm{method}=\mathrm{rosenbrock_dae},\mathrm{abserr}=1.{10}^{-10},\mathrm{relerr}=1.{10}^{-10}\right)$
 ${\mathrm{dsol4}}{≔}{\mathbf{proc}}\left({\mathrm{x_rosenbrock_dae}}\right)\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{...}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathbf{end proc}}$ (13)

The period of interest is t=0..180, but the CO2 concentration shows some interesting behavior in the first 0.5s:

 > ${\mathrm{plots}}_{\mathrm{odeplot}}\left(\mathrm{dsol4},\left[t,\mathrm{CO2}\left(t\right)\right],0..0.5\right)$ The values at t=180 show perfect comparison to the 6-digit accurate benchmark results (CO2(180)=1.20383e-3, ZLA(180)=1.70801e-2):

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

And we plot the CO2 and ZLA concentrations as a function of time over the entire range of interest (adjusting CO2 by a factor of 10 for scale):

 > ${\mathrm{plots}}_{\mathrm{odeplot}}\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)$ 