Al–Mohy A H and Higham N J (2011) Computing the action of the matrix exponential, with an application to exponential integrators SIAM J. Sci. Statist. Comput. 33(2) 488-511

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

Parameters

Compulsory Input Parameters

1: $m$ – int64int32nag_int scalar: $m$ , the number of columns of the matrix $B$ .

Constraint: $m \geq 0$ .
2: $a (lda :)$ – double array: The first dimension of the array a must be at least $n$ .
The second dimension of the array a must be at least $n$ .

The $n$ by $n$ matrix $A$ .
3: $b (ldb :)$ – double array: The first dimension of the array b must be at least $n$ .
The second dimension of the array b must be at least $m$ .

The $n$ by $m$ matrix $B$ .
4: $t$ – double scalar: The scalar $t$ .

Optional Input Parameters

1: $n$ – int64int32nag_int scalar: Default: the first dimension of the arrays a, b. (An error is raised if these dimensions are not equal.)
$n$ , the order of the matrix $A$ .

Constraint: $n \geq 0$ .

Output Parameters

1: $a (lda :)$ – double array: The first dimension of the array a will be $n$ .
The second dimension of the array a will be $n$ .

$A$ is overwritten during the computation.
2: $b (ldb :)$ – double array: The first dimension of the array b will be $n$ .
The second dimension of the array b will be $m$ .

The $n$ by $m$ matrix $e^{t A} B$ .
3: $ifail$ – int64int32nag_int scalar: $ifail = 0$ unless the function detects an error (see Error Indicators and Warnings).

Error Indicators and Warnings

Errors or warnings detected by the function:

Cases prefixed with W are classified as warnings and do not generate an error of type NAG:error_n. See nag_issue_warnings.

W $ifail = 2$: $e^{t A} B$ has been computed using an IEEE double precision Taylor series, although the arithmetic precision is higher than IEEE double precision.

$ifail = - 1$: Constraint: $n \geq 0$ .

$ifail = - 2$: Constraint: $m \geq 0$ .

$ifail = - 4$: Constraint: $lda \geq n$ .

$ifail = - 6$: Constraint: $ldb \geq n$ .

$ifail = - 99$: An unexpected error has been triggered by this routine. Please contact NAG.

$ifail = - 399$: Your licence key may have expired or may not have been installed correctly.

$ifail = - 999$: Dynamic memory allocation failed.

Accuracy

For a symmetric matrix

A

(for which

A^{T} = A

) the computed matrix

e^{t A} B

is guaranteed to be close to the exact matrix, that is, the method is forward stable. No such guarantee can be given for non-symmetric matrices. See Section 4 of Al–Mohy and Higham (2011) for details and further discussion.

Further Comments

The matrix

e^{t A} B

could be computed by explicitly forming

e^{t A}

using nag_matop_real_gen_matrix_exp (f01ec) and multiplying

B

by the result. However, experiments show that it is usually both more accurate and quicker to use nag_matop_real_gen_matrix_actexp (f01ga).

The cost of the algorithm is

O (n^{2} m)

. The precise cost depends on

A

since a combination of balancing, shifting and scaling is used prior to the Taylor series evaluation.

Approximately

n^{2} + (2 m + 8) n

of real allocatable memory is required by nag_matop_real_gen_matrix_actexp (f01ga).

nag_matop_complex_gen_matrix_actexp (f01ha) can be used to compute

e^{t A} B

for complex

A

B

, and

t

. nag_matop_real_gen_matrix_actexp_rcomm (f01gb) provides an implementation of the algorithm with a reverse communication interface, which returns control to the user when matrix multiplications are required. This should be used if

A

is large and sparse.

Example

This example computes

e^{t A} B

, where

A = (\begin{array}{r} 0.7 & - 0.2 & 1.0 & 0.3 \\ 0.3 & 0.7 & 1.2 & 1.0 \\ 0.9 & 0.0 & 0.2 & 0.7 \\ 2.4 & 0.1 & 0.0 & 0.2 \end{array}),

B = (\begin{matrix} 0.1 & 1.2 \\ 1.3 & 0.2 \\ 0.0 & 1.0 \\ 0.4 & - 0.9 \end{matrix})

and

t = 1.2 .

Open in the MATLAB editor: f01ga_example

function f01ga_example


fprintf('f01ga example results\n\n');

a = [0.7,-0.2, 1.0, 0.3;
     0.3, 0.7, 1.2, 1.0;
     0.9, 0.0, 0.2, 0.7;
     2.4, 0.1, 0.0, 0.2];
b = [0.1, 1.2;
     1.3, 0.2;
     0.0, 1.0;
     0.4,-0.9];

t = 1.2;

% Compute exp(ta)b
[a, exptab, ifail] = f01ga(a, b, t);

disp('exp(tA)B');
disp(exptab);

f01ga example results

exp(tA)B
    0.2138    7.6756
    4.9980   11.6051
    0.8307    7.5468
    1.2406    9.7261

PDF version (NAG web site, 64-bit version, 64-bit version)

Chapter Contents

Chapter Introduction

NAG Toolbox