DAE - Maple Help

dsolve/numeric/DAE

find numerical solution of differential-algebraic initial value problems

 Calling Sequence dsolve(daesys, numeric, vars, options)

Parameters

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

Description

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

 'output' = keyword or array 'range' = numeric..numeric 'abserr' = numeric or list 'relerr' = numeric 'stiff' = boolean 'known' = name or list of names 'optimize' = boolean

 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' = numeric 'maxstep' = numeric 'initstep' = numeric 'startinit' = boolean 'events' = list 'event_pre' = keyword 'event_maxiter' = integer 'projection' = boolean 'differential' = boolean 'implicit' = boolean 'parameters' = list of names

 '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.

Examples

As a first example, we consider the problem of modeling the dynamics of a mass on a string of unit length in 2-D Cartesian coordinates (the pendulum problem). We let $r$ be the position of the mass on the string, and $v$ the velocity:

 > $\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]$ (1)
 > $v≔\mathrm{VectorCalculus}\left[\mathrm{diff}\right]\left(r,t\right)$
 ${v}{≔}\left[\begin{array}{c}\frac{{ⅆ}}{{ⅆ}{t}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{x}{}\left({t}\right)\\ \frac{{ⅆ}}{{ⅆ}{t}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{y}{}\left({t}\right)\end{array}\right]$ (2)

Since the mass is on a string of fixed length $l$, we have the constraint:

 > $\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}}$ (3)

