mebdfi - Maple Help

dsolve/numeric/mebdfi

numerical solution of differential-algebraic equations

 Calling Sequence dsolve($\mathrm{daesys}$, numeric, method=mebdfi, vars, options)

Parameters

 daesys - set or list; differential-algebraic system of equation and initial conditions numeric - literal name; instruct dsolve to find a numerical solution method=mebdfi - literal equation; numerical method to use vars - (optional) a set or list of dependent variables for daesys options - (optional) equations of the form keyword = value

Description

 • The dsolve command with the options numeric and method=mebdfi finds a numerical solution of a DAE system using the Modified Extended Backward Differentiation Equation Implicit method. This is a stiff method, so can handle stiff problems efficiently, but should only be used for DAE of index 2 or lower. It can be used to solve regular ODE problems, but use of regular ODE solvers is recommended for that purpose (see dsolve[numeric,IVP]).
 • The following options are available for the mebdfi method:

 'output' = keyword or array 'known' = name or list of names 'startinit' = boolean 'optimize' = boolean 'maxfun' = integer 'abserr' = numeric 'relerr' = numeric 'minstep' = numeric 'maxstep' = numeric 'initstep' = numeric 'maxord' = integer

 'output' and 'known'
 The 'output' option specifies the output from dsolve, and the known option specifies user-defined known functions. These options are discussed in dsolve[numeric].
 'startinit' and 'optimize'
 These options control the method and behavior of the computation, the 'startinit' option is the same as for IVP problems, and is discussed in dsolve[numeric,IVP], while the 'optimize' option is discussed in dsolve[numeric].
 'maxfun'
 This option specifies the maximum number of steps taken to obtain the requested solution. For this direct solver, this is the closest parallel to dsolve[maxfun]. By default this option is disabled.
 'abserr' and 'relerr'
 These options specify the desired accuracy of the solution, and are discussed in dsolve[Error_Control]. The default values for mebdfi are $\mathrm{abserr}=1.×{10}^{-7}$ and $\mathrm{relerr}=1.×{10}^{-6}$.
 'minstep', 'maxstep', and 'initstep'
 These options provide finer control over the step size used in the method, and are also discussed in dsolve[Error_Control]. By default $\mathrm{minstep}$ is determined within the computation, $\mathrm{maxstep}$ is disabled, and $\mathrm{initstep}=\mathrm{relerr}$. When setting $\mathrm{minstep}$ ensure that the value is adequately small to allow for the first few steps, which are computed at a lower order.
 'maxord'
 This option specifies the maximum order of the method used in the course of the computation as an integer between 2 and 8. The default value is $8$, but for difficult problems it may be necessary to specify this as $4$ or $5$. Lower order gives better stability, but is less efficient.
 • The computation may return with an error message corresponding to an error condition of the mebdfi procedure.
 • Results can be plotted using the function odeplot in the plots package.

Examples

