Numerical solution of an ordinary differential equation on a homogeneous manifold
MuPAD® notebooks are not recommended. Use MATLAB® live scripts instead.
MATLAB live scripts support most MuPAD functionality, though there are some differences. For more information, see Convert MuPAD Notebooks to MATLAB Live Scripts.
numeric::odesolveGeometric(f
,t_{0} .. t
,Y_{0}
, <LieGroupAction = LAMBDA
>, <method
>, <RelativeError = tol
>, <Stepsize = h
>, <Alldata = n
>) numeric::odesolveGeometric(t0 .. t
,f
,Y_{0}
, <LieGroupAction = LAMBDA
>, <method
>, <RelativeError = tol
>, <Stepsize = h
>, <Alldata = n
>)
numeric::odesolveGeometric(f, t_0..t, Y_0 )
approximates
the solution of
,
where f(t, Y(t)) returns n×n matrices
and
.
numeric::odesolveGeometric
is a "geometrical
integrator" for ordinary differential equations on homogeneous
manifolds embedded in the space of n×m matrices.
The call numeric::odesolveGeometric(f, t_0..t, Y_0
)
returns a numerical approximation of the solution Y(t) of
the first order differential equation (dynamical system)
with . Here, Y(t) is a curve of n×m matrices ( or vectors in or ). The function f must produce n×n matrices as return values.
The following geometrical feature is preserved by the numerical
solution: If the matrices produced by f lie
in some Lie subalgebra g of
the n×n matrices,
then, within the numerical working precision, the approximation produced
by numeric::odesolveGeometric
stays on the homogeneous
manifold
,
where G is
the matrix Lie group of g.
As an introductory example, consider the ODE
,
where Y is
a vector in
and f
produces
skew symmetric matrices. The solution lies on the orbit of the orthogonal
group SO(n) generated
by the skew symmetric matrices through the initial point
.
Here, SO(n) acts
on
by
standard matrix multiplication. The homogeneous manifold given by
the orbit of SO(n) through Y_{0} is
the sphere
.
Using standard numerical schemes, the numerical solution drifts
away from this manifold in the course of the integration. The geometrical
"Lie group" integrator numeric::odesolveGeometric
,
however, produces a numerical solution that stays on this manifold,
preserving the invariants of the group action. In this case, the invariant
is
preserved numerically. See Example 1.
With Y(t) = G(t) Y_{0}, the matrix ODE
is solved on the space of the complex n×n matrices (1_{n, n} is the identity matrix). Following MuntheKaas [1], the ansatz reduces a time step for the ode above to the solution of the matrix ODE
,
where
and [u, f]
= u f  f u is
the commutator on the Lie algebra of n×n matrices.
In each step, the ODE for u is
solved numerically in a classical way by the RungeKutta scheme specified
by the parameter method
. Finally, numeric::odesolveGeometric
performs
the time step
by
computing
.
If the matrices produced by f(t, Y) lie in a Lie subalgebra g of the n×n matrices, then the numerical solution u also lies in g. The matrix is an element of the corresponding Lie group, and Y = G Y_{0} lies on the orbit of the Lie group through the initial value Y_{0}. Thus, the geometrical invariants of the homogeneous manifold are preserved in the course of the numerical integration.
The input data t_{0} and t must
not contain symbolic objects which cannot be converted to floating
point values via float
.
Numerical expressions such as
,
etc.
are accepted.
The initial condition Y_{0} defines the space in which the homogeneous manifold containing the solution is embedded.
If Y_{0} is
a list with n entries
or a 1dimensional array array(1..n)
, then the
solution Y(t) consists
of vectors from a submanifold of
or
,
respectively.
If Y_{0} is
specified as a 2dimensional array(1..n, 1..m)
or
as a matrix of the corresponding dimension generated by the function matrix
, then the solution Y(t) consists
of matrices from a submanifold of the space
of n×m matrices.
Internally, 2dimensional n×m arrays
are used to represent the points on the manifold where m =
1 for vectors in
or
.
It is recommended to specify Y_{0} in
the form array(1..n, 1..m)
in order to avoid the
overhead of internal conversions.
The "vector field" f
defining
the dynamical system
must
be represented by a procedure with two input parameters: the scalar
time t and
the matrix or vector Y.
Internally, f is
called with real floatingpoint values t and
matrices/vectors Y of
the same domain type
as the initial condition Y_{0}.
The procedure f
has to return an n×n matrix
either as an array(1..n, 1..n)
or as a corresponding
matrix object of category Cat::Matrix
(generated
by the function matrix
).
It is recommended that the procedure returns an array of the
type array(1..n, 1..n)
. This avoids the overhead
of internal conversions.
The return value of f
may contain numerical
expressions such as π,
etc.
However, all values must be convertible to real or complex floating
point numbers by float
.
Autonomous systems, where f(t, Y) does
not depend on t,
must be represented by a procedure with two arguments t
and Y
,
too.
The optional arguments method
, RelativeError
= tol
, and Stepsize = h
determine how
the ODE
is
solved. They correspond to the methods of the classical ODE solver numeric::odesolve
.
The numerical precision is controlled by the global variable DIGITS
:
an adaptive control of the step size keeps local relative discretization
errors below
,
unless a different tolerance is specified via the option RelativeError
= tol
. The error control may be switched off by specifying
a fixed Stepsize = h
.
Note: Only local errors are controlled by the adaptive mechanism. No control of the global error is provided! 
With Y := t > numeric::odesolveGeometric(f, t_0..t,
Y_0)
, the numerical solution can be repesented by a MuPAD^{®} function:
the call Y(t)
will start the numerical integration
from t_{0} to t.
Classical integration preserves the geometrical invariants up
to the relative precision of the solution whereas the geometrical
integrator preserves the invariants independent of tol up
to the working precision set by DIGITS
: departure from the homogeneous
manifold is a pure roundoff effect.
numeric::odesolveGeometric
is useful when
a tolerance tol much
largerthan
is
specified by RelativeError = tol
. For small tolerances,
you may consider to use the classical solver numeric::odesolve
instead.
Since classical integration is significantly faster, larger
values of DIGITS
and
smaller tolerances for the discretization error may be used in numeric::odesolve
. Depending
on the concrete problem, this may lead to better results than produced
by numeric::odesolveGeometric
.
The function is sensitive to the environment variable DIGITS
,
which determines the numerical working precision.
We consider the initial value problem
for with fixed parameters , J_{2} = 1, J_{3} = 2. Writing this ODE as
,
it is clear that the solution is restricted to the orbit of the orthogonal group SO(3) through the initial point (f_{1} produces skew symmetric matrices). The invariant of this action is , i.e., the solution is restricted to a sphere. Writing the ODE as
,
it is clear that the solution is also restricted to the orbit of the "Jorthogonal" group SO(J, 3) through the initial point. This group consists of matrices G satisfying G^{T} J G = J, where J = diag(J_{1}, J_{2}, J_{3}). The invariant of this group action is , i.e., the solution is restricted to an ellipsoid.
We consider the first representation and compute a numerical solution that is restricted to a sphere:
f1 := proc(t, Y) begin array(1..3, 1..3, [ [ 0 , J3*Y[3], J2*Y[2]], [ J3*Y[3], 0 , J1*Y[1]], [J2*Y[2], J1*Y[1], 0 ]]) end_proc: J1 := 1/2: J2 := 1: J3 := 2: tol := 10^(2): Gsolve:= (f, t0_t, Y0) > numeric::odesolveGeometric(f, t0_t, Y0, RelativeError = tol): Y(0) := [1.0, 1.0, 1.0]; Y(1) := Gsolve(f1, 0..1, Y(0)); Y(2) := Gsolve(f1, 1..2, Y(1)); Y(3) := Gsolve(f1, 2..3, Y(2)); Y(4) := Gsolve(f1, 3..4, Y(3)); Y(5) := Gsolve(f1, 4..5, Y(4))
The invariant H_{1} is
preserved numerically up to the working precision set by DIGITS
:
H1 := Y > Y[1]^2 + Y[2]^2 + Y[3]^2: H1(Y(i))  H1(Y(0)) $ i = 1..5
The invariant H_{2} is
only preserved within the relative precision of the solution set by
the option RelativeError = tol
:
H2 := Y > J1*Y[1]^2 + J2*Y[2]^2 + J3*Y[3]^2: H2(Y(i))  H2(Y(0)) $ i = 1..5
Now, we solve the ODE using the second representation:
f2 := proc(t, Y) begin array(1..3, 1..3, [ [ 0 , J2*Y[3], J3*Y[2]], [J1*Y[3], 0 , J3*Y[1]], [ J1*Y[2], J2*Y[1], 0 ]]) end_proc: Y(0) := [1.0, 1.0, 1.0]; Y(1) := Gsolve(f2, 0..1, Y(0)); Y(2) := Gsolve(f2, 1..2, Y(1)); Y(3) := Gsolve(f2, 2..3, Y(2)); Y(4) := Gsolve(f2, 3..4, Y(3)); Y(5) := Gsolve(f2, 4..5, Y(4))
Now, the invariant H_{2} is
preserved to working precision, whilst H_{1} is
only preserved to the tolerance specified by RelativeError
= tol
:
H2(Y(i))  H2(Y(0)) $ i = 1..5
H1(Y(i))  H1(Y(0)) $ i = 1..5
delete J1, J2, J3, Gsolve, f1, f2, Y, H1, H2:
We consider the "Toda lattice equations"
,
with a_{0} = a_{n} = 0. Introducing the tridiagonal n×n matrices
,
these equations can be encoded by the matrix ODE
.
The solution Y(t) is
known to be "isospectral", i.e., the eigenvalues of Y(t) do
not depend on the time parameter t.
As mentioned in the description of the option LieGroupAction
,
the solution of this type of matrix ODE is given by the group action Y(t)
= G(t) Y(0) G(t)^{
1} = G(t) Y(0) G(t)^{T},
where G(t) are
orthogonal matrices (note that f(Y) is
skew symmetric). The eigenvalues of the matrices Y(t) are
invariants of the group action.
The exact dynamics also preserves the tridiagonal form of the matrices. The numerical dynamics, however, fills in further elements. The following vector field f ignores alle elements outside the central bands:
f := proc(t, Y) local i, r; begin r := array(1..n, 1..n, [[0 $ n] $ n]); for i from 1 to n  1 do r[i + 1, i] := Y[i, i + 1]; r[i, i + 1] := Y[i, i + 1]; end_for; return(r) end_proc:
In the following, the initial value Y(0) is
specified by a matrix generated by the function matrix
. Consequently, both arguments G and Y are
passed to the Lie group action LAMBDA
as corresponding
matrices. They can be multiplied by the multiplication operator *
:
LAMBDA:= proc(G, Y) begin G*Y*(G::dom::transpose(G)) end_proc:
We define the initial value:
n := 3: Y(0) := matrix(n, n, [1, 1, 1], Banded)
Now, the dynamics is integrated from t = 0 to t = 1:
tol := 10^(4): Y(1) := numeric::odesolveGeometric(f, 0..1, Y(0), LieGroupAction = LAMBDA, RelativeError = tol)
The invariants of the dynamics are the eigenvalues of the matrices Y(t). They are preserved numerically:
numeric::eigenvalues(Y(0)) = numeric::eigenvalues(Y(1))
For comparison, we also solve the Toda lattice equations by
classical numerics using numeric::odesolve
.
The system is encoded by a vector Y =
[b_{1}, …, b_{n}, a_{1},
…, a_{n  1}] in
:
f := proc(t, Y) local a, b, i; begin b := [Y[i] $ i = 1..n]; a := [Y[n + i] $ i = 1..n1]; [2*a[1]^2, // = d/dt b[1] 2*(a[i1]^2  a[i]^2) $ i = 2..n1, // = d/dt b[i] 2*a[n1]^2, // = d/dt b[n] a[i]*(b[i]  b[i+1]) $ i = 1..n1 // = d/dt a[i] ] end_proc: solution := numeric::odesolve(f, 0..1, [1 $ 2*n  1], RelativeError = tol);
The invariants are only preserved up to the precision of the
solution determined by the tolerance set via RelativeError
= tol
:
Y(1) := array(1..n, 1..n, [[0 $ n] $ n]): for i from 1 to n do Y(1)[i, i] := solution[i]; end_for: for i from 1 to n1 do Y(1)[i, i + 1] := solution[n + i]; Y(1)[i + 1, i] := solution[n + i]; end_for: Y(1)
numeric::eigenvalues(Y(1))
Comparing these data with the previously computed eigenvalues
of the initial condition Y(0),
one sees that the invariants are not preserved numerically to the
working precision determined by DIGITS
.
delete f, LAMBDA, n, Y, tol, solution, i:

A procedure accepting two parameters 

A numerical real value for the initial time 

A numerical real value (the "time") 

The initial condition: a list,
a 1dimensional 

Option, specified as The procedure The default action is the usual matrix multiplication . With this option, the default group action As a group action, the procedure .
On the space of the n×n matrices and returns Y(t) = LAMBDA(G(t), Y_{0}). For the standard group action LAMBDA(G, Y) = G Y, this is the solution of the ODE . For homogeneous manifolds embedded in the n×n matrices,
the group action LAMBDA(G, Y)
= G Y G^{
1} may be considered. For this action,
the curve Y(t)
= LAMBDA(G(t), Y_{0}) returned
by
The procedure If the initial condition Y_{0} is
specified by a list or a 1dimensional Internally, the return value of It is recommended that  

Option, specified as Forces internal numerical RungeKutta steps to use step sizes
with local discretization errors below The internal control mechanism estimates the local relative
discretization error of a RungeKutta step and adjusts the step size
adaptively to keep this error below The default setting of ensures that the local discretization errors are of the same order of magnitude as numerical roundoff. Usually there is no need to use this option to change this setting.
However, occasionally the numerical evaluation of the RungeKutta
steps may be illconditioned or step sizes are so small that the time
parameter cannot be incremented by the step size within working precision.
In such a case, this option may be used to bound the local discretization
error by tol and
use a higher working precision given by Only real numerical values are accepted.
 

Option, specified as Switches off the internal error control and uses a RungeKutta
iteration with constant step size By default, A final step with smaller step size is used to match the end t of the integration interval , if is not an integer.
There is usually no need to invoke this option. However, occasionally the builtin adaptive error control mechanism may fail when integrating close to a singularity. In such a case, this option may be used to customize a control mechanism for the global error by using different step sizes and investigating the convergence of the corresponding results.  

Option, specified as With this option, The integer n controls
the size of the output list. For n =
1, all internal mesh points are returned. This
case may also be invoked by entering the simplified option The output list may be useful to inspect the internal numerical process. Also further graphical processing of the mesh data may be useful. 
The solution Y(t) is
returned as a list or as an array of floatingpoint values. The type
of the result matrix/vector coincides with the type of the input matrix/vector Y_{0}
.
With the option Alldata
, a list of mesh data
is returned.
[1] H. MuntheKaas and A. Zanna: "Numerical integration of differential equations on homogeneous manifolds", in F. Cucker (ed.), Foundations of Computational Mathematics, Springer (1997), pp. 305315.