Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

Chapter Contents
Chapter Introduction
NAG Toolbox

# NAG Toolbox: nag_matop_complex_gen_matrix_frcht_exp (f01kh)

## Purpose

nag_matop_complex_gen_matrix_frcht_exp (f01kh) computes the Fréchet derivative $L\left(A,E\right)$ of the matrix exponential of a complex $n$ by $n$ matrix $A$ applied to the complex $n$ by $n$ matrix $E$. The matrix exponential ${e}^{A}$ is also returned.

## Syntax

[a, e, ifail] = f01kh(a, e, 'n', n)
[a, e, ifail] = nag_matop_complex_gen_matrix_frcht_exp(a, e, 'n', n)

## Description

The Fréchet derivative of the matrix exponential of $A$ is the unique linear mapping $E⟼L\left(A,E\right)$ such that for any matrix $E$
 $eA+E - e A - LA,E = oE .$
The derivative describes the first-order effect of perturbations in $A$ on the exponential ${e}^{A}$.
nag_matop_complex_gen_matrix_frcht_exp (f01kh) uses the algorithms of Al–Mohy and Higham (2009a) and Al–Mohy and Higham (2009b) to compute ${e}^{A}$ and $L\left(A,E\right)$. The matrix exponential ${e}^{A}$ is computed using a Padé approximant and the scaling and squaring method. The Padé approximant is then differentiated in order to obtain the Fréchet derivative $L\left(A,E\right)$.

## References

Al–Mohy A H and Higham N J (2009a) A new scaling and squaring algorithm for the matrix exponential SIAM J. Matrix Anal. 31(3) 970–989
Al–Mohy A H and Higham N J (2009b) Computing the Fréchet derivative of the matrix exponential, with an application to condition number estimation SIAM J. Matrix Anal. Appl. 30(4) 1639–1657
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

## Parameters

### Compulsory Input Parameters

1:     $\mathrm{a}\left(\mathit{lda},:\right)$ – complex array
The first dimension of the array a must be at least ${\mathbf{n}}$.
The second dimension of the array a must be at least ${\mathbf{n}}$.
The $n$ by $n$ matrix $A$.
2:     $\mathrm{e}\left(\mathit{lde},:\right)$ – complex array
The first dimension of the array e must be at least ${\mathbf{n}}$.
The second dimension of the array e must be at least ${\mathbf{n}}$.
The $n$ by $n$ matrix $E$

### Optional Input Parameters

1:     $\mathrm{n}$int64int32nag_int scalar
Default: the first dimension of the arrays a, e. (An error is raised if these dimensions are not equal.)
$n$, the order of the matrix $A$.
Constraint: ${\mathbf{n}}\ge 0$.

### Output Parameters

1:     $\mathrm{a}\left(\mathit{lda},:\right)$ – complex array
The first dimension of the array a will be ${\mathbf{n}}$.
The second dimension of the array a will be ${\mathbf{n}}$.
The $n$ by $n$ matrix exponential ${e}^{A}$.
2:     $\mathrm{e}\left(\mathit{lde},:\right)$ – complex array
The first dimension of the array e will be ${\mathbf{n}}$.
The second dimension of the array e will be ${\mathbf{n}}$.
The Fréchet derivative $L\left(A,E\right)$
3:     $\mathrm{ifail}$int64int32nag_int scalar
${\mathbf{ifail}}={\mathbf{0}}$ unless the function detects an error (see Error Indicators and Warnings).

## Error Indicators and Warnings

Errors or warnings detected by the function:
${\mathbf{ifail}}=1$
The linear equations to be solved for the Padé approximant are singular; it is likely that this function has been called incorrectly.
${\mathbf{ifail}}=2$
${e}^{A}$ has been computed using an IEEE double precision Padé approximant, although the arithmetic precision is higher than IEEE double precision.
${\mathbf{ifail}}=3$
${\mathbf{ifail}}=-1$
Constraint: ${\mathbf{n}}\ge 0$.
${\mathbf{ifail}}=-3$
Constraint: $\mathit{lda}\ge {\mathbf{n}}$.
${\mathbf{ifail}}=-5$
Constraint: $\mathit{lde}\ge {\mathbf{n}}$.
${\mathbf{ifail}}=-99$
${\mathbf{ifail}}=-399$
Your licence key may have expired or may not have been installed correctly.
${\mathbf{ifail}}=-999$
Dynamic memory allocation failed.

## Accuracy

For a normal matrix $A$ (for which ${A}^{\mathrm{H}}A=A{A}^{\mathrm{H}}$) 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), Al–Mohy and Higham (2009a) and Al–Mohy and Higham (2009b) for details and further discussion.

The cost of the algorithm is $O\left({n}^{3}\right)$ and the complex allocatable memory required is approximately $9{n}^{2}$; see Al–Mohy and Higham (2009a) and Al–Mohy and Higham (2009b).
If the matrix exponential alone is required, without the Fréchet derivative, then nag_matop_complex_gen_matrix_exp (f01fc) should be used.
If the condition number of the matrix exponential is required then nag_matop_complex_gen_matrix_cond_exp (f01kg) should be used.
As well as the excellent book Higham (2008), the classic reference for the computation of the matrix exponential is Moler and Van Loan (2003).

## Example

This example finds the matrix exponential ${e}^{A}$ and the Fréchet derivative $L\left(A,E\right)$, where
 $A = 1+0i 2+0i 2+0i 2+i 3+2i 1i+0 1i+0 2+i 3+2i 2+0i 1i+0 2+i 3+2i 3+2i 3+2i 1+i and E = 1i+0 2+0i 2i+0 4+i 3+2i 0i+0 1i+0 0+i 0+2i 0+0i 1i+0 0i+ 1+0i 2+2i 0+3i 1i+ .$
```function f01kh_example

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

% Exponential of matrix A and Frechet derivative of exp(A)E.

a =  [1+ i, 2+ i, 2+ i, 2+i;
3+2i, 1,    1,    2+i;
3+2i, 2+ i, 1,    2+i;
3+2i, 3+2i, 3+2i, 1+i];

e =  [1,    2+ i, 2,    4+i;
3+2i, 0,    1,    0+i;
0+2i, 0+ i, 1,    0;
1+ i, 2+2i, 0+3i, 1];

[expa, lae, ifail] = f01kh(a,e);

[ifail] = x04da('General', ' ', expa, 'exp(A):');
disp(' ');
[ifail] = x04da('General', ' ', lae, 'L_exp(A,E):');

```
```f01kh example results

exp(A):
1          2          3          4
1   -157.9003  -194.6526  -186.5627  -155.7669
-754.3717  -555.0507  -475.4533  -520.1876

2   -206.8899  -225.4985  -212.4414  -186.5627
-694.7443  -505.3938  -431.0611  -475.4533

3   -208.7476  -238.4962  -225.4985  -194.6526
-808.2090  -590.8045  -505.3938  -555.0507

4   -133.3958  -208.7476  -206.8899  -157.9003
-1085.5496  -808.2090  -694.7443  -754.3717

L_exp(A,E):
1          2          3          4
1   1571.5852   778.4238   500.2085   740.7485
-4640.2429 -3719.8308 -3246.0234 -3424.1963

2   1472.7846   731.6608   473.2569   692.0895
-4273.5048 -3432.5961 -2990.9285 -3148.4635

3   1996.4848  1107.9174   782.1266  1031.5808
-4568.8881 -3714.9923 -3249.1926 -3400.8557

4   3327.1347  2015.2763  1514.3130  1873.9421
-5829.0773 -4810.2591 -4234.6812 -4404.0163
```