# Documentation

### This is machine translation

Translated by
Mouseover text to see original. Click the button below to return to the English verison of the page.

# `numeric`::`odesolveGeometric`

Numerical solution of an ordinary differential equation on a homogeneous manifold

MATLAB live scripts support most MuPAD functionality, though there are some differences. For more information, see Convert MuPAD Notebooks to MATLAB Live Scripts.

## Syntax

```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`>)
```

## Description

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

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

## Examples

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

```f := proc(t, Y) local a, b, i; begin 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] ] 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 n-1 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:`

## Parameters

 `f` A procedure accepting two parameters ```(t, Y)``` `t0` A numerical real value for the initial time `t` A numerical real value (the “time”) `Y0` 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

## Options

`LieGroupAction`

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.

`RelativeError`

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.

### Note

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!

`Stepsize`

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.

### Note

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.

`Alldata`

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.

## References

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