pyyeti.expmint.expmint

pyyeti.expmint.expmint(A, h, geti2=False)

Compute the matrix exponential and its integral(s) using Pade approximation.

Parameters:

• A ((M, M) array_like or sparse matrix) – 2D Array or Matrix (sparse or dense) to be exponentiated

• h (scalar) – Time step

• geti2 (bool) – If True, also return the I2 integral below. Useful for first order holds.

Returns:

• E ((M, M) ndarray) – Matrix exponential of A*h: exp(A*h)

• I ((M, M) ndarray) – Integral of exp(A*t) dt from 0 to h

• I2 ((M, M) ndarray; optional) – Integral of exp(A*t)*t dt from 0 to h. Only returned if geti2 is True.

Notes

This routine is modeled after and augments scipy.linalg.expm(). The power series expansions for these matrices are (I = identity):

E = I + A*h + (A*h)**2/2! + (A*h)**3/3! + ...
I1 = h*(I + A*h/2 + (A*h)**2/3! + (A*h)**3/4! + ...)
I2 = h*h*(I/2 + A*h/3 + (A*h)**2/(4*2!) + (A*h)**3/(5*3!) + ...)

If A is non-singular, the exact solutions for I1 and I2 are:

E = exp(A*h)
I1 = inv(A)*(E-I)
I2 = inv(A)*(E*h-I1)

The Pade approximants for those power series are used for E and I1. If necessary, the ‘squaring and scaling’ method will be used such that a 13th order Pade approximation will be accurate. See references [1], [2], and [3] for more information.

For I2, a Pade approximation is used if a 3rd, 5th, 7th or 9th order is accurate. If not, A is checked for singularity. If non-singular, the exact solution show above is used. If A is singular, a the power series is used directly until it converges (and a warning message is printed about using a finer time step.)

References

[1]

Awad H. Al-Mohy and Nicholas J. Higham (2009) “A New Scaling and Squaring Algorithm for the Matrix Exponential.” SIAM Journal on Matrix Analysis and Applications. 31 (3). pp. 970-989. ISSN 1095-7162

[2]

Nicholas J. Higham (2005) “The Scaling and Squaring Method for the Matrix Exponential Revisited.” SIAM Journal on Matrix Analysis and Applications. Vol 26, No. 4, pp 1179-1193.

[3]

David Westreich (1990) “A Practical Method for Computing the Exponential of a Matrix and its Integral.” Communications in Applied Numerical Methods, Vol 6, 375-380.

Examples

>>> from pyyeti import expmint
>>> import numpy as np
>>> import scipy.linalg as la
>>> np.set_printoptions(4)
>>> a = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
>>> e, i, i2 = expmint.expmint(a, .05, True)
>>> e
array([[ 1.0996,  0.1599,  0.2202],
       [ 0.3099,  1.3849,  0.46  ],
       [ 0.5202,  0.61  ,  1.6998]])
>>> i
array([[ 0.052 ,  0.0034,  0.0048],
       [ 0.0067,  0.0583,  0.01  ],
       [ 0.0114,  0.0133,  0.0651]])
>>> i2
array([[ 0.0013,  0.0001,  0.0002],
       [ 0.0002,  0.0015,  0.0003],
       [ 0.0004,  0.0005,  0.0018]])