Next Article in Journal
Numerical and Evolutionary Optimization 2020
Previous Article in Journal
Odd Exponential-Logarithmic Family of Distributions: Features and Modeling
Previous Article in Special Issue
An Efficient Numerical Scheme Based on Radial Basis Functions and a Hybrid Quasi-Newton Method for a Nonlinear Shape Optimization Problem
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

On Some Numerical Methods for Solving Large Differential Nonsymmetric Stein Matrix Equations

Laboratory of Engineering Sciences for Energy, National School of Applied Sciences of El Jadida, University Chouaib Doukkali, Av. des Facultés, El Jadida 24000, Morocco
*
Author to whom correspondence should be addressed.
Math. Comput. Appl. 2022, 27(4), 69; https://doi.org/10.3390/mca27040069
Submission received: 30 June 2022 / Revised: 5 August 2022 / Accepted: 9 August 2022 / Published: 12 August 2022
(This article belongs to the Special Issue Computing and Green Technology)

Abstract

:
In this paper, we propose a new numerical method based on the extended block Arnoldi algorithm for solving large-scale differential nonsymmetric Stein matrix equations with low-rank right-hand sides. This algorithm is based on projecting the initial problem on the extended block Krylov subspace to obtain a low-dimensional differential Stein matrix equation. The obtained reduced-order problem is solved by the backward differentiation formula (BDF) method or the Rosenbrock (Ros) method, the obtained solution is used to build the low-rank approximate solution of the original problem. We give some theoretical results and report some numerical experiments.

1. Preliminaries

In the present paper, we are interested in the numerical solution of large-scale differential nonsymmetric Stein matrix equations on the time interval [ t 0 , T f ] of the form:
X ˙ ( t ) = A X ( t ) B X ( t ) + E F T , t [ t 0 , T f ] , ( D N S E ) ; X ( t 0 ) = X 0 ,
where A and B are real, sparse and square matrices of size n × n and p × p , respectively, and E and F are matrices of size n × r and p × r , respectively, with r n , p . The initial condition is given in a factored form as X 0 = Z 0 Z ˜ 0 T R n × p and the matrices A and B are assumed to be large and sparse.
Differential Stein matrix equations play an important role in many problems in control and filtering theory for discrete-time large-scale dynamical systems and other problems; see [1,2,3,4,5,6,7,8]. For solving the large matrix equations during the last years, several projection methods based on Krylov subspaces have been proposed, for example Sylvester matrix equations, differential Sylvester matrix equations, Riccati and others, see, e.g., [8,9,10,11,12,13,14,15,16]. The main idea employed in these methods is to project the initial problem by using extended Krylov subspaces and then apply the Petrov–Galerkin orthogonality condition.
In this work, we present an extended block Krylov method for solving large differential Stein matrix equations. Our method uses the Petrov–Galerkin projection method and the extended block Arnoldi algorithm. The problem (1) is projected onto small extended Krylov subspaces to obtain low-order differential Stein equations that are solved by a BDF method or Ros method. The approximate solution is then given as a product of matrices with low rank.
We recall the extended block Arnoldi algorithm applied to ( A , V ) where V R n × r . The extended block Krylov subspace associated to ( A , V ) (see [8,9,14,15,17]) is given by:
K m e ( A , V ) = R a n g e { V , A 1 V , A V , A 2 V , A 2 V , , A m 1 V , A m V } ) = K m ( A , V ) + K m ( A 1 , A 1 V ) ,
where
K m ( A , V ) = R a n g e { V , A V , A 2 V , , A m 1 V } .
Algorithm 1 allows us to construct an orthonormal matrix V m = V 1 , V 2 , , V m that is a basis of the block extended Krylov subspace K m e ( A , V ) . Let T m , A the restriction of the matrix A to the block extended Krylov subspace K m e ( A , V ) is given by T m , A = V m T A V m . Then, we have the following relations (see [15])
A V m = V m T m , A + V m + 1 T m + 1 , m A E m T .
where E m = [ 0 2 r × 2 ( m 1 ) r , I 2 r ] T and T m , A = [ T i , j A ] .
Algorithm 1: The extended block Arnoldi algorithm (EBA)
 1.  Inputs: A an n × n matrix, V an n × r matrix and m an integer.
 2.  Compute the QR decomposition of [ V , A 1 V ] , i.e., [ V , A 1 V ] = V 1 Λ ;
 3.  Set V 0 = ;
 4.  For j = 1 , 2 , . . . , m
    (a)  Set V j ( 1 ) : first r columns of V j ; V j ( 2 ) : second r columns of V j
    (b)   V j = V j 1 , V j ; V ^ j + 1 = A V j ( 1 ) , A 1 V j ( 2 ) ;
    (c)  Orthogonalize V ^ j + 1 i.e.,
      *  for i = 1 , 2 , , j
      *   H i , j = V i T V ^ j + 1 ;
      *   V ^ j + 1 = V ^ j + 1 V i H i , j ;
      *  End for
    (d)  Compute the Q R decomposition of V ^ j + 1 , i.e., V ^ j + 1 = V j + 1 H j + 1 , j ;
 5.  End For.
