Abstract
In this paper, two new versions of the Schur algorithm for computing the matrix exponential of an complex matrix are presented. Instead of the Schur form, these algorithms use the Jordan–Schur form of . The Jordan–Schur form is found by less computation and it is determined more reliable than the reduction to Jordan form since it is obtained using only unitary similarity transformations. In contrast to the known methods, the diagonal blocks of the matrix exponential are obtained by using finite Taylor series. This improves the accuracy and avoids the decisions made about the termination of the series expansion. The off-diagonal blocks of the exponential are determined by modifications of the Schur–Parlett or Schur–Fréchet method, which takes advantage of the Jordan–Schur form of the matrix. The numerical features of the new algorithms are discussed, revealing their advantages and disadvantages in comparison with the other methods for computing the matrix exponential. Computational experiments show that using the new algorithms, the matrix exponential is determined in certain cases with higher accuracy than some widely used methods, however, at the price of an increase in the computational cost which is of order . It is shown that the Jordan–Schur algorithms for computing the matrix exponential are appropriate for matrices with multiple eigenvalues and are especially efficient in cases of large Weyr characteristics.
1. Introduction
The computation of the matrix exponential of a real or complex matrix is an important task in the solution of linear differential equations, obtaining the time response and discretization of control systems, computing the integral of the exponential, and several other problems arising in applied mathematics. As it is well known [Moler and Van Loan [1, 2]], there exist at least nineteen “dubious” methods for computing the matrix exponential and there is no single “best” method for solving this important numerical problem. Among the good methods for computing the exponential that works in many cases, is the scaling and squaring method proposed in [Ward [3]], [Van Loan [4]] and improved in [Higham [5]]. An efficient method in this group combines the computation of the matrix exponential by using Padé approximations with scaling and squaring to improve the accuracy; see, for instance, [Higham [5]], [Al-Mohy and Higham [6]], and [Dieci and Papini [7]]. Alternative methods of the same group are based on the Taylor series, see [Caliari and Zivcovich [8]] and [Sastre et al. [9]].
The methods to compute the matrix exponential of non-normal matrices that are based on the Schur form of the matrix [Golub and Van Loan [10]; Ch. 9] and [Higham [11]; Ch. 10] are frequently preferred due to the numerical stability of the QR algorithm for computing the Schur form. Among the efficient methods for this aim is the Schur–Parlett algorithm [Davis and Higham [12]], which has comparable performance with the scaling and squaring method. The numerical properties of this algorithm are investigated in [Higham [13]; Ch. 9], see also Algorithm 10.23 in [Higham [11]; Ch. 10]. Another opportunity to compute the off-diagonal blocks of the matrix exponential by using the Schur form is provided by the Schur–Fréchet algorithm developed and analyzed in [Najfeld and Havel [14]], [Kenney and Laub [15]], and [Higham [13]; Ch. 9].
The difficulty with the Schur–Parlett and Schur–Fréchet methods comes from the necessity to evaluate the exponential of a triangular matrix which may lead to large errors in some cases. Due to the occasional failures of these methods, there is a need for better versions that are more robust from the numerical point of view. In most cases, the failures occur in the presence of multiple eigenvalues of associated with complicated Jordan structures, so it is natural to make an attempt to use the Jordan form of the original matrix. However, determining the Jordan form of a matrix requires the use of nonunitary transformations which may cause numerical instability. More importantly, the computing of the Jordan form is an ill-posed computational problem since infinitely small perturbations of may lead to abrupt changes in the Jordan structure, see [Golub and Wilkinson [16]] and [Demmel [17]]. That is why it seems advantageous to implement some other form of which carries information about the Jordan structure of but can be obtained by well-conditioned similarity transformations and is relatively insensitive to perturbations of the original matrix. Such form is the Jordan–Schur form of a square matrix, described in [Bai et al. [18]; Sect. 2.5], [Kågström and Wiberg [19]], and [Elmroth et al. [20]]. The Jordan–Schur can be considered as an intermediate form between Schur and Jordan form, and it is possible to compute it by a minor modification of a part of the algorithm of Kågström and Ruhe [Kågström and Ruhe [21, 22]] intended for finding the numerical Jordan form of a matrix. It should be noted that the reduction to the Jordan–Schur form is done by using only unitary similarity transformations and thus it is more reliable numerically than the reduction to the Jordan form. The use of the Jordan–Schur form allows for deriving more accurate algorithms for computing the matrix exponential, however at the price of an increase in the volume of computational work.
In this paper, two new versions of the Schur method for computing the matrix exponential are presented. Instead of the Schur form, these algorithms use the Jordan–Schur form of a given complex matrix, which leads to some advantages in computing the exponential. In particular, the exponentials of the diagonal blocks are determined by using the finite Taylor series which improves the accuracy, and the off-diagonal blocks are found by modification of the Schur–Parlett or Schur–Fréchet method that takes advantage of the Jordan–Schur form of the matrix. It is confirmed by numerical examples that the new method may produce more accurate results than some widely used implementations of the scaling and squaring and Schur method.
The paper is organized as follows. The Weyr and Segre characteristics, along with the Jordan–Schur form are introduced briefly in Section 2. The new algorithms for computing the matrix exponential, based on the Jordan–Schur form are described in Section 3 and their numerical properties are discussed in Section 4. Four numerical experiments illustrating the behavior of the proposed algorithms are presented in Section 5, and some conclusions are drawn in Section 6.
A notation will be the following. is the set of complex numbers; is the transposed of ; is the Hermitian transpose (the complex conjugate transpose) of ; is the unit matrix; is the th block of the th power of the block matrix ; is an -by- block matrix; is the th eigenvalue of ; , and are the spectral norm, the Frobenius norm, and the infinity norm of , respectively; is the dimension of the subspace ; is the null space of ; is the rank of ; is a perturbation of ; is the condition number of with respect to the inversion in the Frobenius norm; is the exponential of ; is the condition number of the matrix exponential; is the Kronecker product of and ; is the separation between square matrices and ; is the condition number of the Sylvester equation.
2. The Jordan–Schur Form
The Jordan–Schur form is closely related to Jordan and Weyr canonical forms [Horn and Johnson [23], sect. 3.4]. To present this form, one needs to introduce the Segre and Weyr characteristics associated with multiple eigenvalues.
2.1. Segre and Weyr Characteristics
Consider the general case when a matrix has distinct eigenvalues with algebraic multiplicities.respectively.
Denote by,andThe Segre and Weyr characteristics of [Horn and Johnson [23]; p. 170], [Shapiro [24]], respectively, are associated with , where and . The elements of the Segre characteristic are equal to the sizes of the Jordan blocks associated with in the Jordan form of . The quantity is equal to the number of the linearly independent eigenvectors of corresponding to , i.e., to the number of the Jordan blocks associated with . In turn, the elements of the Weyr characteristic satisfy the relationships.
Note that,is the index of nilpotency of , i.e., the smallest positive integer for which.
Let be a partition of the positive integer (i.e., ) and let denotes the number of the integers , which are greater or equal to . Then, the integers form a partition of , known as a conjugate partition of . The relationship between the Segre and Weyr characteristics of associated with an eigenvalue with multiplicity , is revealed by the fact that these characteristics are conjugated partitions of [Horn and Johnson [23]; p. 172]. Thus, the Weyr characteristic determines uniquely the Segre characteristic of associated with and vice versa. Note that, the largest element (the size of the largest Jordan block) of the Segre characteristic associated with is equal to the number of the elements of the corresponding Weyr characteristic. In turn, the first element of the Weyr characteristic (the geometric multiplicity of the eigenvalue ) is equal to the number of the elements of the Segre characteristic (the number of Jordan blocks) associated with .
2.2. Reduction to Jordan–Schur Form
Using unitary similarity transformations with a matrix , the matrix is reduced to the upper block triangular form (the Jordan–Schur form).where the diagonal blocksare in the so-called staircase form and the off-diagonal blocks have no special form. Thus, the structure of the Jordan–Schur form is determined by the Weyr characteristics associated with . These characteristics can be used to find easily the corresponding Segre characteristics of .
A possible way to determine the Jordan–Schur form of a matrix is to implement a modification of the algorithm of Kågström and Ruhe [Kågström and Ruhe [21, 22]] intended to find the Jordan form of a matrix. This algorithm is based on the idea that using a sequence of unitary transformations it is possible to find the dimensions of the null spaces (and hence the elements of the Weyr characteristic) without computing the matrix powers , and consequently, the size of the Jordan blocks associated with multiple eigenvalues. As a result, the matrix is reduced to a staircase form and the corresponding algorithm is called a staircase algorithm because of the shape of the resulting matrix (for a detailed history of this algorithm, see [Edelman and Ma [25]]).
The reduction of an -th order complex matrix to Jordan–Schur form can be carried out by a minor modification of the algorithm of Kågström and Ruhe [Petkov [26]].
Concerning the computational cost, the reduction to the Jordan–Schur form requires more operations than the reduction to the Schur form. This is mainly due to the reduction of the diagonal blocks to the staircase form, which needs a large volume of operations if these blocks have a high index of nilpotence, i.e., if their Jordan forms involve high-order blocks. According to the estimates given in [Kågström and Ruhe [22]], in addition to the reduction into Schur form, the sorting of the eigenvalues, the clustering of multiple eigenvalues, and the reduction of the diagonal blocks of the matrix to staircase form require where one flop consists of one addition and one multiplication of floating point numbers. (note that, one flop in complex arithmetic is equivalent to four ordinary flops and the corresponding number of operations should be multiplied by four.) For (nonderogatory but defective matrix with one Jordan block), this volume of work is of order and for (derogatory but nondefective matrix with scalar Jordan blocks) it is approximately . This way, for matrices with large Jordan blocks, the reduction to Jordan–Schur form requires a volume of work of order . Some increase of the efficiency of the reduction to staircase form is proposed, for instance, in [Guglielmi et al. [27]].
Example 1. Consider the matrixThe following matrix has the Jordan form:The eigenvalue has a Weyr characteristic equal to and Segre characteristic equal to (2). The eigenvalue has a Weyr characteristic equal to and the same Segre characteristic.
Using a unitary similarity transformation with a matrix (all results are rounded to eight decimal digits),The matrix is reduced by the modification of the algorithm of Kågström and Ruhe to the Jordan–Schur form.The reduction is carried out with a relative backward error.
3. Description of the Algorithms Proposed
The two algorithms proposed in this paper are modifications of the Schur–Parlett and Schur–Fréchet methods for computing the matrix exponential which instead of the Schur form uses the Jordan–Schur form. The implementation of the Jordan–Schur formin which the diagonal blocks have distinct eigenvalues , simplifies the corresponding computations and increases the accuracy. More specifically, the benefit of using the Jordan–Schur form is twofold.
Firstly, the diagonal blocks have the following structure:where to simplify the notation, the Weyr characteristic associated with is denoted by
The corresponding Segre characteristic byand . Due to the structure of these blocks, their exponentials can be found by finite Taylor series accessing only the nonzero blocks .
Secondly, the Jordan–Schur form is in upper block triangular form with the numerical multiple eigenvalues already clustered into blocks. As shown in the following theory, this facilitates significantly the usage of the Schur–Parlett and Schur–Fréchet method. Moreover, due to the special structure of the matrices and , the solution of the Sylvester matrix equations arising in determining the off-diagonal blocks of by both methods, can be carried out more efficiently in comparison with the original methods.
3.1. Determine the Exponentials of the Diagonal Blocks
Denote by the nilpotent part of a diagonal block, which has an index of nilpotence equal to . To determine the exponentials of the diagonal blocks, one exploits the following relationship:
Thus, the problem of computing reduces to the evaluation of the exponential of a nilpotent matrix with a staircase structure. The Taylor series expansion,used to determine the exponential contains only the first powers of (if one doesn’t count ). Note that, only in the extreme case of one Jordan block associated with of order when and (nonderogatory but defective matrix). In the other extreme case of (derogatory but nondefective matrix), the index of nilpotence is equal to 1 and the exponential is found directly from
For , one has that so that the raising of to power is not necessary. In the generic case, it is fulfilled that which allows obtaining the exponential with a few terms of the Taylor series expansion. Apart from increasing efficiency, this helps avoid some numerical difficulties related to the decision about the termination of the expansion. Since the nilpotent matrix is in staircase form, it is possible to achieve further reduction of the number of computational operations in evaluating the powers of this matrix. Denoting by the th power of , one has that
Hence, the power of the nilpotent matrix with staircase structure can be evaluated from for using only the nonzero blocks of and of . In computing the powers of the nilpotent matrix, one can use the scaling and squaring technique which makes the norm of less than one. The scaling is necessary in case of the large norm of and small values of the index of nilpotence when the norm of powers decreases quickly with which is related to cancelations in the floating point arithmetic and the introduction of large rounding errors. The scaling and squaring exploit the relationship where and is chosen so that . Hence,
Note that, the squaring of the matrix can be carried out taking into account that this matrix is upper triangular.
Under the assumptions for a matrix with equal eigenvalue multiplicities and equal elements of the Weyr characteristic , one obtains the following overestimate of the number of operations necessary to compute the exponentials of diagonal blocks,
The maximum volume of operations is obtained for . Hence, for (this is the case of a single Jordan block of size ), the estimated volume of operations is about and for it is equal to . Note that, the computation of the first powers of a full triangular matrix requires .
3.2. Determining the Off-Diagonal Blocks by the Schur–Parlett Method
In [Parlett [28]], Parlett developed a recurrent method for finding the off-diagonal elements of a matrix function using the upper triangular (Schur) form of a given matrix , where is unitary. The method is based on the commutativity relation and requires knowing the diagonal elements of . Since is upper triangular, so is which allows to calculate a superdiagonal of at a time starting with . A disadvantage of the Parlett recurrence is that it breaks down in the presence of equal diagonal elements (eigenvalues) and of and can produce large errors in case of close eigenvalues. That is why the idea is to use the block-diagonal form of where the diagonal blocks have distinct eigenvalues. In such case, the off-diagonal blocks of can be computed by using the following block-recurrence.
The Schur–Parlett method [Davis and Higham [12]] is as follows:
The recurrence (24) represents a Sylvester matrix equation with respect to . To have a well-conditioned Sylvester equation (24), the eigenvalues of the blocks and should be well separated for each . Similarly to the scalar case, it is necessary to know in advance the diagonal blocks . The computation of the exponential of the diagonal blocks can be carried out efficiently by using Padé approximations combined with scaling and squaring [Higham [5]], [Al-Mohy and Higham [6]].
The use of the Jordan–Schur form of facilitates significantly the usage of Parlett block recurrence (24). The exponential of the Jordan–Schur form is determined aswhere the diagonal blocks are already evaluated.
The computation of the off-diagonal blocks of the matrix exponential implementing the Parlett block recurrence involves the solution of the Sylvester equations.where the blocks and are of the form (15), are the superdiagonal blocks of the Jordan–Schur form (14) and
Due to the specific structure of the blocks and , equation (26) will be called the staircase Sylvester equation. It can be solved efficiently taking into account the structure of the diagonal blocks and .
Consider the staircase Sylvester equationwhere for simplification we use the following notation.and . It is assumed that and have eigenvalues and , respectively, with the associated Weyr characteristics.
The staircase Sylvester equation (28) can be solved by the Bartels–Stewart algorithm [Bartels and Stewart [29]]. However, in the given case, it is advantageous to exploit the staircase structure of the matrices and using Algorithm 1 described in [Petkov [26]]. Briefly, this algorithm works as follows. Partitioning the right-hand side and the unknown matrix aswhere , the blocks of the matrix are obtained as follows:
If the elements and of the Weyr characteristics are large, then this algorithm finds at once a large block instead of a block of size . Hence, in the cases of multiple eigenvalues and large elements of the Weyr characteristics, the staircase algorithm is very efficient. In fact, the solution of the staircase Sylvester equation (28) in which the matrix is , and the matrix is by the staircase algorithm requires
For large indices of nilpotency, all elements of the Weyr characteristics are equal to 1 and the staircase algorithm reduces to the Bartels–Stewart algorithm which requires
By decreasing the index of nilpotency and increasing the elements of the Weyr characteristics, the computational work associated with the staircase algorithm decreases.
Finally, having computed , the exponential of is determined aswhere is the matrix of the unitary transformation of to .
The volume of the computational work associated with the determination of the off-diagonal blocks of is estimated as follows. Assuming again equal eigenvalue multiplicities and equal Weyr elements of the characteristics , and using the estimate (33), the necessary computational work for solving all Sylvester equations of the type (26) is estimated as
The maximum of this volume of work is obtained for and for large it is approximately equal to .
Example 2. Consider again the matrixwhose Jordan–Schur form was given in Example 1. This form has a diagonal blockassociated with , a diagonal block,associated with , and a superdiagonal block,The exponential of the first diagonal block is computed by the algorithm presented in Section 3.1 with scaling equal to . The result isThe exponential of the second diagonal block, evaluated with scaling , isThe solution of the staircase Sylvester (24), carried out by the algorithm described above, isThe computed matrix exponential isThe exponential is computed with a relative error.
3.3. Determining the Off-Diagonal Blocks Using the Schur–Fréchet Method
Consider a block-triangular matrix of the formand assume that the blocks and are scaled so that their norms are less than 1. The exponential of has the following form:
Suppose that one has already calculated and . Then, the block can be computed by the following algorithm, implementing an 8th-order Padé approximation to the function . (See [Higham [5]] for a detailed exposition of the corresponding method.)
The Schur–Fréchet method [Kenney and Laub [15]].(1)Let be the zeros of the polynomial (2)Let be the zeros of the polynomial (3)Set (4)Evaluate the “Sylvester cascade” (5)(6)
In exact arithmetic, this algorithm produces an approximation satisfying [Kenney and Laub [15]]. Note that, the zeros and are purely imaginary. The algorithm requires the solution of eight Sylvester equations.
Thus, the Schur–Fréchet method determines recursively an approximation of the off-diagonal block of the exponential. The recursion involves the solution of the Sylvester equations.with matrices , and the initial approximation . Note that, is a single diagonal block of order (the multiplicity of the eigenvalue ) but the matrix may have on its diagonal several blocks corresponding to the distinct eigenvalues of multiplicities .
The key observation here is that the matrices , , , and , respectively, are in Jordan–Schur form. This allows for solving equation (48) efficiently taking into account the structure of these matrices. More specifically, in simplified notation, one has to solve the Sylvester equation.where and is a staircase matrix with Weyr characteristic associated with , the diagonal block corresponds to a distinct numerically multiple eigenvalue and is a staircase matrix with Weyr characteristic associated with .
Representing the matrices and as
One obtains the system of staircase Sylvester equations.
It is important that (51) is a staircase Sylvester equation of the same type as (28).
Having an algorithm to solve (48), it is possible to find an off-diagonal block of the matrix exponential of utilizing the Schur–Fréchet method.
Consider now, an matrixwhich is in Jordan–Schur form and assumes that the exponentials of the diagonal blocks are already determined. Using the Schur–Fréchet method, the off-diagonal blocks ofcan be evaluated recursively in the following way. First one can find the block using the blocks and . Then, using the diagonal blocks and
It is possible to evaluate the off-diagonal blocks and so on, finding ultimately the off-diagonal blocks in the first row of .
The volume of computational work necessary to determine the off-diagonal blocks by the Schur–Fréchet method is estimated as follows. Assuming as in the case of the Schur–Parlett method, equal eigenvalue multiplicities , equal elements of the Weyr characteristics , and using the estimate (33), the necessary computational work for solving eight times the Sylvester equation of the type (48) by Schur–Fréchet algorithm is estimated as
The maximum of this volume of work is obtained for and for large it is approximately equal to . Thus, the Schur–Fréchet method requires four times more work than the Schur–Parlett method provided the Jordan–Schur form is known. Note that, the reduction to this form may require a volume of work proportional to , so that the difference between the number of operations for the two methods may be negligible in the total cost.
Example 3. Consider again the matrix used in Examples 1 and 2. Using the Schur–Fréchet algorithm implementing the Jordan–Schur form, the matrix exponential is computed with a relative error.
4. Numerical Considerations
The error in the evaluation of the matrix exponential depends on the conditioning of the exponential and the round-off properties of the numerical algorithm used.
The sensitivity of the matrix exponential to perturbations in is studied by several authors, see, for instance, [Van Loan [30]], [Kenney and Laub [31]], and [Mathias [32]]. The relative condition number of the matrix exponential that characterizes its sensitivity to small perturbations of , is defined as [Higham [11]; Ch. 3].
This condition number is determined by using the Fréchet derivative of the exponential whose computation can be carried out by Algorithm 10.27 presented in [Higham [13]; Ch. 10], or by the algorithms given in [Mathias [33]] and [Al-Mohy and Higham [34]].
4.1. The Reduction to Jordan–Schur Form
As mentioned earlier, the reduction to the Jordan–Schur form is numerically stable due to the implementation of unitary transformations. However, it should be emphasized that the numerical stability does not guarantee that the clustering of numerically multiple eigenvalues and the determination of the size of the Jordan blocks are always done in an optimal way. The extensive experiments with this algorithm show that it works reliably except in the case of a very ill-conditioned system of eigenvectors and principal vectors when any clustering method will have serious difficulties in grouping the eigenvalues.
4.2. Computing the Off-Diagonal Blocks by the Schur–Parlett Method
In the present version of the Schur–Parlett algorithm, the exponentials of the diagonal blocks are found by using the finite term Taylor series combined with scaling and squaring. As concluded by Higham [Higham [11]; Ch. 10], no clear general conclusions can be drawn about the stability of the scaling and squaring algorithm. Nevertheless, this algorithm is generally considered stable and is implemented, for instance, in the MATLAB™ function expm.
The next step of the algorithm is the Parlett block recursion which requires the solution of a Sylvester equation in the form (26) with respect to the unknown off-diagonal block . In the case of close eigenvalues of the blocks and , this equation can be ill-conditioned. More specifically, let the Sylvester equation be written aswhere and let . Then, the relative change in the solution of this equation due to small relative perturbations and , is given by [Higham [35]; Ch. 16].whereis the relative condition number of the Sylvester equation andis the maximum relative perturbation in the data. Now, if the separation,of the matrices and [Stewart [36]] is small, then the condition number of the Sylvester equation is large and the sensitivity of this equation can be large even in the case of a well-conditioned matrix exponential. (This follows from the fact that ). Hence, the Jordan–Schur algorithm implementing the Schur–Parlett method may behave unstably in the presence of blocks in the Jordan–Schur form that has a small separation between them. Note that, this instability is common for all methods implementing the Parlett block recurrence and is not specific to the use of the Jordan–Schur form itself.
4.3. Computing the Off-Diagonal Blocks by the Schur–Fréchet Method
The main computational task in the Schur–Fréchet method is the solution of the system of staircase Sylvester equations (the “Sylvester cascade”) (48). Unlike the Sylvester equation (26), solved by the Schur–Parlett method, the sensitivity of (48) does not depend on the closeness of the eigenvalues of Jordan–Schur blocks. Moreover, the separation of matrices and in (51) is always of order 1 which leads to a low condition number of the Sylvester equation. For this reason, the Schur–Fréchet algorithm almost always produces more accurate results than the Schur–Parlett algorithm. The Schur–Fréchet algorithm is especially appropriate in cases when the blocks of the Jordan–Schur form have close eigenvalues.
4.4. Computational Complexity
The number of computational operations in complex flops required by the Schur–Parlett and Schur–Fréchet algorithm is summarized in Tables 1 and 2. The estimates of the computation cost are obtained under the assumptions .
In Table 1, we give the approximate number of operations necessary for reduction into Jordan–Schur form and computation of the diagonal blocks of the exponential. These computations are common for both methods and the corresponding volume of work is given for the two extreme cases: Derogatory but nondefective matrix with scalar Jordan blocks Nonderogatory but defective matrix with a single Jordan block
As seen from Table 1, the number of operations in these cases varies between and .
In Table 2, we give the number of operations for the Schur–Parlett and Schur–Fréchet algorithm including the extreme cases and . The third column of the table contains the total sum of operations necessary for reduction into the Jordan–Schur form, computing the exponentials of the diagonal blocks and computing the off-diagonal blocks by the corresponding method. Clearly, the work necessary for the reduction into Jordan–Schur form dominates the total operation count which is approximately the same for both methods. In the most difficult case , both methods require flops.
5. Numerical Experiments
In assessing the numerical behavior of the algorithms for computing the matrix exponential, the authors commonly use the test matrices included in the sets described in [Higham [37]], [Fasi and Higham [38]], and [Wright [39]]. Many of the matrices from these sets have distinct eigenvalues so that their Jordan–Schur form coincides with the Schur form. In such cases, the Jordan–Schur algorithms for computing the exponential do not have some advantages over the Schur algorithms. That is why in this section we present four experiments involving matrices with different Jordan structures so that the Jordan–Schur form differs significantly from the Schur form.
The computations in the paper are carried out with MATLAB™ Version 9.9 (R2020b) [40] using IEEE double precision arithmetic. The tolerances and , used in the reduction to Jordan–Schur form, are both taken equal to . The relative condition number of the matrix exponential is determined by the function expm_cond from the matrix function toolbox [Higham [13]]. The computations are done using the Schur–Parlett algorithm, the Schur–Fréchet algorithm, the function expm in MATLAB™ implementing Padé approximations combined with scaling and squaring [Al-Mohy and Higham [6]], and the function funm in MATLAB™ which implements the Schur–Parlett algorithm [Davis and Higham [12]].
Experiment 1. The first experiment demonstrates the behavior of the two versions of the Jordan–Schur algorithm for a family of problems with nonderogatory matrices and with increasing conditioning of the exponential. In the given case, the 8th-order matrix is constructed aswhereis the desired Jordan form and is a nonsingular matrix. The matrix is taken aswhere the matrix is chosen asThe parameter varies between 0 and 30 and the condition number of the matrix exponential varies between and , respectively. An accurate approximation of the exact exponential of is determined as , where is computed using the corresponding analytical expression.
The results of this experiment are shown in Figure 1. The Schur–Parlett and the Schur–Fréchet algorithm produce almost identical results, close to the results obtained by the function funm.
Experiment 2. This experiment illustrates the difficulties with the Schur–Parlett version of the algorithm in case of a small separation between the diagonal blocks of the Jordan–Schur form.
The matrix in the experiment is constructed as in Experiment 1,with the matrix taken asNote that, the matrices and have integer entries. The condition number with respect to the inversion of the similarity transformation matrix is . The Jordan form is chosen with two blocks. The first block is of 5th order and its eigenvalue varies between 3 and 3.9, while the second block is of 3rd order with an eigenvalue which is constant and equal to 4. As in Example 1, the exact exponential is obtained using the analytical expression of the exponential . With the increasing of , the separation between the two blocks and decreases quickly, which leads to large errors in the solution of Sylvester (26) appearing in the block-recurrence. The relative error in computing the exponential along with the error bound where is the MATLAB™ function eps, , as a function of is shown in Figure 2.
The experiment presented shows the potential instability of the Schur–Parlett method and its Jordan–Schur version in the case of blocks with close eigenvalues. The function funm also produces large errors, as a result of implementing the Schur–Parlett method. For such matrices, it is appropriate to use the Schur–Fréchet algorithm, the Padé approximation combined with scaling and squaring, or the alternative methods that are based on the Taylor series [Caliari and Zivcovich [8]].
Experiment 3. This experiment illustrates the properties of the algorithms proposed in the case of matrices with defective eigenvalues and large elements of the Weyr characteristic. The matrix is of even order and is taken aswhere the Jordan form.Consists of blocks with equal eigenvalues and the nonsingular transformation matrix is constructed as proposed in [Bavely and Stewart [41]],where are elementary reflections that are orthogonal and symmetric matrices [Stewart [42]]. The condition number of with respect to the inversion is controlled by the variable and is equal to . Note that, the matrices and are obtained with small rounding errors. An accurate approximation of the exact exponential of is determined from , where is evaluated using its analytical expression. In the given case, we take and the parameter controlling the conditioning of is chosen equal to 1.25. The Segre characteristic of associated with is equal to and the corresponding Weyr characteristic is .
The relative errors of the computed exponential for are shown in Figure 3 for the four methods used earlier. The function expm exhibits the worst performance and for the relative error of the computed result is equal to . At the same time, the relative errors of the exponential computed by the two versions of the Jordan–Schur algorithm are identical and for are approximately equal to . Even for , the relative errors of these algorithms are smaller than 1. The good results for the Jordan–Schur algorithms are due to the simple structure of the Jordan–Schur form.which facilities very much the computation of the exponential. The function funm also produces accurate results for this example.
Experiment 4. In this experiment, the matrix is constructed in the same way as in Experiment 3 but in the Jordan form.It consists of blocks of order 2 whose eigenvalues are equal to and one single Jordan block that has same eigenvalue. In the given case, we take and the parameter controlling the conditioning of the nonsingular transformation matrix is chosen equal to 1.2.
The relative error of the computed exponential for different is shown in Figure 4. In the given case, the Jordan–Schur algorithms and the function expm produce results that are bounded by the quantity while the function funm exhibits instability for . Note that, the Segre characteristic of associated with is and the corresponding Weyr characteristic is .
6. Conclusions
The two algorithms proposed in this paper take advantage of using the Jordan–Schur form for finding the matrix exponential in the case of multiple eigenvalues. The algorithms are more accurate than some known methods for computing the matrix exponential in certain cases at the price of increasing the computational cost of the overall algorithm due to the reduction to the Jordan–Schur form. The numerical experiments confirm that the Jordan–Schur algorithms for computing the matrix exponential are appropriate for matrices with multiple eigenvalues and are especially efficient in cases of large elements of the Weyr characteristics.
Data Availability
The datasets generated during the current study are available from the author on reasonable request.
Conflicts of Interest
The author declares that there are no conflicts of interest.
Supplementary Materials
M-files implementing the algorithms and test examples are available at https://www.dropbox.com/sh/wxns1j0eab8flmu/AAD_01G6xqn8bUiW_bV3dgnma?dl=0. (Supplementary Materials)