Main Content

expmv

Matrix exponential multiplied by vector

Since R2023b

    Description

    example

    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.

    example

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

    example

    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

    collapse all

    Find the matrix exponential times a vector: etAb.

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

    F = expmv(A,b,t)
    F = 3×1 complex
    
       2.3796 + 0.2840i
       1.0000 + 1.3272i
       2.3364 + 0.0000i
    
    

    To form etA 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
    

    Find the matrix exponential times a vector: etAb, where A is a magic square with an order that is determined by the length of vector 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 etAb.

    t = 0.03;
    b = [1; 2; 4; 5];
    F = expmv(A,b,t)
    F = 4×1
    
        6.1135
        7.4895
        9.4533
       10.2221
    
    

    Find the matrix exponential times a vector: eAb, where 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.
    

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

    ddtX(t)=AX(t),

    the solution is a matrix exponential times a vector

    X(t)=etAX0.

    To compute X(t) efficiently for a range of time values 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")

    Figure contains an axes object. The axes object contains 3 objects of type line. These objects represent x1, x2, x3.

    Input Arguments

    collapse all

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

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

    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

    Exponential scaling factor, specified as a real scalar.

    Data Types: single | double

    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

    See Also

    | | | | |