DAE Example Worksheet

All Products Maple MapleSim

Home : Support : Online Help : Mathematics : Differential-algebraic Equations : DAE Example Worksheet

Introduction to Numerical Differential-Algebraic Equation Solvers

Differential Equations and Differential-Algebraic Equations (DAEs)

Differential Equations - ODE

Differential equations describe the behavior of a system through rates of change.

Simplest Example:

>	$restart &colon;$

>	$\frac{&DifferentialD;}{&DifferentialD; t} y (t) = 0$

$\frac{&DifferentialD;}{&DifferentialD; t} y (t) = 0$

(1.1.1)

>	$dsolve (\{, y (0) = 1\})$

$y (t) = 1$

(1.1.2)

This can extend to systems with more than one dependent variable.

>	$dsys := \{\frac{&DifferentialD;}{&DifferentialD; t} x (t) = y (t), \frac{&DifferentialD;}{&DifferentialD; t} y (t) = - x (t), x (0) = 0, y (0) = 1\}$

$dsys ≔ \{\frac{&DifferentialD;}{&DifferentialD; t} x (t) = y (t), \frac{&DifferentialD;}{&DifferentialD; t} y (t) = - x (t), x (0) = 0, y (0) = 1\}$

(1.1.3)

>	$dsolve (dsys)$

$\{x (t) = \sin (t), y (t) = \cos (t)\}$

(1.1.4)

It is not always possible to find the closed form solution of an ODE problem. This creates a need for solvers that compute a numerical approximation to the solution.

Differential Algebraic Equations - DAE

A DAE is a system that is a combination of differential equations and constraints.

One of the simplest examples is the pendulum system, where the m*g (mass times gravitational constant) term has been replaced by Pi^2.

>	$dsys := [\frac{{&DifferentialD;}^{2}}{&DifferentialD; t^{2}} x (t) = - 2 λ (t) x (t), \frac{{&DifferentialD;}^{2}}{&DifferentialD; t^{2}} y (t) = - 2 λ (t) y (t) - π^{2}, {x (t)}^{2} + {y (t)}^{2} = 1]$

$dsys ≔ [\frac{{&DifferentialD;}^{2}}{&DifferentialD; t^{2}} x (t) = - 2 λ (t) x (t), \frac{{&DifferentialD;}^{2}}{&DifferentialD; t^{2}} y (t) = - 2 λ (t) y (t) - π^{2}, {x (t)}^{2} + {y (t)}^{2} = 1]$

(1.2.1)

To display derivatives and functions using a simpler notation, use the PDEtools[declare] command.

>	$PDEtools [declare] (\{x (t), y (t), θ (t), r (t), λ (t)\}, prime = t, quiet) &semi;$ $dsys$

$[x'' = - 2 λ x, y'' = - 2 λ y - π^{2}, x^{2} + y^{2} = 1]$

(1.2.2)

The typical solution approach for this problem (as presented in textbooks) is to apply a change of variables to polar co-ordinates, which would proceed as follows.

>	$tr := \{x (t) = r (t) \sin (θ (t)), y (t) = - r (t) \cos (θ (t))\} &semi;$ $nsys ≔ PDEtools [dchange] (tr, dsys, [r (t), θ (t)])$

$tr ≔ \{x = r \sin (θ), y = - r \cos (θ)\}$

$nsys ≔ [(r'') \sin (θ) + 2 (r') (θ') \cos (θ) + r (θ'') \cos (θ) - r {(θ')}^{2} \sin (θ) = - 2 λ r \sin (θ), - (r'') \cos (θ) + 2 (r') (θ') \sin (θ) + r (θ'') \sin (θ) + r {(θ')}^{2} \cos (θ) = 2 λ r \cos (θ) - π^{2}, r^{2} {\sin (θ)}^{2} + r^{2} {\cos (θ)}^{2} = 1]$

(1.2.3)

The last equation implies r=1.

>	$simplify (nsys [- 1], trig)$

$r^{2} = 1$

(1.2.4)

Evaluate and remove lambda(t) as follows, giving an ODE in theta(t).

>	$nsys := \binom{nsys}{} \| \binom{}{r (t) = 1} &semi;$ $isolate (nsys [1], λ (t)) &semi;$

$nsys ≔ [(θ'') \cos (θ) - {(θ')}^{2} \sin (θ) = - 2 λ \sin (θ), (θ'') \sin (θ) + {(θ')}^{2} \cos (θ) = 2 λ \cos (θ) - π^{2}, {\sin (θ)}^{2} + {\cos (θ)}^{2} = 1]$

