expmv

Matrix exponential multiplied by vector

Since R2023b

Syntax

F = expmv(A,b,t)

F = expmv(A,b)

F = expmv(A,b,tvals)

Description

F = expmv(A,b,t) calculates the product of a matrix exponential and a vector. This function computes expm(t*A)*b without explicitly forming the matrix exponential expm(t*A). Here, A is an n-by-n square matrix, b is an n-by-1 column vector, and t is a scalar. The output F is an n-by-1 column vector.

A can also be a function handle, where A(flag,X) must accept an input matrix X, create a square matrix Am with an order that is equal to the column length of X, and return the following values for the indicated flag:

• A("real",X) returns logical 1 (true) if the square matrix Am is real and returns logical 0 (false) otherwise.

• A("notransp",X) returns Am*X.

• A("transp",X) returns Am'*X.

F = expmv(A,b) computes expm(A)*b without explicitly forming the matrix exponential expm(A). This syntax is the same as expmv(A,b,1).

F = expmv(A,b,tvals) computes expm(t*A)*b without explicitly forming expm(t*A) for each element t in tvals. Here, tvals is a vector of length m and the output F is an n-by-m matrix such that each column of F is F(:,j) = expm(tvals(j)*A)*b for j = 1:m.

For an equally spaced tvals, such as a time vector with a fixed time step, expmv uses an efficient algorithm that reuses information based on the spacing pattern of tvals. MATLAB^® uses the isuniform function to check whether tvals is equally spaced.

Examples

Matrix Exponential Times Vector

Find the matrix exponential times a vector: ${\mathit{e}}^{\mathit{t}\mathit{A}}\overrightarrow{\mathit{b}}$.

First, create the scalar t, the square matrix A, and the column vector b.

t = 0.25;
A = [1 1i 0; 0 0 2i; 0 0 -1];
b = [2; 1; 3];

Compute the matrix exponential times a vector without explicitly forming ${\mathit{e}}^{\mathit{t}\mathit{A}}$.

F = expmv(A,b,t)

F = 3×1 complex

   2.3796 + 0.2840i
   1.0000 + 1.3272i
   2.3364 + 0.0000i

To form ${\mathit{e}}^{\mathit{t}\mathit{A}}$ explicitly, use expm.

C = expm(t*A)

C = 3×3 complex

   1.2840 + 0.0000i   0.0000 + 0.2840i  -0.0628 + 0.0000i
   0.0000 + 0.0000i   1.0000 + 0.0000i   0.0000 + 0.4424i
   0.0000 + 0.0000i   0.0000 + 0.0000i   0.7788 + 0.0000i

Show that the result of expm(t*A)*b is equal to expmv(A,b,t) within round-off error.

F2 = C*b;
norm(F - F2)

ans = 
8.8991e-16

Matrix Exponential of Magic Square Times Vector

Find the matrix exponential times a vector: ${\mathit{e}}^{\mathit{t}\mathit{A}}\overrightarrow{\mathit{b}}$, where $\mathit{A}$ is a magic square with an order that is determined by the length of vector $\overrightarrow{\mathit{b}}$.

Create a function in a file named Afunction.m. This function accepts an input matrix X and creates a magic square Am with an order that is equal to the column length of X. This function also accepts an input flag and returns the following values for the indicated flag:

• Afunction("real",X) returns logical 1 (true) if the square matrix Am is real and returns logical 0 (false) otherwise.

• Afunction("notransp",X) returns Am*X.

• Afunction("transp",X) returns Am'*X.

function Y = Afunction(flag,X)
n = size(X,1);
Am = magic(n);
switch flag
    case "real"
        Y = isreal(Am);
    case "notransp"
        Y = Am*X;
    case "transp"
        Y = Am'*X;
end
end

Create a function handle to Afunction.

A = @(flag,X) Afunction(flag,X);

Use the function handle as an input to expmv to compute ${\mathit{e}}^{\mathit{t}\mathit{A}}\overrightarrow{\mathit{b}}$.

t = 0.03;
b = [1; 2; 4; 5];
F = expmv(A,b,t)

F = 4×1

    6.1135
    7.4895
    9.4533
   10.2221

Matrix Exponential of Sparse Matrix Times Vector

Find the matrix exponential times a vector: ${\mathit{e}}^{\mathit{A}}\overrightarrow{\mathit{b}}$, where $\mathit{A}$ is a sparse matrix.

Create a tridiagonal sparse square matrix A and a column vector b.

A = gallery("tridiag",10000,1,-2,1);
b = (1:size(A,1))';

First, compute the matrix exponential of A explicitly and multiply it by b. Use a pair of tic and toc calls to measure the time required for this computation.

tic
F = expm(A)*b;
toc

Elapsed time is 18.533127 seconds.

For comparison, use expmv to compute the matrix exponential of A times b. This computation is significantly faster because expmv does not compute the matrix exponential explicitly.

tic
F = expmv(A,b);
toc

Elapsed time is 0.040429 seconds.

Solution of System of Ordinary Differential Equations

Given a system of ordinary differential equations (ODEs) that has the form

$$\frac{\mathrm{d}}{\mathrm{dt}}\overrightarrow{\mathit{X}}\left(\mathit{t}\right)=\mathit{A}\overrightarrow{\mathit{X}}\left(\mathit{t}\right)$$,

the solution is a matrix exponential times a vector

$\overrightarrow{\mathit{X}}\left(\mathit{t}\right)={\mathit{e}}^{\mathit{t}\mathit{A}}\overrightarrow{{\mathit{X}}_0}$.

To compute $\overrightarrow{\mathit{X}}\left(\mathit{t}\right)$ efficiently for a range of time values $\mathit{t}$ that are equally spaced with a fixed time step, you can use expmv.

For example, find the solution of the following system of ODEs for 100 evenly spaced time values in the interval [0,10].

X0 = [5; -1; 1];
A = [0 5 0; 0 -5 2; 0 2 -1];
tvals = linspace(0,10,100);
F = expmv(A,X0,tvals);

Plot the solution as a function of time.

plot(tvals,F)
legend("x1","x2","x3")

Input Arguments

A — Input matrix

square matrix | function handle

Input matrix, specified as a square matrix or a function handle.

Data Types: single | double | function_handle
Complex Number Support: Yes

b — Input vector

column vector

Input vector, specified as a column vector. If the size of A is n-by-n, then the size of b must be n-by-1.

Data Types: single | double
Complex Number Support: Yes

t — Exponential scaling factor

real scalar

Exponential scaling factor, specified as a real scalar.

Data Types: single | double

tvals — Exponential scaling factor or time vector

real vector

Exponential scaling factor or time vector, specified as a real vector.

Data Types: single | double

References

[1] Al-Mohy, Awad H., and Nicholas J. Higham. “Computing the Action of the Matrix Exponential, with an Application to Exponential Integrators.” SIAM Journal on Scientific Computing 33, no. 2 (January 2011): 488–511. https://doi.org/10.1137/100788860.

Version History

Introduced in R2023b