MAPLE API

Diff

The Diff command extends of the MAPLE built-in diff command to work with function as derivation argument. For example, if you try to differentiate the following function with the built-in diff command like:

> fun := sin(x(t)):
> diff(sin(x(t)), x(t));

an error is raised because MAPLE expect a symbol and not a function as second argument. The Diff command instead can be used to differentiate a function with respect to a function argument:

> fun := sin(x(t)):
> Diff(sin(x(t)),x(t));

    cos(x(t))

this returning the correct answer.

Jacobian

This command is the natural extension of the Diff command to compute the Jacobian matrix of a function (refer to the Diff command documentation for more details). For example, if one wants to compute the Jacobian matrix of the following function/map fun with respect to the function/parameters pars the syntax is the following:

> Jacobian(fun, pars);

A realistic example is the following:

> fun  := <x^2+y(t), y(t)*cos(x*y(t))>:
> pars := [x, y(t)]:
> Jacobian(fun, pars);

    [       2*x                         1               ]
    [-y(t)^2*sin(x*y(t))  cos(x*y(t))-y(t)*x*sin(x*y(t))]

KernelBuild

Give a (potentially rectangular) matrix \(\mathbf{E}\), the command KernelBuild computes the matrices \(\mathbf{K}\) and \(\mathbf{L}\) such that:

\[\begin{split}\begin{array}{l} \mathbf{K}\mathbf{E} = \mathbf{0} \\[0.5em] \mathbf{L}\mathbf{E} \quad \textrm{is full rank} \\[0.5em] \mathbf{M} = \left[\begin{array}{c} \mathbf{L} \\ \mathbf{K} \end{array}\right] \quad \textrm{is non singular}. \end{array}\end{split}\]

This command also returns the rank of the matrix \(\mathbf{E}\). The algorithm use Gaussian elimination so that \(\mathbf{E}\) can contain symbolic expressions.

A usage example is the following:

> K, L, r := KernelBuild(E);

where E can be a rectangular matrix.

DAE_separate_algebraic

Given a DAE in the form