$λ = \frac{- (θ'') \cos (θ) + {(θ')}^{2} \sin (θ)}{2 \sin (θ)}$

(1.2.5)

>	$eval (nsys [2],) &semi;$

$(θ'') \sin (θ) + {(θ')}^{2} \cos (θ) = \frac{(- (θ'') \cos (θ) + {(θ')}^{2} \sin (θ)) \cos (θ)}{\sin (θ)} - π^{2}$

(1.2.6)

>	$pode ≔ simplify (isolate (, \frac{{&DifferentialD;}^{2}}{&DifferentialD; t^{2}} θ (t)), trig)$

$pode ≔ θ'' = - \sin (θ) π^{2}$

(1.2.7)

This is the more commonly known pendulum equation.

As an approximation for small displacements, we can approximate sin(theta) ~ theta giving us the simple harmonic oscillator equation:

>	$\frac{{&DifferentialD;}^{2}}{&DifferentialD; t^{2}} θ (t) + π^{2} θ (t) = 0$

$θ'' + π^{2} θ = 0$

(1.2.8)

This is appropriate for the pendulum equation, but, in general, such a coordinate transform may not exist. Hence, there is the need for numerical solution methods for DAEs.

Structure of DAE

DAEs have an interesting property.

$dsys$

$[x'' = - 2 λ x, y'' = - 2 λ y - π^{2}, x^{2} + y^{2} = 1]$

(1.3.1)

Consider our constraint:

>	$c1 ≔ dsys [- 1]$

$c1 ≔ x^{2} + y^{2} = 1$

(1.3.2)

Differentiate:

>	$c2 := \frac{\frac{\partial}{\partial t} c1}{2}$

$c2 ≔ x (x') + y (y') = 0$

(1.3.3)

Differentiate, eliminate, and isolate:

>	$c3 := eval (\frac{\partial}{\partial t} c2, dsys [1 .. 2]) &semi;$ $c3 := factor (isolate (c3, λ (t)))$

$c3 ≔ {(x')}^{2} - 2 x^{2} λ + {(y')}^{2} + y (- 2 λ y - π^{2}) = 0$

$c3 ≔ λ = - \frac{- {(x')}^{2} - {(y')}^{2} + y π^{2}}{2 (x^{2} + y^{2})}$

(1.3.4)

One more differentiation and elimination yields an ODE for lambda(t).

>	$d4 := eval (\frac{\partial}{\partial t} c3, dsys [1 .. 2])$

$d4 ≔ λ' = - \frac{4 (x') λ x - 2 (y') (- 2 λ y - π^{2}) + (y') π^{2}}{2 (x^{2} + y^{2})} + \frac{(- {(x')}^{2} - {(y')}^{2} + y π^{2}) (2 x (x') + 2 y (y'))}{2 {(x^{2} + y^{2})}^{2}}$

(1.3.5)

The above process of successive differentiation until we obtain an ODE for each variable is called index reduction.
The number of required differentiations is called the gear index. It sometimes corresponds to the difficulty in obtaining a numerical solution of the DAE problem. (A greater index corresponds to a more difficult problem.)

Issues:

We obtained 3 algebraic relations for lambda, x, x', y, y': c1,c2,c3. Initial conditions must satisfy these—as a result, we really only have 2 degrees of freedom in the initial conditions for a problem that would normally have 5.
This restricts the initial conditions. They must be checked for consistency.
Note: Constraints are typically nonlinear, and cannot normally be solved in terms of the free initial conditions.

The algebraic constraints must hold along the entire solution.

Solution Methods

Method of Index Reduction

As previously shown, this involves successive differentiation of system until the DAE becomes an ODE system.

Advantages

•	Allows us to obtain all constraints directly

•	Allows use of a standard ODE solver to obtain the solution

Disadvantages

•	Does not enforce the constraints, other than at the initial point (resulting in constraint drift)

Example:

>	$IRdsys := \binom{\{dsys [1], dsys [2], d4\}}{} \| \binom{}{g = π^{2}} &semi;$ $IRcons ≔ \binom{[c1, c2, c3]}{} \| \binom{}{g = π^{2}}$

$IRdsys ≔ \{x'' = - 2 λ x, y'' = - 2 λ y - π^{2}, λ' = - \frac{4 (x') λ x - 2 (y') (- 2 λ y - π^{2}) + (y') π^{2}}{2 (x^{2} + y^{2})} + \frac{(- {(x')}^{2} - {(y')}^{2} + y π^{2}) (2 x (x') + 2 y (y'))}{2 {(x^{2} + y^{2})}^{2}}\}$

$IRcons ≔ [x^{2} + y^{2} = 1, x (x') + y (y') = 0, λ = - \frac{- {(x')}^{2} - {(y')}^{2} + y π^{2}}{2 (x^{2} + y^{2})}]$

(2.1.1)

Find suitable initial conditions:

>	$ICs := \{y (0) = - 1, D (x) (0) = \frac{1}{10}\}$

$ICs ≔ \{y (0) = −1, D (x) (0) = \frac{1}{10}\}$

(2.1.2)

From the first constraint, we have:

>	$\binom{IRcons [1]}{} \| \binom{}{y (t) = - 1}$

$x^{2} + 1 = 1$