Now we want to construct the DAE system using the Euler-Lagrange formulation, so we compute the kinetic energy $T$ and the potential energy $V$ as:

 > $T≔\frac{m}{2}\mathrm{VectorCalculus}\left[\mathrm{DotProduct}\right]\left(v,v\right)$
 ${T}{≔}\frac{{m}{}\left({\left(\frac{{ⅆ}}{{ⅆ}{t}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{x}{}\left({t}\right)\right)}^{{2}}{+}{\left(\frac{{ⅆ}}{{ⅆ}{t}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{y}{}\left({t}\right)\right)}^{{2}}\right)}{{2}}$ (4)
 > $V≔mgy\left(t\right)$
 ${V}{≔}{m}{}{g}{}{y}{}\left({t}\right)$ (5)

where $m$ is the mass, and $g$ is the gravitational constant (we will use 9.8). The Lagrangian and modified Lagrangian are given by:

 > $L≔T-V$
 ${L}{≔}\frac{{m}{}\left({\left(\frac{{ⅆ}}{{ⅆ}{t}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{x}{}\left({t}\right)\right)}^{{2}}{+}{\left(\frac{{ⅆ}}{{ⅆ}{t}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{y}{}\left({t}\right)\right)}^{{2}}\right)}{{2}}{-}{m}{}{g}{}{y}{}\left({t}\right)$ (6)
 > $\mathrm{ML}≔L-\mathrm{\lambda }\left(t\right)\mathrm{con}$
 ${\mathrm{ML}}{≔}\frac{{m}{}\left({\left(\frac{{ⅆ}}{{ⅆ}{t}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{x}{}\left({t}\right)\right)}^{{2}}{+}{\left(\frac{{ⅆ}}{{ⅆ}{t}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{y}{}\left({t}\right)\right)}^{{2}}\right)}{{2}}{-}{m}{}{g}{}{y}{}\left({t}\right){-}{\mathrm{\lambda }}{}\left({t}\right){}\left({{x}{}\left({t}\right)}^{{2}}{+}{{y}{}\left({t}\right)}^{{2}}{-}{{l}}^{{2}}\right)$ (7)

We can then construct the Euler-Lagrange formulation via:

 > $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\{{-}{2}{}{\mathrm{\lambda }}{}\left({t}\right){}{x}{}\left({t}\right){-}{m}{}\left(\frac{{{ⅆ}}^{{2}}}{{ⅆ}{{t}}^{{2}}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{x}{}\left({t}\right)\right){,}{{l}}^{{2}}{-}{{x}{}\left({t}\right)}^{{2}}{-}{{y}{}\left({t}\right)}^{{2}}{,}{-}{m}{}{g}{-}{2}{}{y}{}\left({t}\right){}{\mathrm{\lambda }}{}\left({t}\right){-}{m}{}\left(\frac{{{ⅆ}}^{{2}}}{{ⅆ}{{t}}^{{2}}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{y}{}\left({t}\right)\right)\right\}$ (8)
 > $\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}}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{x}{}\left({t}\right){=}{-}\frac{{2}{}{\mathrm{\lambda }}{}\left({t}\right){}{x}{}\left({t}\right)}{{m}}$ (9)
 > $\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}}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{y}{}\left({t}\right){=}{-}\frac{{m}{}{g}{+}{2}{}{y}{}\left({t}\right){}{\mathrm{\lambda }}{}\left({t}\right)}{{m}}$ (10)
 > $\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}}$ (11)

Now we have the equations of motion for the pendulum. Next, we need to determine consistent initial conditions. To do so, we must identify any hidden constraints of the system. These are easy to find, as we have only one constraint.

 > $\mathrm{Dcon}≔\mathrm{diff}\left(\mathrm{EQc},t\right)$
 ${\mathrm{Dcon}}{≔}{-}{2}{}{x}{}\left({t}\right){}\left(\frac{{ⅆ}}{{ⅆ}{t}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{x}{}\left({t}\right)\right){-}{2}{}{y}{}\left({t}\right){}\left(\frac{{ⅆ}}{{ⅆ}{t}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{y}{}\left({t}\right)\right)$ (12)
 > $\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}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{x}{}\left({t}\right)\right)}^{{2}}{+}\frac{{4}{}{{x}{}\left({t}\right)}^{{2}}{}{\mathrm{\lambda }}{}\left({t}\right)}{{m}}{-}{2}{}{\left(\frac{{ⅆ}}{{ⅆ}{t}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{y}{}\left({t}\right)\right)}^{{2}}{+}\frac{{2}{}{y}{}\left({t}\right){}\left({m}{}{g}{+}{2}{}{y}{}\left({t}\right){}{\mathrm{\lambda }}{}\left({t}\right)\right)}{{m}}$ (13)

Our initial conditions must satisfy $\mathrm{EQc}$, $\mathrm{Dcon}$, and $\mathrm{DDcon}$ at the initial point, leaving only 2 degrees of freedom for the conditions. So for a pendulum starting at the minimum value of $y\left(0\right)=-l$ having an initial horizontal velocity of $\mathrm{D}\left(x\right)\left(0\right)=\mathrm{vx}$, we get:

 > $\mathrm{sys}≔\left\{y\left(0\right)=-l,\mathrm{D}\left(x\right)\left(0\right)=\mathrm{vx}\right\}\phantom{\rule[-0.0ex]{0.3em}{0.0ex}}\mathbf{union}\phantom{\rule[-0.0ex]{0.3em}{0.0ex}}\mathrm{eval}\left(\mathrm{convert}\left(\left\{\mathrm{DDcon},\mathrm{Dcon},\mathrm{EQc}\right\},\mathrm{D}\right),t=0\right)$
 ${\mathrm{sys}}{≔}\left\{{-}{2}{}{x}{}\left({0}\right){}{\mathrm{D}}{}\left({x}\right){}\left({0}\right){-}{2}{}{y}{}\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{{2}{}{y}{}\left({0}\right){}\left({m}{}{g}{+}{2}{}{y}{}\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\}$ (14)
 > $\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({g}{}{l}{+}{{\mathrm{vx}}}^{{2}}\right)}{{2}{}{{l}}^{{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\}$ (15)

So we consider the above with a pendulum of unit length $l=1$ having unit mass $m=1$ and an initial horizontal velocity of $\mathrm{vx}=\frac{1}{10}$, giving us the DAE system and initial conditions:

 > $\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}}}\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){=}{-}{9.8}{-}{2}{}{y}{}\left({t}\right){}{\mathrm{\lambda }}{}\left({t}\right)\right\}$ (16)
 > $\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){=}{-1}{,}{\mathrm{D}}{}\left({x}\right){}\left({0}\right){=}\frac{{1}}{{10}}{,}{\mathrm{D}}{}\left({y}\right){}\left({0}\right){=}{0}\right\}$ (17)