\[\mathbf{E}(\mathbf{x},t) \mathbf{x}' = \mathbf{g}(\mathbf{x},t)\]

using KernelBuild transform the DAE to

\[\begin{split}\left\{\begin{array}{rcl} \mathbf{E}_1(\mathbf{x},t) \mathbf{x}' &=& \mathbf{g}_1(\mathbf{x},t) \\[1em] \mathbf{0} &=& \mathbf{g}_2(\mathbf{x},t) \end{array}\right.\end{split}\]

separating the algebraic part into \(\mathbf{g}_2(\mathbf{x},t)\).

Usage:

> E1, G1, G2, r := DAE_separate_algebraic( E, G ): # r = rank or E

Notice that the routine return also the rank of the matrix \(\mathbf{E}\).

There is also a function DAE_separate_algebraic_bis which do the same job when the DAE is passed as a list of differential equations. In this case you must also pass the list of differential variables to transform (internally) to the form \(\mathbf{E}(\mathbf{x},t) \mathbf{x}' = \mathbf{g}(\mathbf{x},t)\)

> E1, G1, G2, r := DAE_separate_algebraic_bis( EQNS, DVARS ): # r = rank or E

DAE_make_semi_explicit

Given DAE passed as a list of differential equations build a new DAE in semi explicit form. The user musty pass

• A list the the DAE system

• A list with the variables (functions) of the DAE

> ODE, DVARS, AVARS, ALG := DAE_make_semi_explicit( DAE, vars )

After the reduction you have

• ODE the ODE part \(x' = f(x,y)\) of the DAE

• DVARS the list of function that appers as derivative \(x(t)\)

• AVARS the list of function that DO NOT appers as derivative \(y(t)\)

• ALG the algebraic part \(0 = g(x,y)\) of the DAE

In the process of semi-explicit formation some new variable may be created. Moreover ALG part can contain trivial equations that can be manually solved by the user.

For exmaple the Pendulum DAE

\[\begin{split}\left\{ \begin{array}{l} x' = u \\ y' = v \\ u' + \lambda x = 0 \\ v' + \lambda y = -mg \\ x^2+y^2=1 \end{array} \right.\end{split}\]

is transformed to

ODE:

\[\begin{split}\left\{ \begin{array}{l} x' = u \\ y' = v \\ u' = \dot{u} \\ u' = \dot{v} \end{array} \right.\end{split}\]

ALG

\[\begin{split}\left\{ \begin{array}{l} \dot{u} + \lambda x = 0 \\ \dot{v} + \lambda y + mg = 0 \\ x^2+y^2-1 = 0 \end{array} \right.\end{split}\]

DVARS

\[[ x(t), y(t), u(t), v(t) ]\]

AVARS

\[[ \dot{u}(t), \dot{v}(t), \lambda(t) ]\]

For a non trivial usare of DAE_make_semi_explicit lokk at DAE-toolbox-usare-2.mw

DAE_reduce_index_by_1

Given a DAE in the form (you che put in this form using DAE_separate_algebraic)

\[\begin{split}\left\{\begin{array}{rcl} \mathbf{E}_1(\mathbf{x},t) \mathbf{x}' &=& \mathbf{g}_1(\mathbf{x},t) \\[1em] \mathbf{0} &=& \mathbf{a}_1(\mathbf{x},t) \end{array}\right.\end{split}\]

Tranform to a new one

\[\begin{split}\left\{\begin{array}{rcl} \mathbf{E}_2(\mathbf{x},t) \mathbf{x}' &=& \mathbf{g}_2(\mathbf{x},t) \\[1em] \mathbf{0} &=& \mathbf{a}_2(\mathbf{x},t) \end{array}\right.\end{split}\]

That has index reduced by one. The command usage is the following

> E2, G2, A2, r := DAE_reduce_index_by_1( E1, G1, A1, Dvars );

where

• E1 is the matrix \(\mathbf{E}_1(\mathbf{x},t)\)

• G1 is the vector \(\mathbf{g}_1(\mathbf{x},t)\)

• A1 is the vector \(\mathbf{a}_1(\mathbf{x},t)\) of the algebraic constraints

• Dvars is the list of the differential variable \(\mathbf{x}'(t)\)

and

• E2 is the matrix \(\mathbf{E}_2(\mathbf{x},t)\)

• G2 is the vector \(\mathbf{g}_2(\mathbf{x},t)\)

• A2 is the vector \(\mathbf{a}_2(\mathbf{x},t)\) of the NEW algebraic constraints

• r the rank of the output matrix \(\mathbf{E}_2(\mathbf{x},t)\)

If the reduced DAE is an ODE A2 is empty and r is equal to the number of equations.

Library has also the functions:

• DAE_reduce_index_by_1_full( E, G, Dvars ) Do not need to previously separate algebraic part, is done internally.

• DAE_reduce_index_by_1_full2proc( EQS, Dvars ) Do not need to put in the form E x’ = G is done internally.

DAE3_get_ODE_and_invariants

Given an index-3 DAE of the form

\[\begin{split}\mathrm{DAE}: \left\{ \begin{array}{l} \mathbf{q}' = \mathbf{v} \\[0.5em] \mathbf{M}(\mathbf{q},\mathbf{v},t) \mathbf{v}' + \mathbf{\Phi}_q(\mathbf{q},t)^T\boldsymbol{\lambda} = \mathbf{f}(\mathbf{q},\mathbf{v},t) \\[0.5em] \mathbf{\Phi}(\mathbf{q},t) = \mathbf{0} \end{array} \right.\end{split}\]

Trasform to semi-explicit DAE

\[\begin{split}\mathrm{ODE}: \left\{ \begin{array}{l} \mathbf{q}' = \mathbf{v} \\[0.5em] \mathbf{v}' = \dot{\mathbf{v}} \end{array} \right. \qquad \mathrm{ALG}: \left\{ \begin{array}{l} \mathbf{M}(\mathbf{q},\mathbf{v},t) \dot{\mathbf{v}} + \mathbf{\Phi}_q(\mathbf{q},t)^T\boldsymbol{\lambda} = \mathbf{f}(\mathbf{q},\mathbf{v},t) \\[0.5em] \mathbf{\Phi}(\mathbf{q},t) = \mathbf{0} \end{array} \right.\end{split}\]

Then build first and second derivative of the constraints \(\mathbf{\Phi}(\mathbf{q},t)\):

\[\mathbf{a}(\mathbf{q},\mathbf{v},t)=\dfrac{\mathrm{d}}{\mathrm{d}t}\mathbf{\Phi}(\mathbf{q},t) = \mathbf{\Phi}_q(\mathbf{q},t)\mathbf{v}+ \mathbf{\Phi}_t(\mathbf{q},t)\]

and

\[\dfrac{\mathrm{d}}{\mathrm{d}t}\mathbf{a}(\mathbf{q},\mathbf{v},t)= \mathbf{\Phi}_q(\mathbf{q},t)\dot{\mathbf{v}}-\mathbf{g}(\mathbf{q},\mathbf{v},t)\]

where

\[\mathbf{g}(\mathbf{q},\mathbf{v},t) = -\dfrac{\mathrm{d}}{\mathrm{d}t}\mathbf{a}(\mathbf{q},\mathbf{v},t)|_{\mathbf{v}=\mathrm{fixed}}\]

USAGE:

res := DAE3_get_ODE_and_invariants( Mass, Phi, f, qvars, vvars, lvars )

where

Parameter correspondence

Mass	\(\mathbf{M}(\mathbf{q},\mathrm{v},t)\)
Phi	\(\mathbf{\Phi}(\mathbf{q},t)\)
f	\(\mathbf{f}(\mathbf{q},\mathbf{v},t)\)
qvars	\(\mathbf{q}=[q_1(t),q_2(t),\ldots,q_n(t)]\)
vvars	\(\mathbf{v}=[v_1(t),v_2(t),\ldots,v_n(t)]\)
lvars	\(\boldsymbol{\lambda}=[\lambda_1(t),\ldots,\lambda_m(t)]\)

the result res is a maple table that contains

Table contents

res["PVARS"]	The position states \([q_1(t),q_2(t),\ldots,q_n(t)]\)
res["VVARS"]	The velocity states \([v_1(t),v_2(t),\ldots,v_n(t)]\)
res["LVARS"]	The Lagrange multipliers \([\lambda_1(t),\lambda_2(t),\ldots,\lambda_m(t)]\)
res["VDOT"]	The added algebraic states \([\dot{v}_1(t),\dot{v}_2(t),\ldots,\dot{v}_n(t)]\)
res["ODE_RHS"]	The r.h.s for the ODE part (complete)
res["ODE_POS"]	The r.h.s for the ODE part: position equations
res["ODE_VEL"]	The r.h.s for the ODE part: velocity equations
res["Phi"]	The vector ot the constraints \(\mathbf{\Phi}(\mathbf{q},t)\)
res["Phi_P"]	Partial derivative of the constraints \(\mathbf{\Phi}_q(\mathbf{q},t)\)
res["A"]	\(\mathbf{a}(\mathbf{q},\mathbf{v},t)=\mathbf{\Phi}_q(\mathbf{q},t)\dot{\mathbf{v}}-\mathbf{b}(\mathbf{q},\mathbf{v},t)\)
res["A_rhs"]	\(-\mathbf{\Phi}_t(\mathbf{q},t)\)
res["g"]	\(\mathbf{g}(\mathbf{q},\mathbf{v},t)\)
res["bigVAR"]	\([\dot{v}_1(t),\dot{v}_2(t),\ldots,\dot{v}_n(t),\lambda_1(t),\lambda_2(t),\ldots,\lambda_m(t)]\),
res["bigM"]	\(\left[\begin{array}{cc}\mathbf{M}(\mathbf{q},\mathbf{v},t) & \mathbf{\Phi}_q(\mathbf{q},t)^T \\ \mathbf{\Phi}_q(\mathbf{q},t) & \mathbf{0}\end{array}\right]\)
res["bigRHS"]	\(\left[\begin{array}{c}\mathbf{f}(\mathbf{q},\mathbf{v},t) \\ \mathbf{b}(\mathbf{q},\mathbf{v},t)\end{array}\right]\)

DAE3_get_ODE_and_invariants_full

The extended version of the call DAE3_get_ODE_and_invariants

res := DAE3_get_ODE_and_invariants_full( Mass, Phi, f, qvars, vvars, lvars )

return the same table of DAE3_get_ODE_and_invariants with in addition

Table contents

res["bigETA"]	\(\boldsymbol{\eta}(\mathbf{q},\mathbf{v},\boldsymbol{\mu},t)=\mathbf{M}(\mathbf{q},\mathbf{v},t)\boldsymbol{\mu}\) where \(\boldsymbol{\mu}=[\mu_1,\mu_2,\ldots,\mu_n]^T\)
res["JbigETA"]	\([\boldsymbol{\eta}_{\mathbf{q}}(\mathbf{q},\mathbf{v},\boldsymbol{\mu},t),\boldsymbol{\eta}_{\mathbf{v}}(\mathbf{q},\mathbf{v},\boldsymbol{\mu},t)]\)
res["JbigRHS"]	\(\left[\begin{array}{cc}\mathbf{f}_{\mathbf{q}}(\mathbf{q},\mathbf{v},t) & \mathbf{f}_{\mathbf{v}}(\mathbf{q},\mathbf{v},t) \\ \mathbf{b}_{\mathbf{q}}(\mathbf{q},\mathbf{v},t) & \mathbf{b}_{\mathbf{v}}(\mathbf{q},\mathbf{v},t) \end{array}\right]\)

DAE3_get_ODE_and_invariants_stab

The extended version of the call DAE3_get_ODE_and_invariants_full

res := DAE3_get_ODE_and_invariants_full( Mass, Phi, f, qvars, vvars, lvars )

return the same table of DAE3_get_ODE_and_invariants_full with in addition the stabilized constraints with Baumgarte technique.

Table contents

res["h"]	\(\mathbf{h}(\mathbf{q},\mathbf{v},t)=\mathbf{g}(\mathbf{q},\mathbf{v},t)-2\eta\mathbf{a}(\mathbf{q},\mathbf{v},t)-\omega^2\mathbf{\Phi}(\mathbf{q},t)\)
res["bigRHS_stab"]	\(\left[\begin{array}{c}\mathbf{f}(\mathbf{q},\mathbf{v},t) \\ \mathbf{h}(\mathbf{q},\mathbf{v},t)\end{array}\right]\)
res["JbigRHS_stab"]	The jacobian of res["bigRHS_stab"]

F_TO_MATLAB

JF_TO_MATLAB