Double pendulum (index-reduced to index-2):

 > $\mathrm{dsys}≔\left\{\left(\mathrm{x1}\left(t\right)-\mathrm{x2}\left(t\right)\right)\mathrm{diff}\left(\mathrm{x1}\left(t\right)-\mathrm{x2}\left(t\right),t\right)+\left(\mathrm{y1}\left(t\right)-\mathrm{y2}\left(t\right)\right)\mathrm{diff}\left(\mathrm{y1}\left(t\right)-\mathrm{y2}\left(t\right),t\right),\mathrm{x1}\left(t\right)\mathrm{diff}\left(\mathrm{x1}\left(t\right),t\right)+\mathrm{y1}\left(t\right)\mathrm{diff}\left(\mathrm{y1}\left(t\right),t\right),\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),\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{dsys}}{≔}\left\{\left({\mathrm{x1}}{}\left({t}\right){-}{\mathrm{x2}}{}\left({t}\right)\right){}\left(\frac{{ⅆ}}{{ⅆ}{t}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{\mathrm{x1}}{}\left({t}\right){-}\frac{{ⅆ}}{{ⅆ}{t}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{\mathrm{x2}}{}\left({t}\right)\right){+}\left({\mathrm{y1}}{}\left({t}\right){-}{\mathrm{y2}}{}\left({t}\right)\right){}\left(\frac{{ⅆ}}{{ⅆ}{t}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{\mathrm{y1}}{}\left({t}\right){-}\frac{{ⅆ}}{{ⅆ}{t}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{\mathrm{y2}}{}\left({t}\right)\right){,}{\mathrm{x1}}{}\left({t}\right){}\left(\frac{{ⅆ}}{{ⅆ}{t}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{\mathrm{x1}}{}\left({t}\right)\right){+}{\mathrm{y1}}{}\left({t}\right){}\left(\frac{{ⅆ}}{{ⅆ}{t}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{\mathrm{y1}}{}\left({t}\right)\right){,}\frac{{{ⅆ}}^{{2}}}{{ⅆ}{{t}}^{{2}}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{\mathrm{x2}}{}\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}}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{\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}}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{\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}}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{\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\}$ (1)
 > $\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){=}{-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\}$ (2)

Default solution:

 > $\mathrm{dsol1}≔\mathrm{dsolve}\left(\mathrm{dsys}\phantom{\rule[-0.0ex]{0.3em}{0.0ex}}\mathbf{union}\phantom{\rule[-0.0ex]{0.3em}{0.0ex}}\mathrm{ics},\mathrm{numeric},\mathrm{method}=\mathrm{mebdfi}\right)$
 ${\mathrm{dsol1}}{≔}{\mathbf{proc}}\left({\mathrm{x_mebdfi}}\right)\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{...}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathbf{end proc}}$ (3)
 > $\mathrm{t1}≔\mathrm{time}\left(\right):$
 > $\mathrm{dsol1}\left(10\right)$
 $\left[{t}{=}{10.}{,}{\mathrm{λ1}}{}\left({t}\right){=}{8.05374866836774572}{,}{\mathrm{λ2}}{}\left({t}\right){=}{-6.40533775902082869}{,}{\mathrm{x1}}{}\left({t}\right){=}{-0.355531961724009915}{,}\frac{{ⅆ}}{{ⅆ}{t}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{\mathrm{x1}}{}\left({t}\right){=}{-1.25978921249117670}{,}{\mathrm{x2}}{}\left({t}\right){=}{-0.492864788890956418}{,}\frac{{ⅆ}}{{ⅆ}{t}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{\mathrm{x2}}{}\left({t}\right){=}{-1.05303335398026632}{,}{\mathrm{y1}}{}\left({t}\right){=}{-0.934664111102251338}{,}\frac{{ⅆ}}{{ⅆ}{t}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{\mathrm{y1}}{}\left({t}\right){=}{0.479204587580475905}{,}{\mathrm{y2}}{}\left({t}\right){=}{-0.453894585514793392}{,}\frac{{ⅆ}}{{ⅆ}{t}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{\mathrm{y2}}{}\left({t}\right){=}{0.538264833117563346}\right]$ (4)
 > $\mathrm{time}\left(\right)-\mathrm{t1}$
 ${0.917}$ (5)

Lower order for stability:

 > $\mathrm{dsol2}≔\mathrm{dsolve}\left(\mathrm{dsys}\phantom{\rule[-0.0ex]{0.3em}{0.0ex}}\mathbf{union}\phantom{\rule[-0.0ex]{0.3em}{0.0ex}}\mathrm{ics},\mathrm{numeric},\mathrm{method}=\mathrm{mebdfi},\mathrm{maxord}=4\right)$
 ${\mathrm{dsol2}}{≔}{\mathbf{proc}}\left({\mathrm{x_mebdfi}}\right)\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{...}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}{\mathbf{end proc}}$ (6)
 > $\mathrm{t2}≔\mathrm{time}\left(\right):$
 > $\mathrm{dsol2}\left(10\right)$
 $\left[{t}{=}{10.}{,}{\mathrm{λ1}}{}\left({t}\right){=}{8.05124380921260041}{,}{\mathrm{λ2}}{}\left({t}\right){=}{-6.42065568374894546}{,}{\mathrm{x1}}{}\left({t}\right){=}{-0.357252010506513173}{,}\frac{{ⅆ}}{{ⅆ}{t}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{\mathrm{x1}}{}\left({t}\right){=}{-1.25072569463780869}{,}{\mathrm{x2}}{}\left({t}\right){=}{-0.492112228851291800}{,}\frac{{ⅆ}}{{ⅆ}{t}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{\mathrm{x2}}{}\left({t}\right){=}{-1.04698380844525074}{,}{\mathrm{y1}}{}\left({t}\right){=}{-0.934008348789023168}{,}\frac{{ⅆ}}{{ⅆ}{t}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{\mathrm{y1}}{}\left({t}\right){=}{0.478394293227407508}{,}{\mathrm{y2}}{}\left({t}\right){=}{-0.452537634291989843}{,}\frac{{ⅆ}}{{ⅆ}{t}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{\mathrm{y2}}{}\left({t}\right){=}{0.535462511444270373}\right]$ (7)
 > $\mathrm{time}\left(\right)-\mathrm{t2}$
 ${1.253}$ (8)

References

 Cash, J.R. "The Integration of stiff IVP in ODE using modified extended BDF." Computers and Mathematics with Applications. Vol. 9. (1983): 645-657.
 Cash, J.R., and Considine, S. "An MEBDF code for stiff IVP." ACM Trans Math Software. 1992: 142-158.