NAG Library Function Document

nag_real_gen_matrix_exp (f01ecc)

+− Contents

1 Purpose

2 Specification

3 Description

4 References

5 Arguments

6 Error Indicators and Warnings

7 Accuracy

8 Further Comments

+− 9 Example

9.1 Program Text

9.2 Program Data

9.3 Program Results

1 Purpose

nag_real_gen_matrix_exp (f01ecc) computes the matrix exponential,

e^{A}

, of a real

n

n

matrix

A

2 Specification

#include <nag.h>

#include <nagf01.h>

void	nag_real_gen_matrix_exp (Nag_OrderType order, Integer n, double a[], Integer pda, NagError *fail)

3 Description

e^{A}

is computed using a Padé approximant and the scaling and squaring method described in Higham (2005) and Higham (2008).

A

has a full set of eigenvectors

V

then

A

can be factorized as

A = V D V^{- 1},

where

D

is the diagonal matrix whose diagonal elements,

d_{i}

, are the eigenvalues of

A

e^{A}

is then given by

e^{A} = V e^{D} V^{- 1},

where

e^{D}

is the diagonal matrix whose

i

th diagonal element is

e^{d_{i}}

Note that

e^{A}

is not computed this way as to do so would, in general, be unstable.

4 References

Higham N J (2005) The scaling and squaring method for the matrix exponential revisited SIAM J. Matrix Anal. Appl. 26(4) 1179–1193

Higham N J (2008) Functions of Matrices: Theory and Computation SIAM, Philadelphia, PA, USA

Moler C B and Van Loan C F (2003) Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later SIAM Rev. 45 3–49

5 Arguments

1: order – Nag_OrderTypeInput

On entry: the order argument specifies the two-dimensional storage scheme being used, i.e., row-major ordering or column-major ordering. C language defined storage is specified by

order = Nag_RowMajor

. See Section 3.2.1.3 in the Essential Introduction for a more detailed explanation of the use of this argument.

Constraint:

order = Nag_RowMajor

or Nag_ColMajor.

2: n – IntegerInput

On entry:

n

, the order of the matrix

A

Constraint:

n \geq 0

3: a[ $pda \times n$ ] – doubleInput/Output

Note: the

(i, j)

th element of the matrix

A

is stored in

$a [(j - 1) \times pda + i - 1]$ when $order = Nag_ColMajor$ ;
$a [(i - 1) \times pda + j - 1]$ when $order = Nag_RowMajor$ .

On entry: the

n

n

matrix

A

On exit: the

n

n

matrix exponential

e^{A}

4: pda – IntegerInput

On entry: the stride separating row or column elements (depending on the value of order) in the array a.

Constraints:

if $order = Nag_ColMajor$ , $pda \geq \max (1, n)$ ;
if $order = Nag_RowMajor$ , $pda \geq n$ .

5: fail – NagError *Input/Output

The NAG error argument (see Section 3.6 in the Essential Introduction).

6 Error Indicators and Warnings

NE_ALLOC_FAIL: Dynamic memory allocation failed.
NE_BAD_PARAM: On entry, argument $〈value〉$ had an illegal value.
NE_INT: On entry, $n = 〈value〉$ .
Constraint: $n \geq 0$ .
NE_INT_2: On entry, $pda = 〈value〉$ and $n = 〈value〉$ .
Constraint: $pda \geq n$ .
NE_INTERNAL_ERROR: An internal error has occurred in this function. Check the function call and any array sizes. If the call is correct then please contact NAG for assistance.; An internal error occurred when processing the Padé approximant. Please contact NAG.
NE_SINGULAR: The linear equations to be solved are nearly singular and the Padé approximant probably has no correct figures; it is likely that this function has been called incorrectly.; The linear equations to be solved for the Padé approximant are singular; it is likely that this function has been called incorrectly.
NW_SOME_PRECISION_LOSS: The arithmetic precision is higher than that used for the Padé approximant computed matrix exponential.

7 Accuracy

For a normal matrix

A

(for which

A^{T} A = A A^{T}

) the computed matrix,

e^{A}

, is guaranteed to be close to the exact matrix, that is, the method is forward stable. No such guarantee can be given for non-normal matrices. See Section 10.3 of Higham (2008) for details and further discussion.

For discussion of the condition of the matrix exponential see Section 10.2 of Higham (2008).

8 Further Comments

The cost of the algorithm is

O (n^{3})

; see Algorithm 10.20 of Higham (2008).

As well as the excellent book cited above, the classic reference for the computation of the matrix exponential is Moler and Van Loan (2003).

9 Example

This example finds the matrix exponential of the matrix

A = (\begin{array}{r} 1 & 2 & 2 & 2 \\ 3 & 1 & 1 & 2 \\ 3 & 2 & 1 & 2 \\ 3 & 3 & 3 & 1 \end{array}) .

NAG Library Function Documentnag_real_gen_matrix_exp (f01ecc)