Numerical solution of an ordinary differential equation on a homogeneous manifold

MuPAD® notebooks will be removed in a future release. 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, t0 .. t, Y0, <LieGroupAction = LAMBDA>, <method>, <RelativeError = tol>, <Stepsize = h>, <Alldata = n>)
numeric::odesolveGeometric(t0 .. t, f, Y0, <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 Y0 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) Y0, the matrix ODE

is solved on the space of the complex n×n matrices (1n, n is the identity matrix). Following Munthe-Kaas [1], the ansatz reduces a time step for the ode above to the solution of the matrix ODE


where and [u, f] = uf - fu 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 Runge-Kutta 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 = GY0 lies on the orbit of the Lie group through the initial value Y0. Thus, the geometrical invariants of the homogeneous manifold are preserved in the course of the numerical integration.

The input data t0 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 Y0 defines the space in which the homogeneous manifold containing the solution is embedded.

If Y0 is a list with n entries or a 1-dimensional array array(1..n), then the solution Y(t) consists of vectors from a submanifold of or , respectively.

If Y0 is specified as a 2-dimensional 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, 2-dimensional n×m arrays are used to represent the points on the manifold where m = 1 for vectors in or . It is recommended to specify Y0 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 floating-point values t and matrices/vectors Y of the same domain type as the initial condition Y0.

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.


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

Environment Interactions

The function is sensitive to the environment variable DIGITS, which determines the numerical working precision.


Example 1

We consider the initial value problem

for with fixed parameters , J2 = 1, J3 = 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 (f1 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 “J-orthogonal” group SO(J, 3) through the initial point. This group consists of matrices G satisfying GTJG = J, where J = diag(J1, J2, J3). 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    ]])
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 H1 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 H2 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    ]])

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 H2 is preserved to working precision, whilst H1 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:

Example 2

We consider the “Toda lattice equations”


with a0 = an = 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;
  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];

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)

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 = [b1, …, bn, a1, …, an - 1] in :

f := proc(t, Y) 
local a, b, i;
  b := [Y[i] $ i = 1..n];
  a := [Y[n + i] $ i = 1..n-1];
  [-2*a[1]^2,                           // = d/dt b[1]
    2*(a[i-1]^2 - a[i]^2) $ i = 2..n-1, // = d/dt b[i]
    2*a[n-1]^2,                         // = d/dt b[n]
    a[i]*(b[i] - b[i+1]) $ i = 1..n-1   // = d/dt a[i]

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];
for i from 1 to n-1 do
  Y(1)[i, i + 1] := solution[n + i];
  Y(1)[i + 1, i] := solution[n + i];



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 (t, Y)


A numerical real value for the initial time


A numerical real value (the “time”)


The initial condition: a list, a 1-dimensional array(1..n), a 2-dimensional array(1..n, 1..m), or an n×m matrix of category Cat::Matrix with numerical entries



Option, specified as LieGroupAction = LAMBDA

The procedure LAMBDA = proc(G, Y) ... end_proc defines the action of the group element G (an n×n matrix) on the point Y on the homogeneous manifold (an n×m matrix or an n dimensional vector). This procedure must return a corresponding point (a matrix or a vector).

The default action is the usual matrix multiplication .

With this option, the default group action LAMBDA: of the n×n matrices G acting on the n×m matrices or n dimensional vectors Y by left multiplication may be replaced by other group actions.

As a group action, the procedure LAMBDA must satisfy LAMBDA(1n, n, Y) = Y and


numeric::odesolveGeometric computes the solution of the matrix ODE

On the space of the n×n matrices and returns Y(t) = LAMBDA(G(t), Y0).

For the standard group action LAMBDA(G, Y) = GY, this is the solution of the ODE .

For homogeneous manifolds embedded in the n×n matrices, the group action LAMBDA(G, Y) = GYG- 1 may be considered. For this action, the curve Y(t) = LAMBDA(G(t), Y0) returned by numeric::odesolveGeometric is the solution of the ODE . Cf. Example 2.

LAMBDA(G, Y) is called with n×m matrices or n dimensional vectors Y of the same domain type as the initial condition Y0. If Y0 is a matrix generated by the function matrix, then also the n×n matrix G is passed to LAMBDA as such a matrix object. In all other cases, G is passed as a 2-dimensional array(1..n, 1..n).

The procedure LAMBDA should return a 2-dimensional array(1..n, 1..m) or a corresponding matrix of category Cat::Matrix.

If the initial condition Y0 is specified by a list or a 1-dimensional array(1..n), LAMBDA may also return a corresponding list or array.

Internally, the return value of LAMBDA is converted to a 2-dimensional array(1..n, 1..m) where m = 1 if a list or a 1-dimensional array is returned.

It is recommended that LAMBDA returns a 2-dimensional array(1..n, 1..m) in order to avoid the overhead of internal conversions.


Option, specified as RelativeError = tol

Forces internal numerical Runge-Kutta steps to use step sizes with local discretization errors below tol. This tolerance must be a numerical real value . The default tolerance is .

The internal control mechanism estimates the local relative discretization error of a Runge-Kutta step and adjusts the step size adaptively to keep this error below tol.

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 Runge-Kutta steps may be ill-conditioned 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 DIGITS.

Only real numerical values are accepted.


Usually, the global error of the numeric approximation returned by numeric::odesolveGeometric is larger than the local errors bounded by tol. Although the result is displayed with DIGITS decimal places, one should not expect that all of them are correct. The relative precision of the final result is tol at best!


Option, specified as Stepsize = h

Switches off the internal error control and uses a Runge-Kutta iteration with constant step size h. h must be a numerical real value.

By default, numeric::odesolveGeometric uses an adaptive step size control mechanism in the numerical iteration. The option Stepsize = h switches off this adaptive mechanism and uses the Runge-Kutta method specified by method (or the default method DOPRI78) with constant step size h.

A final step with smaller step size is used to match the end t of the integration interval , if is not an integer.


With this option, there is no automatic error control! Depending on the problem and on the order of the method, the result may be a poor numerical approximation of the exact solution.

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 Alldata = n

With this option, numeric::odesolveGeometric returns a list of numerical mesh points [[t0, Y0], [t1, Y1], …, [t, Y(t)]] generated by the internal Runge-Kutta iteration.

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 Alldata, which is equivalent to Alldata = 1. For n > 1, only each n-th mesh point is stored in the list. The list always contains the initial point [t0, Y0] and the final point [t, Y(t)]. For n ≤ 0, only the data [[t0, Y0], [t, Y(t)]] are returned.

The output list may be useful to inspect the internal numerical process. Also further graphical processing of the mesh data may be useful.

Return Values

The solution Y(t) is returned as a list or as an array of floating-point values. The type of the result matrix/vector coincides with the type of the input matrix/vector Y0.

With the option Alldata, a list of mesh data is returned.


[1] H. Munthe-Kaas and A. Zanna: “Numerical integration of differential equations on homogeneous manifolds”, in F. Cucker (ed.), Foundations of Computational Mathematics, Springer (1997), pp. 305-315.

See Also

MuPAD Functions

MuPAD Graphical Primitives