(2.1.3)

>	$ICs := ICs \cup \{x (0) = 0\} &colon;$

From the second constraint, we have:

>	$eval (\binom{convert (IRcons [2], D)}{} \| \binom{}{t = 0}, ICs)$

$- D (y) (0) = 0$

(2.1.4)

>	$ICs := ICs \cup \{D (y) (0) = 0\} &colon;$

And finally from the third constraint, we have:

>	$eval (\binom{convert (IRcons [3], D)}{} \| \binom{}{t = 0}, ICs)$

$λ (0) = \frac{1}{200} + \frac{π^{2}}{2}$

(2.1.5)

>	$ICs := ICs \cup \{\}$

$ICs ≔ \{λ (0) = \frac{1}{200} + \frac{π^{2}}{2}, x (0) = 0, y (0) = −1, D (x) (0) = \frac{1}{10}, D (y) (0) = 0\}$

(2.1.6)

Now, we can obtain a numerical solution for the index reduced system ignoring the constraints.

(Note: We set maxfun=0 to remove the limit on the number of function evaluations.)

>	$dsn := dsolve (IRdsys \cup ICs, numeric, maxfun = 0)$

$dsn ≔ proc (x_rkf45) ... end proc$

(2.1.7)

Consider the value at t=1000:

>	$t1000a := dsn (1000)$

$t1000a ≔ [t = 1000., λ = 4.72577908209305, x = 0.0145466461773924, x' = 0.0883768189980660, y = −1.04481889187418, y' = 0.00110333926104828]$

(2.1.8)

Now consider the constraints at t=1000:

>	$evalf [15] (eval (map (a \to (rhs - lhs) (a), convert (IRcons, D)), convert (t1000a [2 .. - 1], D)))$

$[−0.0918581217322063, −0.000132796612158576, 1.22948907055331 \times 10^{−7}]$

(2.1.9)

This indicates that the constraints are not satisfied.
This is one of the disadvantages of this approach.

Note: Problems arise when small changes to the problem dictate large changes in the solution at a later time (instability).

The Maple DAE extension solvers (rkf45_dae and rosenbrock_dae) are modifications to ODE solvers, and can operate in this mode:

>	$dsn := dsolve (\{op (dsys)\} \cup ICs, numeric, differential = true, projection = false, maxfun = 0)$

$dsn ≔ proc (x_rkf45_dae) ... end proc$

(2.1.10)

>	$t1000a2 := dsn (1000)$

$t1000a2 ≔ [t = 1000., λ = 4.71144373570398, x = −0.0213094042255820, x' = 0.0742773565891539, y = −1.04753334672762, y' = −0.00164448903842064]$

(2.1.11)

>	$evalf [15] (eval (map (a \to (rhs - lhs) (a), convert (IRcons, D)), convert (t1000a2 [2 .. - 1], D)))$

$[−0.0977802032148257, −0.000139850889707689, 1.07473340449360 \times 10^{−7}]$

(2.1.12)

Method of Index Reduction with Removal

This approach is similar to the Method of Index Reduction, but involves the removal of variables that can be solved for.
Advantages and disadvantages are the same as those for the basic index reduction method, but this approach usually provides a better solution.

>	$IRRdsys := eval (\{dsys [1], dsys [2]\}, c3) &semi;$ $IRRcons ≔ [c1, c2]$

$IRRdsys ≔ \{x'' = \frac{(- {(x')}^{2} - {(y')}^{2} + y π^{2}) x}{x^{2} + y^{2}}, y'' = \frac{(- {(x')}^{2} - {(y')}^{2} + y π^{2}) y}{x^{2} + y^{2}} - π^{2}\}$

$IRRcons ≔ [x^{2} + y^{2} = 1, x (x') + y (y') = 0]$

(2.2.1)

>	$IRRICs := remove (has, ICs, λ)$

$IRRICs ≔ \{x (0) = 0, y (0) = −1, D (x) (0) = \frac{1}{10}, D (y) (0) = 0\}$

(2.2.2)

>	$dsn := dsolve (IRRdsys \cup IRRICs, numeric, maxfun = 0)$

$dsn ≔ proc (x_rkf45) ... end proc$

(2.2.3)

>	$t1000b := dsn (1000)$

$t1000b ≔ [t = 1000., x = 0.0126433163397750, x' = 0.0919440494596145, y = −0.998759928311559, y' = 0.00116624938372828]$

(2.2.4)

Examine the two remaining relevant constraints:

>	$evalf [15] (eval (map (a \to (rhs - lhs) (a), convert (IRRcons, D)), convert (t1000b [2 .. - 1], D)))$

$[0.00231875215102262, 2.32544800803330 \times 10^{−6}]$

(2.2.5)

We see that the O(1) error for the general index reduction method has been reduced to an O(10^(-2)) error by using the index reduction with removal method.