We can then obtain the solution as:

 > $\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{dini},\mathrm{numeric}\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}}$ (18)
 > $\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}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{x}{}\left({t}\right){=}{0.000564539406239432}{,}{y}{}\left({t}\right){=}{-0.999489811490763}{,}\frac{{ⅆ}}{{ⅆ}{t}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{y}{}\left({t}\right){=}{0.0000180384378593043}\right]$ (19)
 > $\mathrm{dsol1}\left(1\right)$
 $\left[{t}{=}{1.}{,}{\mathrm{\lambda }}{}\left({t}\right){=}{4.90499905723651}{,}{x}{}\left({t}\right){=}{0.000360891517794992}{,}\frac{{ⅆ}}{{ⅆ}{t}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{x}{}\left({t}\right){=}{-0.0999937504894365}{,}{y}{}\left({t}\right){=}{-0.999999934755133}{,}\frac{{ⅆ}}{{ⅆ}{t}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{y}{}\left({t}\right){=}{-0.0000360794101890703}\right]$ (20)

Solution with rosenbrock_dae:

 > $\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{dini},\mathrm{numeric},\mathrm{method}=\mathrm{rosenbrock_dae}\right)$
 ${\mathrm{dsol2}}{≔}{\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}}$ (21)
 > $\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}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{x}{}\left({t}\right){=}{0.000564626004686083}{,}{y}{}\left({t}\right){=}{-0.999489813543334}{,}\frac{{ⅆ}}{{ⅆ}{t}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{y}{}\left({t}\right){=}{0.0000180482728717227}\right]$ (22)
 > $\mathrm{dsol2}\left(1\right)$
 $\left[{t}{=}{1.}{,}{\mathrm{\lambda }}{}\left({t}\right){=}{4.90499906207811}{,}{x}{}\left({t}\right){=}{0.000360968463263150}{,}\frac{{ⅆ}}{{ⅆ}{t}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{x}{}\left({t}\right){=}{-0.0999936779452163}{,}{y}{}\left({t}\right){=}{-0.999999936691318}{,}\frac{{ⅆ}}{{ⅆ}{t}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{y}{}\left({t}\right){=}{-0.0000360897908606990}\right]$ (23)

Solution with mebdfi:

 > $\mathrm{dsol3}≔\mathrm{dsolve}\left(\mathrm{dsys}\phantom{\rule[-0.0ex]{0.3em}{0.0ex}}\mathbf{union}\phantom{\rule[-0.0ex]{0.3em}{0.0ex}}\mathrm{dini},\mathrm{numeric},\mathrm{method}=\mathrm{mebdfi}\right)$
 ${\mathrm{dsol3}}{≔}{\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}}$ (24)
 > $\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}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{x}{}\left({t}\right){=}{0.000565199414553933343}{,}{y}{}\left({t}\right){=}{-0.999489820801930273}{,}\frac{{ⅆ}}{{ⅆ}{t}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{y}{}\left({t}\right){=}{0.0000180581765179766362}\right]$ (25)
 > $\mathrm{dsol3}\left(1\right)$
 $\left[{t}{=}{1.}{,}{\mathrm{\lambda }}{}\left({t}\right){=}{4.90499581413915475}{,}{x}{}\left({t}\right){=}{0.000361181747066772996}{,}\frac{{ⅆ}}{{ⅆ}{t}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{x}{}\left({t}\right){=}{-0.0999926673290293111}{,}{y}{}\left({t}\right){=}{-0.99999993019115307}{,}\frac{{ⅆ}}{{ⅆ}{t}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{y}{}\left({t}\right){=}{-0.0000360957473722881540}\right]$ (26)

Now consider a similar problem as above, but in addition add a second mass supported from the first by another string, this one of length 1/2 (the double pendulum). The system can be obtained and solved as:

 > $\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}}}\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){,}{\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}}}\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\}$ (27)
 > $\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\}$ (28)
 > $\mathrm{dsold}≔\mathrm{dsolve}\left(\mathrm{dsysd}\phantom{\rule[-0.0ex]{0.3em}{0.0ex}}\mathbf{union}\phantom{\rule[-0.0ex]{0.3em}{0.0ex}}\mathrm{ics},\mathrm{numeric}\right)$
 ${\mathrm{dsold}}{≔}{\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}}$ (29)

The trajectory of the second mass can be plotted via:

 > $\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)$

Note that we did not specify initial conditions for $\mathrm{λ1}\left(t\right)$ and $\mathrm{λ2}\left(t\right)$ as we were using an extension method, and the system was sufficiently linear in the lambda.