Throughout the paper, we use the following notations. The Frobenius inner product of the matrices X and Y is defined by X , Y F = t r ( X T Y ) , where t r ( Z ) denotes the trace of a square matrix Z. The associated norm is the Frobenius norm denoted by F . The Kronecker product A B = [ a i , j B ] R n p × n p where A = [ a i , j ] . This product satisfies the properties: ( A B ) ( C D ) = ( A C B D ) and ( A B ) T = A T B T .
The rest of the paper is organized as follows. In Section 2, we derive the low-rank approximate solutions method and give some theoretical results. In Section 3, we apply the backward differentiation formula method and Rosenbrock method for solving reduced order problem. Section 4 is devoted to some numerical examples.

2. Low-Rank Approximate Solutions Method

In this section, we project the large differential equation by using the extended block Krylov subspaces K m e ( A , E ) and K m e ( B T , F ) to obtain low-rank approximate solutions of Equation (1).
We apply the extended block Arnoldi algorithm (Algorithm 1) to the pairs ( A , E ) and ( B T , F ) , respectively, and to obtain orthonormal bases V 1 , . . . , V m and W 1 , . . . , W m and we have
T m , A = V m T A V m and T m , B = W m T B T W m ,
where
V m = [ V 1 , V 2 , . . . , V m ] and W m = [ W 1 , W 2 , . . . , W m ] .
We then consider approximate solutions of the large differential Stein matrix Equation (1) that have the low-rank form
X m ( t ) = V m Y m ( t ) W m T .
The matrix Y m ( t ) is determined from the following Petrov–Galerkin orthogonality condition:
V m T R m ( t ) W m = 0 , t [ t 0 , T f ] ,
where R m ( t ) is the residual
R m ( t ) = X ˙ m ( t ) A X m ( t ) B + X m ( t ) E F T .
Then, from (4) and (3), we obtain the low dimensional differential Stein equation
Y ˙ m ( t ) = T m , A Y m ( t ) T m , B T Y m ( t ) + E ˜ m F ˜ m T , Y m ( t 0 ) = V m X 0 W m T ,
where E ˜ m = V m T E and F ˜ m = W m T F .
In the following result, we derive an expression for computation of the norm of the residual R, without having to compute matrix products with the large matrices A and B. This result allows us to reduce the cost in the proposed method when checking if the residual norm is less than some fixed tolerance.
Theorem 1. 
Let X m ( t ) = V m Y m ( t ) W m T be the approximation obtained at step m by the extended block Arnoldi method, and Y m ( t ) the exact solution of low dimensional differential nonsymmetric Stein Equation (5). Then, the residual R m ( t ) associated to X m ( t ) satisfies the relation
R m ( t ) F 2 = T m + 1 , m A Y ¯ m , 1 ( t ) T m , B T F 2 + T m + 1 , m B Y ¯ m , 2 ( t ) T m , A T F 2 + T m + 1 , m A Y m , m ( t ) ( T m + 1 , m B ) T F ,
where Y ¯ m , 1 ( t ) is the 2 r × 2 m r matrix corresponding to the last 2 r rows of Y m ( t ) , Y ¯ m , 2 ( t ) = E m T Y m ( t ) T and Y m , m ( t ) = E m T Y m ( t ) E m .
Proof. 
Using the relations (4) and (2), we have
R m ( t ) = V m Y ˙ m ( t ) W m T A V m Y m ( t ) W m T B + V m Y m ( t ) W m T E F T .
By using the following relations
A V m = V m + 1 T m , A T m + 1 , m A E m T , W m T B = T m , B T E m ( T m + 1 , m B ) T W m + 1 T , V m = V m + 1 I 2 r m 0 2 r × 2 r , W m T = I 2 r m 0 2 r × 2 r W m + 1 T ,
we have R m ( t ) = V m + 1 I 2 r m 0 2 r × 2 r Y ˙ m ( t ) I 2 r m 0 2 r × 2 r W m + 1 T V m + 1 T m , A T m + 1 , m A E m T
Y m ( t ) T m , B T E m ( T m + 1 , m B ) T W m + 1 T + V m + 1 I 2 r m 0 2 r × 2 r Y m ( t ) I 2 r m 0 2 r × 2 r W m + 1 T E F T , since E F T = V m + 1 E ˜ m F ˜ m T 0 2 m r × 2 r 0 2 r × 2 m r 0 2 r × 2 m r W m + 1 T , so
R m ( t ) = V m + 1 Y ˙ m ( t ) T m , A Y m ( t ) T m , B T + Y m ( t ) E ˜ m F ˜ m T T m , A Y m ( t ) E m ( T m + 1 , m B ) T T m + 1 , m A E m T Y m ( t ) T m , B T T m + 1 , m A E m T Y m ( t ) E m ( T m + 1 , m B ) T W m + 1 T .
Now, as Y m ( t ) exact solution of the low dimensional differential Stein equation
Y ˙ m ( t ) = T m , A Y m ( t ) T m , B T Y m ( t ) + E ˜ m F ˜ m T ,
so
R m ( t ) = V m + 1 0 T m , A Y m ( t ) E m ( T m + 1 , m B ) T T m + 1 , m A E m T Y m ( t ) T m , B T T m + 1 , m A E m T Y m ( t ) E m ( T m + 1 , m B ) T W m + 1 T .
and since V m + 1 and W m + 1 are orthonormal, we obtain
R m ( t ) F 2 = T m + 1 , m A Y ¯ m , 1 ( t ) T m , B T F 2 + T m + 1 , m B Y ¯ m , 2 ( t ) T m , A T F 2 + T m + 1 , m A Y m , m ( t ) ( T m + 1 , m B ) T F .
where Y ¯ m , 1 ( t ) = E m T Y m ( t ) , Y ¯ m , 2 ( t ) = E m T Y m ( t ) T and Y m , m ( t ) = E m T Y m ( t ) E m .    □
The approximate solution X m ( t ) can be given as a product of two matrices of low rank. Consider the singular value decomposition of the 2 m r × 2 m r matrix:
Y m ( t ) = Y ˜ 1 Σ Y ˜ 2 T ,
where Σ is the diagonal matrix of the singular values of Y m sorted in decreasing order. Let Y 1 , l and Y 2 , l be the 2 m r × l matrices of the first l columns of Y ˜ 1 and Y ˜ 2 , respectively, corresponding to the l singular values of magnitude greater than some tolerance. We obtain the truncated singular value decomposition:
Y m U 1 , l Σ l U 2 , l T ,
where Σ l = diag [ σ 1 , , σ l ] . Setting Z 1 , m = V m U 1 , l Σ l 1 / 2 , and Z 2 , m = W m U 2 , l Σ l 1 / 2 , it follows that:
X m Z 1 , m ( t ) Z 2 , m T ( t ) .
This is very important for large problems when one does not need to compute and store the approximation X m at each iteration.
The following result shows that the approximation X m ( t ) is an exact solution of a perturbed differential Stein equation.
Theorem 2. 
Let X m ( t ) be the approximate solution given by (2). Then, we have
X ˙ m ( t ) = ( A F m , A ) X m ( t ) ( B F m , B ) X m ( t ) + E F T ,
where
F m , A = V m + 1 T m + 1 , m A E m T V m T , F m , B = W m E m ( T m + 1 , m B ) T W m + 1 T .
Proof. 
By multiplying the (5) left by V m and right by W m T , we obtain
V m Y ˙ m ( t ) W m T = V m T m , A Y m ( t ) T m , B T W m T V m Y m ( t ) W m T + V m E ˜ m F ˜ m T W m T .
Using relationships
A V m = V m T m , A + V m + 1 T m + 1 , m A E m T , W m T B = T m , B T W m T + E m ( T m + 1 , m B ) T W m + 1 T ,
we find
X ˙ m ( t ) = [ A V m V m + 1 T m + 1 , m A E m T ] Y m ( t ) [ W m T B E m ( T m + 1 , m B ) T W m + 1 T ] X m ( t ) + E F T = [ A V m + 1 ( T m + 1 , m A E m T V m T ] X m ( t ) [ B W m E m ( T m + 1 , m B ) T W m + 1 T ] X m ( t ) + E F T .
So
X ˙ m ( t ) = ( A F m , A ) X m ( t ) ( B F m , B ) X m ( t ) + E F T ,
where F m , A = V m + 1 T m + 1 , m A V m T and F m , B = W m ( T m + 1 , m B ) T W m + 1 T .    □
The next result states that the error E m ( t ) = X ( t ) X m ( t ) satisfies also a differential Stein matrix equation.
Theorem 3. 
Let X ( t ) , the exact solution of (1), and X m ( t ) be the low-rank approximate solution obtained at step m. The error E m ( t ) satisfies the following differential Stein equation
E ˙ m ( t ) = A E m ( t ) B E m ( t ) R m ( t ) .
Proof. 
According to (1) and (4) we have
E ˙ m ( t ) = X ˙ ( t ) X ˙ m ( t ) = A X ( t ) B X ( t ) + E F T A X m ( t ) B + X m ( t ) E F T R m ( t ) = A ( X ( t ) X m ( t ) ) B ( X ( t ) X m ( t ) ) R m ( t ) = A E m ( t ) B E m ( t ) R m ( t ) .
   □
Next, we give an upper bound for the norm of the error X X m by using the 2-logarithmic norm defined by μ 2 ( A ) = 1 2 λ max ( A + A T ) . The logarithmic norm satisfies the following relation
e t A 2 e μ 2 ( A ) t ; t 0 .
Theorem 4. 
At step m of the extended block Arnoldi process, we have the following upper bound for the norm of the error,
E m ( t ) F e ( t t 0 ) μ 2 ( A ) E m ( t 0 ) F + max τ [ t 0 , t ] R m ( τ ) F e ( t t 0 ) μ 2 ( A ) 1 μ 2 ( A ) .
Proof. 
We notice that the differential Stein Equation (9) is equivalent to the linear ordinary differential equation
e r r ˙ m ( t ) = A e r r m ( t ) b m ( t ) , e r r 0 = v e c ( E m ( t 0 ) ) ,
where
A = B T A I n p , e r r m ( t ) = v e c ( E m ( t ) ) , b m ( t ) = v e c ( R m ( t ) ) ,
Then, the error e r r m ( t ) can be expressed in the integral form as follows
e r r m ( t ) = e ( t t 0 ) A e r r 0 t 0 t e ( t τ ) A b m ( τ ) d τ , t [ t 0 , T f ] .
By using e t A 2 e μ 2 ( A ) t , we have
e r r m ( t ) 2 e ( t t 0 ) A e r r 0 t 0 t e ( t τ ) A b m ( τ ) d τ 2 e ( t t 0 ) A e r r 0 2 + t 0 t e ( t τ ) A b m ( τ ) d τ 2 e ( t t 0 ) μ 2 ( A ) e r r 0 2 + t 0 t e ( t τ ) μ 2 ( A ) b m ( τ ) 2 d τ e ( t t 0 ) μ 2 ( A ) e r r 0 2 + max τ [ t 0 , t ] b m ( τ ) 2 t 0 t e ( t τ ) μ 2 ( A ) d τ e ( t t 0 ) μ 2 ( A ) e r r 0 2 + max τ [ t 0 , t ] b m ( τ ) 2 e ( t t 0 ) μ 2 ( A ) 1 μ 2 ( A ) e ( t t 0 ) μ 2 ( A ) e r r 0 2 + max τ [ t 0 , t ] v e c ( R m ( τ ) ) 2 e ( t t 0 ) μ 2 ( A ) 1 μ 2 ( A ) .
As v e c ( E m ( t ) ) 2 = E m ( t ) F , so
E m ( t ) F e ( t t 0 ) μ 2 ( A ) E m ( t 0 ) F + max τ [ t 0 , t ] R m ( τ ) F e ( t t 0 ) μ 2 ( A ) 1 μ 2 ( A ) .
   □

3. Solving the Projected Differential Stein Matrix Equation

3.1. Rosenbrock Method

In this section, we are applying the Ros method [18] for solving the projected differential Stein matrix equation. We can write the low dimensional nonsymmetric differential Stein Equation (5) in the following form
Y m ˙ ( t ) = S ( Y m ( t ) ) , t [ t 0 , T f ] Y m ( t 0 ) = Y m ( 0 ) ,
where S ( Y m ) = T m , A Y m T m , B T Y m + E ˜ m F ˜ m T and Y m ( 0 ) = V m X 0 W m T .
The approximation Y m , k + 1 of Y m ( t k + 1 ) obtained at step k + 1 by Ros method is given by
Y k + 1 = Y k + 3 2 P 1 + 1 2 P 2 ,
where P 1 and P 2 solve the following Stein matrix equations
Λ A P 1 Λ B T P 1 + S ( Y k ) = 0 ,
and
Λ A P 2 Λ B T P 2 + S ( Y k + P 1 ) 2 h P 1 = 0 ,
where
Λ A = γ T m , A + 1 h + 1 γ I s i z e ( T m , A ) , Λ B = γ T m , B + 1 h + 1 γ I s i z e ( T m , B ) ,
with γ = 1 / 2 .
The Ros algorithm for solving the reduced differential Stein matrix Equation (5) is summarized in Algorithm 2.
Algorithm 2: The 2-Rosenbrock method for solving the reduced NDSE (12)
 Input: T m , A , T m , B , E ˜ m , F ˜ m , t 0 , T f .
 1.  Choose h.
 2.  ompute: N = T f t 0 h
 3.  Compute: Λ A
 4.  Compute: Λ B
 5.  For k = 1 : N + 1
    (a)  Compute P 1 by solving Stein matrix Equation (14).
    (b)  Compute P 2 by solving Stein matrix Equation (15).
    (c)  Compute Y k + 1 by (13).
 6.  End For k.
 Output: Y T f .

3.2. Backward Differentiation Formula Method

In this section, we present a BDF method [19] for solving, at each step m of the extended block Arnoldi process, the low dimensional differential Stein matrix Equation (5).
At each time t k , let Y m , k be the approximation of Y m ( t k ) , where Y m is a solution of (5). Then, the new approximation Y m , k + 1 of Y m ( t k + 1 ) obtained at step k + 1 by BDF method is defined by the implicit relation
Y m , k + 1 = i = 0 p 1 α i Y m , k i + h β S ( Y m , k + 1 ) ,
where h = t k + 1 t k is the step size, α i and β are the coefficients of the BDF method as listed in the Table 1.
The approximate Y m , k + 1 solves the following matrix equation
Y m , k + 1 + h β ( T m , A Y m , k + 1 T m , B T Y m , k + 1 + E ˜ m F ˜ m T ) + i = 0 p 1 α i Y m , k i = 0 .
We assume that at each time t k , the approximation Y m , k is factorized as a low-rank product Y m , k U ˜ m , k V ˜ m , k T , where U ˜ m , k R n × m k and V ˜ m , k R p × m k , with m k n , p . We define
T ˜ m , A = h β T m , A h β I 2 m , T ˜ m , B = h β T m , B T + h β I 2 m , E m , k + 1 = [ h β E ˜ m T , α 0 U ˜ m , k T , . . . , α p 1 U ˜ m , k + 1 p T ] T , F m , k + 1 = [ h β F ˜ m T , α 0 V ˜ m , k T , . . . , α p 1 V ˜ m , k + 1 p T ] T .
Then, we obtain the following matrix Stein equation:
T ˜ m , A Y m , k + 1 T ˜ m , B Y m , k + 1 + E m , k F m , k T = 0 .
The low-rank approximate solutions method by extended block Arnoldi algorithm for the large differential nonsymmetric stein matrix equations is summarized as follows in Algorithm 3.
Algorithm 3: The low-rank extended block Arnoldi method for DNSE
 Input: X 0 , choose a tolerance t o l > 0 , a maximum number of m m a x iterations.
 1.  For m = , 2 , 3 , . . . , m m a x do
 2.  Update V m , T m , A , by Algorithm 1 (EBA) applied to ( A , E )
 3.  Update W m , T m , B , by Algorithm 1 (EBA) applied to ( B T , F )
 4.  Solve the low-dimensional problem (5) by BDF method or Ros method.
 5.  If R m ( t ) F < t o l , stop.
 6.  End If;
 7.  End For (m);
 8.  Using (7), the approximate solution X m is given by X m Z 1 , m Z 2 , m T .
 Output: X m .

4. Numerical Experiments

In this section, we give some numerical examples of large nonsymmetric differential Stein matrix equations. All the experiments were performed on a computer of Intel Core i5 at 1.3 GHz and 4 GB of RAM. The Algorithm 3 were coded in Matlab2014. In all of the examples, the coefficients of the matrices E and F were random values uniformly distributed on [ 0 , 1 ] . The stopping criterion used for EBA-BDF method and EBA-Ros was R ( X m ) F < 10 10 or a maximum of m m a x = 40 iterations was achieved.

4.1. Example 1

In this first example, the matrices A and B are obtained from the centered finite difference discretization of the operators:
L A ( u ) = Δ u + f 1 ( x , y ) u x + , f 2 ( x , y ) u y + f ( x , y ) u ,
L B ( u ) = Δ u + g 1 ( x , y ) u x + g 2 ( x , y ) u y + g ( x , y ) u ,
on the unit square [ 0 , 1 ] × [ 0 , 1 ] with homogeneous Dirichlet boundary conditions. The number of inner grid points in each direction was n 0 and p 0 for the operators L A and L B , respectively. The matrices A and B were obtained from the discretization of the operator L A and L B with the dimensions n = n 0 2 and p = p 0 2 , respectively. The discretization of the operator L A ( u ) and L B ( u ) yields matrices extracted from the Lyapack package [20] using the command fdm_2d_matrix and denoted as A = fdm( n 0 ,’f_1(x,y)’,’f_2(x,y)’,’f(x,y)’). In this example, n = 10,000 and p = 10,000, respectively, and are named as A = fdm ( n 0 , f 1 ( x , y ) , f 2 ( x , y ) , f ( x , y ) ) and B = fdm ( p 0 , g 1 ( x , y ) , g 2 ( x , y ) , g ( x , y ) ) with f 1 ( x , y ) = e x y , f 2 ( x , y ) = sin ( x y ) , f ( x , y ) = y 2 , g 1 ( x , y ) = 100 e x , g 2 ( x , y ) = 12 x y and g ( x , y ) = x 2 + y 2 . For this experiment, we used r = 2 .
In Figure 1, we plotted the Frobenius norms of the residuals versus the number of iterations for the EBA-BDF and the EBA-Ros method.
In Table 2, we compared the performances of the EBA-BDF method and the EBA-Ros. For both methods, we listed the residual norms, the maximum number of iteration and the corresponding execution time.

4.2. Example 2

For the second set of experiments, we use the matrices a d d 32 , p d e 2961 , and t h e r m a l from the Harwell Boeing Collection [21]. The tolerance was set to 10 8 for the stop test on the residual. For the EBA-BDF and EBA-Ros methods, we used a constant timestep h = 0.1 ,
In Figure 2, the matrices A = a d d 32 and B = p d e 2961 with dimensions n = 4960 and p = 2961 , respectively, and r = 2 . We plotted the Frobenius norms of the residuals R m ( T f ) F at final time T f versus the number of extended block Arnoldi iterations for the EBA-BDF and EBA-Ros method.
In Figure 3, the matrices A = fdm ( x + 10 y 2 , 2 x 2 + y 2 , y 2 x 2 ) and B = t h e r m a l with dimensions n = 10,000 and p = 3456 , respectively, h = 0.2 and r = 3 . We plotted the Frobenius norms of the residuals R m ( T f ) F at final time T f versus the number of extended block Arnoldi iterations for the EBA-BDF and EBA-Ros method.
In Table 3, we list the Frobenius residual norms at final time T f = 2 and the corresponding CPU time for each method.
The numerical results are promising, showing the effectiveness of the proposed methods.

5. Conclusions

We presented, in this paper, a new iterative method for computing numerical solutions for large-scale differential Stein matrix equations with low-rank right-hand sides. This approach is based on projection onto extended block Krylov subspaces with a Galerkin method. The approximate solutions are given as products of two low-rank matrices and allow for saving memory for large problems. The numerical experiments show that the proposed extended block Krylov-based method is effective for large and sparse matrices.

Author Contributions

Conceptualization, L.S., E.M.S. and H.T.A.; methodology, E.M.S.; software, L.S.; validation, L.S., E.M.S. and H.T.A.; formal analysis, L.S.; investigation, E.M.S.; writing—original draft preparation, L.S., E.M.S. and H.T.A.; writing—review and editing, L.S., E.M.S. and H.T.A. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Acknowledgments

The authors should express their deep-felt thanks to the anonymous referees for their encouraging and constructive comments, which improved this paper.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Van Dooren, P.M. Gramian Based Model Reduction of Large-Scale Dynamical Systems. Available online: https://citeseerx.ist.psu.edu/viewdoc/download;jsessionid=E53840A73CA1503CFAAF2F4EE3009768?doi=10.1.1.29.2500&rep=rep1&type=pdf (accessed on 29 June 2022).
  2. Zhou, B.; Lam, J.; Duan, G.R. On Smith-type iterative algorithms for the Stein matrix equation. Appl. Math. Lett. 2009, 22, 1038–1044. [Google Scholar] [CrossRef]
  3. Zhou, B.; Lam, J.; Duan, G.R. Toward solution of matrix equation X = Af(X) B + C. Linear Algebra Appl. 2011, 435, 1370–1398. [Google Scholar] [CrossRef]
  4. Bouhamidi, A.; Jbilou, K. A note on the numerical approximate solutions for generalized Sylvester matrix equations with applications. Appl. Math. Comput. 2008, 206, 687–694. [Google Scholar] [CrossRef]
  5. Calvetti, D.; Levenberg, N.; Reichel, L. Iterative methods for XAXB = C. J. Comput. Appl. Math. 1997, 86, 73–101. [Google Scholar] [CrossRef]
  6. Datta, B. Numerical Methods for Linear Control Systems; Academic Press: Cambridge, MA, USA, 2004; Volume 1. [Google Scholar]
  7. Datta, B.N.; Datta, K. Theoretical and computational aspects of some linear algebra problems in control theory. Comput. Comb. Methods Syst. Theory 1986, 177, 201–212. [Google Scholar]
  8. Bentbib, A.H.; Jbilou, K.; Sadek, E.M. On some extended block Krylov based methods for large scale nonsymmetric Stein matrix equations. Mathematics 2017, 5, 21. [Google Scholar] [CrossRef]
  9. Agoujil, S.; Bentbib, A.H.; Jbilou, K.; Sadek, E.M. A minimal residual norm method for large-scale Sylvester matrix equations. Electron. Trans. Numer. Anal. 2014, 43, 45–59. [Google Scholar]
  10. Sadek, E.; Bentbib, A.H.; Sadek, L.; Talibi Alaoui, H. Global extended Krylov subspace methods for large-scale differential Sylvester matrix equations. J. Appl. Math. Comput. 2020, 62, 157–177. [Google Scholar] [CrossRef]
  11. Sadek, L.; Talibi Alaoui, H. Numerical methods for solving large-scale systems of differential equations. Ric. Mat. 2021, 1–18. [Google Scholar] [CrossRef]
  12. Sadek, L.; Talibi Alaoui, H. The extended block Arnoldi method for solving generalized differential Sylvester equations. J. Math. Model. 2020, 8, 189–206. [Google Scholar]
  13. Sadek, L.; Alaoui, H.T. Application of MGA and EGA algorithms on large-scale linear systems of ordinary differential equations. J. Comput. Sci. 2022, 62, 101719. [Google Scholar] [CrossRef]
  14. Bentbib, A.; Jbilou, K.; Sadek, E.M. On some Krylov subspace based methods for large-scale nonsymmetric algebraic Riccati problems. Comput. Math. Appl. 2015, 70, 2555–2565. [Google Scholar] [CrossRef]
  15. Heyouni, M. Extended Arnoldi methods for large low-rank Sylvester matrix equations. Appl. Numer. Math. 2010, 60, 1171–1182. [Google Scholar] [CrossRef]
  16. Hached, M.; Jbilou, K. Computational Krylov-based methods for large-scale differential Sylvester matrix problems. Numer. Linear Algebra Appl. 2018, 25, e2187. [Google Scholar] [CrossRef]
  17. Druskin, V.; Knizhnerman, L. Extended Krylov subspaces: Approximation of the matrix square root and related functions. SIAM J. Matrix Anal. Appl. 1998, 19, 755–771. [Google Scholar] [CrossRef]
  18. Rosenbrock, H.H. Some general implicit processes for the numerical solution of differential equations. Comput. J. 1963, 5, 329–330. [Google Scholar] [CrossRef]
  19. Butcher, J.C. Numerical Methods for Ordinary Differential Equations; John Wiley & Sons: Hoboken, NJ, USA, 2016. [Google Scholar]
  20. Penzl, T. LYAPACK A MATLAB Toolbox for Large Lyapunov and Riccati Equations, Model Reduction Problems, and Linear-Quadratic Optimal Control Problems. Available online: http://www.tu-chemintz.de/sfb393/lyapack (accessed on 29 June 2022).
  21. Davis, T. The University of Florida Sparse Matrix Collection, NA Digest, Vol. 97, No. 23. 7 June 1997. Available online: http://www.cise.ufl.edu/research/sparse/matrices (accessed on 10 June 2016).
Figure 1. EBA-BD F method and EBA-Ros method.
Figure 1. EBA-BD F method and EBA-Ros method.
Mca 27 00069 g001
Figure 2. Residual norm vs number m of extended block Arnoldi iterations.
Figure 2. Residual norm vs number m of extended block Arnoldi iterations.
Mca 27 00069 g002
Figure 3. Residual norm vs number m of extended block Arnoldi iterations.
Figure 3. Residual norm vs number m of extended block Arnoldi iterations.
Mca 27 00069 g003
Table 1. Coefficients of the p-step BDF method with p 2 .
Table 1. Coefficients of the p-step BDF method with p 2 .
p β α 0 α 1
BDF11110
BDF22 2 / 3 4 / 3 1 / 3
Table 2. Results for Example 1.
Table 2. Results for Example 1.
Test ProblemMethodIterationsResidual NormTimes (s)
n = 8100 , p = 4900 , r = 2 , h = 0.3 EBA-BDF52.03 × 10 11 0.497
EBA-Ros5 1.40 × 10 11 0.854
n = 10,000, p = 4900 , r = 3 , h = 0.2 EBA-BDF5 4.38 × 10 12 1.025
EBA-Ros5 5.77 × 10 13 1.53
n = 40,000, p = 12,100, r = 4 , h = 0.1 EBA-BDF5 8.27 × 10 12 6.80
EBA-Ros5 8.70 × 10 13 6.85
Table 3. Runtimes in seconds and the residual norms.
Table 3. Runtimes in seconds and the residual norms.
Test ProblemMethodCPU TimeIterations R m ( T f ) F
n = p = 2500 r = 2 and h = 0.2 EBA-BDF 11.01 s16 4.57 × 10 9
A = a d d 32 and B = p d e 2961 EBA-Ros 11.49 s19 6.59 × 10 9
n = 12,100, p = 4900 EBA-BDF 10.77 s4 7.75 × 10 9
p = 8100 , h = 0.1 and r = 3 EBA-Ros 10.78 s5 6.31 × 10 10
B = t h e r m a l , n = 10,000, p = 3456 and r = 3 , h = 0.2 EBA-BDF 6.20 s7 5.87 × 10 9
A = fdm ( x + 10 y 2 , 2 x 2 + y 2 , y 2 x 2 ) EBA-Ros 6.27 s8 2.18 × 10 9
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Sadek, L.; Sadek, E.M.; Talibi Alaoui, H. On Some Numerical Methods for Solving Large Differential Nonsymmetric Stein Matrix Equations. Math. Comput. Appl. 2022, 27, 69. https://doi.org/10.3390/mca27040069

AMA Style

Sadek L, Sadek EM, Talibi Alaoui H. On Some Numerical Methods for Solving Large Differential Nonsymmetric Stein Matrix Equations. Mathematical and Computational Applications. 2022; 27(4):69. https://doi.org/10.3390/mca27040069

Chicago/Turabian Style

Sadek, Lakhlifa, El Mostafa Sadek, and Hamad Talibi Alaoui. 2022. "On Some Numerical Methods for Solving Large Differential Nonsymmetric Stein Matrix Equations" Mathematical and Computational Applications 27, no. 4: 69. https://doi.org/10.3390/mca27040069

APA Style

Sadek, L., Sadek, E. M., & Talibi Alaoui, H. (2022). On Some Numerical Methods for Solving Large Differential Nonsymmetric Stein Matrix Equations. Mathematical and Computational Applications, 27(4), 69. https://doi.org/10.3390/mca27040069

Article Metrics

Back to TopTop