Note: The Maple DAE extension solvers can operate in this way, but also have additional error control capabilities. They do not remove the index 1 equation for lambda(t), but instead use it in the numerical integration, and monitor it for possible error.
As a result, the solution is more accurate than previously.

>	$dsn := dsolve (\{op (dsys), op (ICs)\}, numeric, maxfun = 0, projection = false)$

$dsn ≔ proc (x_rkf45_dae) ... end proc$

(2.2.6)

>	$t1000c := dsn (1000)$

$t1000c ≔ [t = 1000., λ = 4.94437510816666, x = 0.0126469052386010, x' = 0.0919391749877821, y = −0.998759726357413, y' = 0.00116651428951809]$

(2.2.7)

>	$evalf [15] (eval (map (a \to (rhs - lhs) (a), convert (IRRcons, D)), convert (t1000c [2 .. - 1], D)))$

$[0.00231906479435229, 2.32145880546108 \times 10^{−6}]$

(2.2.8)

The error in the constraints is now O(10^(-4)).

Method of Index Reduction with Projection

This approach is similar to the Method of Index Reduction, but at each step the solution is projected onto the constraints of the problem.
The magnitude of the projection is monitored as part of the error control strategy.

Advantages

•	Allows us to obtain all constraints directly

•	Enforces constraints, removing the constraint drift problem

Disadvantages

•	Requires a modified solver to add a projection step and error control.

•	Depending on the stability of system, the projection step may introduce errors in the solution.

The Maple rkf45_dae and rosenbrock_dae solvers operate in this mode by default.

>	$dsn := dsolve (\{op (dsys)\} \cup ICs, numeric, maxfun = 0)$

$dsn ≔ proc (x_rkf45_dae) ... end proc$

(2.3.1)

>	$t1000d := dsn (1000)$

$t1000d ≔ [t = 1000., λ = 4.93952122803426, x = −0.00631755178476623, x' = 0.0981561165511285, y = −0.999980043769625, y' = −0.000620117984924986]$

(2.3.2)

Now the constraints are satisfied well within the tolerance at any solution point.

>	$evalf [15] (eval (map (a \to (rhs - lhs) (a), convert (IRcons, D)), convert (t1000d [2 .. - 1], D)))$

$[6.01946159584088 \times 10^{−10}, 7.39595684608879 \times 10^{−10}, −6.84481804569259 \times 10^{−10}]$

(2.3.3)

Direct Solution Approaches

These methods work directly with the original DAE system, and do not require differentiation of the constraint.
They are typically limited to handling DAEs of index 1, 2, or 3.

Advantages

•	Does not require symbolic manipulation of the system

Disadvantages

•	Does not indicate constraints on the initial data, and simply fails if these are not satisfied.

•	Limited to problems with a maximum index.

•	The core solver is usually more complex than for the other approaches.

Maple has the mebdfi solver (Modified Extended Backward Difference Formula Implicit).
This solver has the ability to fully solve implicit DAE systems.
For the pendulum example, we get:

>	$dsn := dsolve (\{op (dsys)\} \cup ICs, numeric, method = mebdfi, maxfun = 0)$

$dsn ≔ proc (x_mebdfi) ... end proc$

(2.4.1)

>	$t1000e := dsn (1000)$

$t1000e ≔ [t = 1000., λ = 4.93950094847455379, x = −0.00630802673444276138, x' = 0.0979430737998922551, y = −0.999980105346507164, y' = −0.000617845741768820308]$

(2.4.2)

>	$evalf [15] (eval (map (a \to (rhs - lhs) (a), convert (IRcons, D)), convert (t1000e [2 .. - 1], D)))$

$[−2.290093 \times 10^{−9}, −5.921958654 \times 10^{−9}, −3.2170225 \times 10^{−7}]$

(2.4.3)

Important: The constraints are all satisfied within tolerance even though they are not explicitly used in the solution process.

Comparison of All Results to a High Accuracy Solution

The pendulum is an excellent example. We have a transformation that converts it to a single nonlinear ODE in theta, which we can solve to high accuracy using our gear method,. This allows us to compare our DAE solutions.

$pode$

$θ'' = - \sin (θ) π^{2}$

(2.5.1)

>	$tr0 := subs (r (t) = 1, tr)$

$tr0 ≔ \{x = \sin (θ), y = - \cos (θ)\}$

(2.5.2)

>	$tr1 := \frac{\partial}{\partial t} tr0$

$tr1 ≔ \{x' = (θ') \cos (θ), y' = (θ') \sin (θ)\}$

(2.5.3)

>	$eval (\binom{convert (tr0 \cup tr1, D)}{} \| \binom{}{t = 0}, ICs)$

$\{−1 = - \cos (θ (0)), 0 = D (θ) (0) \sin (θ (0)), 0 = \sin (θ (0)), \frac{1}{10} = D (θ) (0) \cos (θ (0))\}$

(2.5.4)

From the last two we get theta(0)=0, and so:

>	$eval (, theta (0) = 0) &semi;$

$\{−1 = −1, 0 = 0, \frac{1}{10} = D (θ) (0)\}$

(2.5.5)

Now obtain a high accuracy solution to this problem using the 'gear' method.

$Digits := trunc (evalhf (Digits)) &colon;$ $dsn := dsolve (\{\frac{{&DifferentialD;}^{2}}{&DifferentialD; t^{2}} θ (t) + π^{2} \sin (θ (t)) = 0, θ (0) = 0, D (θ) (0) = \frac{1}{10}\}, numeric, method = gear, abserr = 1. 10^{-13}, relerr = 1. 10^{-13}) &colon;$ $tt ≔ time () &colon;$ $t1000 ≔ dsn (1000) &semi;$ $time () - tt$

$t1000 ≔ [t = 1000., θ = −0.00629176884678665, θ' = 0.0980270351645652]$

$2.203$

(2.5.6)

Convert to x, y coordinates, and obtain the lambda value:

>	$eval (tr0 \cup tr1, t1000 [2 .. - 1])$

$\{x' = 0.0980250949044482, y' = −0.000616759376763182, x = −0.00629172733550273, y = −0.999980206887684\}$

(2.5.7)

>	$ans := \cup evalf (eval (\{c3\},))$

$ans ≔ \{x' = 0.0980250949044482, y' = −0.000616759376763182, λ = 4.93950917526204, x = −0.00629172733550273, y = −0.999980206887684\}$

(2.5.8)

Compare with the basic index reduction solution, which is expected to be poor:

>	$map (a \to (rhs - lhs) (a), eval (t1000a [2 .. - 1], ans))$

$[−0.213730093168983, 0.0208383735128951, −0.00964827590638218, −0.0448386849864990, 0.00172009863781146]$

(2.5.9)

Compare with index reduction with removal, which is expected to be moderately poor:

>	$map (a \to (rhs - lhs) (a), eval (t1000b [2 .. - 1], ans))$

$[0.0189350436752778, −0.00608104544483369, 0.00122027857612550, 0.00178300876049147]$

(2.5.10)

Compare with index reduction with removal and additional error correction (extension methods with projection=false):

>	$map (a \to (rhs - lhs) (a), eval (t1000c [2 .. - 1], ans))$

$[0.00486593290462611, 0.0189386325741037, −0.00608591991666610, 0.00122048053027168, 0.00178327366628127]$

(2.5.11)

Compare with index reduction with projection (extension method in default mode):

>	$map (a \to (rhs - lhs) (a), eval (t1000d [2 .. - 1], ans))$

$[0.0000120527722291541, −0.0000258244492635012, 0.000131021646680274, 1.63118059393064 \times 10^{−7}, −3.35860816180408 \times 10^{−6}]$

(2.5.12)

Compare with direct method (mebdfi method):

>	$map (a \to (rhs - lhs) (a), eval (t1000e [2 .. - 1], ans))$

$[−8.22678748502170 \times 10^{−6}, −0.0000162993989400338, −0.0000820211045558933, 1.01541177133235 \times 10^{−7}, −1.08636500563820 \times 10^{−6}]$

(2.5.13)

In summary, we see that the Maple default rkf45_dae method and mebdfi method give the best agreement with the solution.

Other Problems

The Car Axis System

This multi-body system models the behavior of car suspension on a bumpy road.
It consists of 4 second order ODEs and two index-3 DAEs.

>	$restart &semi;$ $PDEtools [declare] (prime = t, quiet) &colon;$

Auxiliary variables:

>	$yb := \frac{\sin (10 t)}{10} &colon;$ $xb ≔ \sqrt{1 - {yb}^{2}} &colon;$ $ll ≔ \sqrt{{xl (t)}^{2} + {yl (t)}^{2}} &colon;$ $lr ≔ \sqrt{{(xr (t) - xb)}^{2} + {(yr (t) - yb)}^{2}} &colon;$

DAE system:

$eqs := \{\frac{\frac{{&DifferentialD;}^{2}}{&DifferentialD; t^{2}} xl (t)}{2000} = (\frac{1}{2 ll} - 1) \cdot xl (t) + lam1 (t) \cdot xb + 2 \cdot lam2 (t) \cdot (xl (t) - xr (t)), \frac{\frac{{&DifferentialD;}^{2}}{&DifferentialD; t^{2}} yl (t)}{2000} = (\frac{1}{2 ll} - 1) \cdot yl (t) + lam1 (t) \cdot yb + 2 \cdot lam2 (t) \cdot (yl (t) - yr (t)) - \frac{1}{2000}, \frac{\frac{{&DifferentialD;}^{2}}{&DifferentialD; t^{2}} xr (t)}{2000} = (\frac{1}{2 lr} - 1) \cdot (xr (t) - xb) - 2 \cdot lam2 (t) \cdot (xl (t) - xr (t)), \frac{\frac{{&DifferentialD;}^{2}}{&DifferentialD; t^{2}} yr (t)}{2000} = (\frac{1}{2 lr} - 1) \cdot (yr (t) - yb) - 2 \cdot lam2 (t) \cdot (yl (t) - yr (t)) - \frac{1}{2000}, 0 = xl (t) xb + yl (t) yb, 0 = {(xl (t) - xr (t))}^{2} + {(yl (t) - yr (t))}^{2} - 1\}$

$eqs ≔ \{0 = \frac{xl (t) \sqrt{- {\sin (10 t)}^{2} + 100}}{10} + \frac{yl (t) \sin (10 t)}{10}, 0 = {(xl (t) - xr (t))}^{2} + {(yl (t) - yr (t))}^{2} - 1, \frac{xl''}{2000} = (\frac{1}{2 \sqrt{{xl (t)}^{2} + {yl (t)}^{2}}} - 1) xl (t) + \frac{lam1 (t) \sqrt{- {\sin (10 t)}^{2} + 100}}{10} + 2 lam2 (t) (xl (t) - xr (t)), \frac{xr''}{2000} = (\frac{1}{2 \sqrt{{(xr (t) - \frac{\sqrt{- {\sin (10 t)}^{2} + 100}}{10})}^{2} + {(yr (t) - \frac{\sin (10 t)}{10})}^{2}}} - 1) (xr (t) - \frac{\sqrt{- {\sin (10 t)}^{2} + 100}}{10}) - 2 lam2 (t) (xl (t) - xr (t)), \frac{yl''}{2000} = (\frac{1}{2 \sqrt{{xl (t)}^{2} + {yl (t)}^{2}}} - 1) yl (t) + \frac{lam1 (t) \sin (10 t)}{10} + 2 lam2 (t) (yl (t) - yr (t)) - \frac{1}{2000}, \frac{yr''}{2000} = (\frac{1}{2 \sqrt{{(xr (t) - \frac{\sqrt{- {\sin (10 t)}^{2} + 100}}{10})}^{2} + {(yr (t) - \frac{\sin (10 t)}{10})}^{2}}} - 1) (yr (t) - \frac{\sin (10 t)}{10}) - 2 lam2 (t) (yl (t) - yr (t)) - \frac{1}{2000}\}$

(3.1.1)

Initial conditions and variables:

>	$ini := \{xl (0) = 0, D (xl) (0) = - \frac{1}{2}, yl (0) = \frac{1}{2}, D (yl) (0) = 0, xr (0) = 1, D (xr) (0) = - \frac{1}{2}, yr (0) = \frac{1}{2}, D (yr) (0) = 0, lam1 (0) = 0, lam2 (0) = 0\} &semi;$ $vrs := [xl (t), yl (t), xr (t), yr (t), lam1 (t), lam2 (t)] &colon;$

$ini ≔ \{lam1 (0) = 0, lam2 (0) = 0, xl (0) = 0, xr (0) = 1, yl (0) = \frac{1}{2}, yr (0) = \frac{1}{2}, D (xl) (0) = - \frac{1}{2}, D (xr) (0) = - \frac{1}{2}, D (yl) (0) = 0, D (yr) (0) = 0\}$

(3.1.2)

Exact solution (solution computed via gear with parameter removal and 1e-22 accuracy):

$yans := [0.0493455784275240921227327851163, - 0.0770583684035920842786873048316, 0.496989460230008106764570177300, 0.00744686659206841657568543179878, 1.04174252488542611518584788701, 0.0175568157535417366330123387591, 0.373911027265365819359535268815, 0.770341043779601063189353416447, - 0.00473688659085332651484983657657, - 0.00110468033125956583962543762007] &colon;$

>	$tt := time () &colon;$ $ds1 := dsolve (\{op (eqs), op (ini)\}, vrs, numeric) &semi;$ $time () - tt$

$ds1 ≔ proc (x_rkf45_dae) ... end proc$

$1.125$

(3.1.3)

Solution with default tolerances:

>	$tt := time () &colon;$ $dsol1 := ds1 (3) &semi;$ $time () - tt &semi;$

$dsol1 ≔ [t = 3., xl (t) = 0.0493455779014815, xl' = −0.0770590607488527, yl (t) = 0.496989454931791, yl' = 0.00743988527830008, xr (t) = 1.04174242180064, xr' = 0.0175614035187829, yr (t) = 0.373910194846429, yr' = 0.770371407475209, lam1 (t) = −0.00473687080377615, lam2 (t) = −0.00110467301491489]$

$0.110$

(3.1.4)

Compare:

>	$evalf [16] (map (rhs, dsol1 [2 .. - 1]) - yans)$

$[−5.2604258 \times 10^{−10}, −6.9234526065 \times 10^{−7}, −5.2982167 \times 10^{−9}, −6.981313768339 \times 10^{−6}, −1.03084791 \times 10^{−7}, 4.58776524117 \times 10^{−6}, −8.324189364 \times 10^{−7}, 0.0000303636956078, 1.5787077181 \times 10^{−8}, 7.316344672 \times 10^{−9}]$

(3.1.5)

Solution with tight tolerances:

>	$tt := time () &colon;$ $ds2 := dsolve (\{op (eqs), op (ini)\}, vrs, numeric, abserr = 1. 10^{-11}, relerr = 1. 10^{-11}, maxfun = 0) &semi;$ $time () - tt &semi;$ $tt ≔ time () &colon;$ $dsol2 := ds2 (3) &semi;$ $time () - tt &semi;$ $evalf [16] (map (rhs, dsol2 [2 .. - 1]) - yans)$

$ds2 ≔ proc (x_rkf45_dae) ... end proc$

$0.047$

$dsol2 ≔ [t = 3., xl (t) = 0.0493455784275599, xl' = −0.0770583684067605, yl (t) = 0.496989460230369, yl' = 0.00744686656072501, xr (t) = 1.04174252488285, xr' = 0.0175568158504510, yr (t) = 0.373911027244694, yr' = 0.770341044422276, lam1 (t) = −0.00473688659029054, lam2 (t) = −0.00110468033098021]$

$0.750$

$[3.582 \times 10^{−14}, −3.16846 \times 10^{−12}, 3.606 \times 10^{−13}, −3.1343407 \times 10^{−11}, −2.575 \times 10^{−12}, 9.690928 \times 10^{−11}, −2.06719 \times 10^{−11}, 6.426745 \times 10^{−10}, 5.62785 \times 10^{−13}, 2.79361 \times 10^{−13}]$

(3.1.6)

Trajectory of x1,y1:

>	$plots [odeplot] (ds1, [xl (t), yl (t)], 0 .. 3, numpoints = 500)$

Trajectory of x2,y2:

>	$plots [odeplot] (ds1, [xr (t), yr (t)], 0 .. 3, numpoints = 500)$

The Double Pendulum

The double pendulum is a pendulum with an additional weight attached to the first pendulum bob by a string.
It is a system with 4 second order ODEs and two index-3 DAEs.

(x1(t),y1(t)) are the coordinates of the first pendulum bob, and (x2(t),y2(t)) are the coordinates of the second, which is attached to the first.

>	$restart &semi;$ $PDEtools [declare] (prime = t, quiet) &colon;$

Set masses, lengths of strings, and g:

>	$m1 ≔ 1 &colon;$ $l1 ≔ 1 &colon;$ $m2 ≔ 1 &colon;$ $l2 ≔ \frac{1}{2} &colon;$ $g ≔ 9.8 &colon;$

DAE system and initial conditions:

$sys := \{m1 \cdot (\frac{{&DifferentialD;}^{2}}{&DifferentialD; t^{2}} x1 (t)) + 2 \cdot λ1 (t) \cdot x1 (t) + 2 \cdot λ2 (t) \cdot (x1 (t) - x2 (t)), m1 \cdot (\frac{{&DifferentialD;}^{2}}{&DifferentialD; t^{2}} y1 (t)) + m1 \cdot g + 2 \cdot λ1 (t) \cdot y1 (t) + 2 \cdot λ2 (t) \cdot (y1 (t) - y2 (t)), m2 \cdot (\frac{{&DifferentialD;}^{2}}{&DifferentialD; t^{2}} x2 (t)) - 2 \cdot λ2 (t) \cdot (x1 (t) - x2 (t)), m2 \cdot (\frac{{&DifferentialD;}^{2}}{&DifferentialD; t^{2}} y2 (t)) + m2 \cdot g - 2 \cdot λ2 (t) \cdot (y1 (t) - y2 (t)), x1 {(t)}^{2} + {y1 (t)}^{2} - {l1}^{2}, {(x1 (t) - x2 (t))}^{2} + {(y1 (t) - y2 (t))}^{2} - {l2}^{2}\} &semi;$ $ics := \{x1 (0) = 0, D (x1) (0) = - \frac{5}{2}, y1 (0) = - l1, D (y1) (0) = 0, x2 (0) = 0, D (x2) (0) = 3, y2 (0) = - l1 - l2, D (y2) (0) = 0\} &colon;$

$sys ≔ \{x2'' - 2 λ2 (t) (x1 (t) - x2 (t)), {(x1 (t) - x2 (t))}^{2} + {(y1 (t) - y2 (t))}^{2} - \frac{1}{4}, {x1 (t)}^{2} + {y1 (t)}^{2} - 1, x1'' + 2 λ1 (t) x1 (t) + 2 λ2 (t) (x1 (t) - x2 (t)), y2'' + 9.8 - 2 λ2 (t) (y1 (t) - y2 (t)), y1'' + 9.8 + 2 λ1 (t) y1 (t) + 2 λ2 (t) (y1 (t) - y2 (t))\}$

(3.2.1)

Obtain the trajectory for the second mass:

>	$dsn := dsolve (sys \cup ics, numeric) &semi;$ $plots [odeplot] (dsn, [x2 (t), y2 (t)], 0 .. 5, numpoints = 1000)$

$dsn ≔ proc (x_rkf45_dae) ... end proc$

Construct an animation of the system, where tf is the final time and nfr is the number of frames:

$tf ≔ 5 &colon;$ $nfr ≔ 201 &colon;$ $pl ≔ NULL &colon;$ $ps ≔ 0.02 &colon;$ $for i to nfr do dsv := dsn (\frac{5 (i - 1)}{nfr - 1}) &colon; x1v, x2v, y1v, y2v := seq (rhs ({dsv}_{j}), j = [4, 6, 8, 10]) &colon; pl := pl, [CURVES ([[0, 0], [x1v, y1v], [x2v, y2v]], COLOR (RGB, 1, 0, 0)), POLYGONS ([[x1v + ps, y1v + ps], [x1v + ps, y1v - ps], [x1v - ps, y1v - ps], [x1v - ps, y1v + ps], [x1v + ps, y1v + ps]], [[x2v + ps, y2v + ps], [x2v + ps, y2v - ps], [x2v - ps, y2v - ps], [x2v - ps, y2v + ps], [x2v + ps, y2v + ps]], COLOR (RGB, 0, 0, 1))] &colon; end do &colon;$

>	$PLOT (ANIMATE (pl), AXESSTYLE (NONE))$

Illustration of Constraint Drift

Consider the simple pendulum system and the worst rkf45_dae solution (differential=true, projection=false):

$restart &semi;$ $PDEtools [declare] (\{x (t), y (t), λ (t)\}, prime = t, quiet) &colon;$ $dsys := [\frac{{&DifferentialD;}^{2}}{&DifferentialD; t^{2}} x (t) = - 2 \cdot λ (t) \cdot x (t), \frac{{&DifferentialD;}^{2}}{&DifferentialD; t^{2}} y (t) = - 2 \cdot λ (t) \cdot y (t) - π^{2}, {x (t)}^{2} + {y (t)}^{2} = 1, y (0) = - 1, D (x) (0) = \frac{1}{10}, D (y) (0) = 0, x (0) = 0, λ (0) = \frac{1}{200} + \frac{1 π^{2}}{2}]$

$dsys ≔ [x'' = - 2 λ x, y'' = - 2 λ y - π^{2}, x^{2} + y^{2} = 1, y (0) = −1, D (x) (0) = \frac{1}{10}, D (y) (0) = 0, x (0) = 0, λ (0) = \frac{1}{200} + \frac{π^{2}}{2}]$

(3.3.1)

>	$dsn := dsolve (dsys, numeric, differential = true, projection = false, maxfun = 0)$

$dsn ≔ proc (x_rkf45_dae) ... end proc$

(3.3.2)

Since this is a pendulum, the trajectory (x(t),y(t)) should trace a portion of a circle repeatedly. Constraint drift causes a deviation this behavior.

>	$plots [odeplot] (dsn, [x (t), y (t)], 0 .. 100, numpoints = 2000, frames = 100)$

The default solution mode works much better:

>	$dsn := dsolve (dsys, numeric, maxfun = 0)$

$dsn ≔ proc (x_rkf45_dae) ... end proc$

(3.3.3)

>	$plots [odeplot] (dsn, [x (t), y (t)], 0 .. 100, numpoints = 2000, frames = 100)$

Return to Index for Example Worksheets

Download Help Document

Maple

Maple Add-Ons

Student Success Platform

Improving Retention Rates

MapleSim

MapleSim Add-Ons

Systems Engineering

Consulting Services

Maple T.A. and Möbius

Education

Industries

Automotive and Aerospace

Robotics

Machine Design & Industrial Automation

Other

Application Areas

Product Pricing

Purchasing

Institutional Student Licensing

Maplesoft Elite Maintenance (EMP)

Support

Product Training

Online Product Help

Webinars & Events

Publications

Content Hubs

Examples & Applications

Community

About Maplesoft

Media Center

User Community

Contact

Online Help

All Products Maple MapleSim

Maple

Powerful math software that is easy to use

Maple Add-Ons

Student Success Platform

Improving Retention Rates

MapleSim

Advanced System Level Modeling

MapleSim Add-Ons

Systems Engineering

Consulting Services

Maple T.A. and Möbius

Education

Industries

Automotive and Aerospace

Robotics

Machine Design & Industrial Automation

Other

Application Areas

Product Pricing

Purchasing

Institutional Student Licensing

Maplesoft Elite Maintenance (EMP)

Support

Product Training

Online Product Help

Webinars & Events

Publications

Content Hubs

Examples & Applications

Community

About Maplesoft

Media Center

User Community

Contact

Online Help

All Products Maple MapleSim