Next Article in Journal
Approximate Closed-Form Solutions for the Maxwell-Bloch Equations via the Optimal Homotopy Asymptotic Method
Next Article in Special Issue
On Impulsive Implicit ψ-Caputo Hybrid Fractional Differential Equations with Retardation and Anticipation
Previous Article in Journal
Applied Geospatial Bayesian Modeling in the Big Data Era: Challenges and Solutions
Previous Article in Special Issue
New Sampling Expansion Related to Derivatives in Quaternion Fourier Transform Domain
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Computational Scheme for the First-Order Linear Integro-Differential Equations Based on the Shifted Legendre Spectral Collocation Method

Department of Mathematics and Statistics, Beijing Institute of Technology, Beijing 100081, China
*
Author to whom correspondence should be addressed.
Mathematics 2022, 10(21), 4117; https://doi.org/10.3390/math10214117
Submission received: 10 October 2022 / Revised: 1 November 2022 / Accepted: 1 November 2022 / Published: 4 November 2022
(This article belongs to the Special Issue Theory and Applications of Fractional Equations and Calculus)

Abstract

:
First-order linear Integro-Differential Equations (IDEs) has a major importance in modeling of some phenomena in sciences and engineering. The numerical solution for the first-order linear IDEs is usually obtained by the finite-differences methods. However, the convergence rate of the finite-differences method is limited by the order of the differences in L 1 space. Therefore, how to design a computational scheme for the first-order linear IDEs with computational efficiency becomes an urgent problem to be solved. To this end, a polynomial approximation scheme based on the shifted Legendre spectral collocation method is proposed in this paper. First, we transform the first-order linear IDEs into an Cauchy problem for consideration. Second, by decomposing the system operator, we rewrite the Cauchy problem into a more general form for approximating. Then, by using the shifted Legendre spectral collocation method, we construct a computational scheme and write it into an abstract version. The convergence of the scheme is proven in the sense of L 1 -norm by employing Trotter-Kato theorem. At the end of this paper, we summarize the usage of the scheme into an algorithm and present some numerical examples to show the applications of the algorithm.

1. Introduction

Integro-Differential Equations (IDEs) play an important role in applications of engineering, mechanics, physics and economics [1,2,3]. Particularly, First-order linear IDEs has a major importance in modeling of some phenomena in sciences and engineering, such as age-structured population dynamics, queuing theory, repairable systems reliability models and etc. [4,5,6,7].
A system whose state at time t is usually described by a vector-valued random variable X ( t ) = ( X 1 ( t ) , X 2 ( t ) , , X n ( t ) ) . For example, in the repairable system, X ( t ) may be a one-dimensional variable taking the value 1 or 0 from the functional state or the failed state, respectively. Alternately, in queueing system, X ( t ) may be a vector of the number of people in servers. If the sojourn time of the system in states i, Y i ( t ) , i = 1 , 2 , , n follows an exponential distribution, the probability distribution of the random variable X ( t ) can be investigated in the framework of Markov processes. By denoting P ( t ) = ( P 1 ( t ) , P 2 ( t ) , , P n ( t ) ) the corresponding probability of the system state at time t, the mathematical models of the system can be obtained according to the Chapman–Kolmogorov forward equations and described as a couple of ordinary differential equations. However, by using the exponential distribution to describe the sojourn time Y ( t ) is not suitable for the general cases. Thus, from the practical point of view, some non-exponential distributions, such as Weibull distribution, Erlang distribution and normal distribution, should be taken account into random process { X ( t ) } t 0 , which yields that { X ( t ) } t 0 is non-Markovian. To overcome this challenge, D.R. Cox [8] first put forward the supplementary variable technique and established the M/G/1 queueing model. After that, the supplementary variable technique was used and widely applied in other fields. Based on this technique, a non-Markovian process { X ( t ) } t 0 in continuous time is made Markovian { X ( t ) , Y ( t ) } t 0 by inclusion of one supplementary variable, and thus, the corresponding mathematical models are described by first-order linear IDEs. For example, Li [9] proposed a warm standby repairable system with priority in use and make a reliability analysis. Huo and Xu [10] studied the software rejuvenation system with periodic impulse. Zheng [11] studied the steady-state probability and reliability of a repairable system with three unites. The interested reader may for instance consult [9,10,11,12] and the references within.
As a result, the dynamic behavior of the queuing system, repairable system, the reliability system etc., can be studied by considering the solution of the corresponding first-order linear IDEs. By denoting λ j , i and b j , i ( · ) , j i , respectively, the transtion rate and the transtion rate function from state j to state i, the general form of the corresponding first-order linear IDEs can be expressed as
d P i ( t ) d t = λ i , i P i ( t ) + j = 1 , j i n λ j , i P j ( t ) + j = n + 1 n + m 0 b j , i ( x ) p j ( t , x ) d x , i = 1 , 2 , , n : p i ( t , x ) x + p i ( t , x ) t = ( λ i , i + b i , i ( x ) ) p i ( t , x ) + j = n + 1 , j i n + m λ j , i p j ( t , x ) , i = n + 1 , , n + m
with boundary conditions
p i ( t , 0 ) = j = 1 n λ j , i P j ( t ) + j = n + 1 , j i n + m 0 b j , i ( x ) p j ( t , x ) d x , i = n + 1 , , n + m .
where
λ i , i = j = 1 , j i n λ j , i , i = 1 , 2 , , n λ i , i = j = n + 1 , j i n + m λ j , i , b i , i ( x ) = j = n + 1 , j i n + m b j , i ( x ) , i = n + 1 , , n + m .
P i ( t ) , i = 1 , 2 , , n denotes the probability of system in state i at time t, p i ( t , x ) , i = n + 1 , , n + m denotes the probability of system in state i with elapsed time x at time t. x may be the age in population dynamics, the elapsed time for serving of server in queueing system or the elapsed time for repairing of component in repairable system, which are assumed to follow a general distribution with transition rate b i , j ( x ) satisfying
b i , j ( x ) 0 , x [ 0 , ) , sup x [ 0 , ) b i , j ( x ) < , 0 M b i , j ( x ) d x < , M [ 0 , ) , 0 b i , j ( x ) d x = . i , j = n + 1 , , n + m , i j .
However, from the view point of practical, the repair time of component, the life span of population or the service time of server can not be extended to infinity. For example, Yuan et al. [13] studied an optimal repair-replacement policy for a cold standby system. Zong et al. [14] considered a deteriorating system with increasing repair times and derived an optimal replacement policy. Therefore, it is reasonable to assume the elapsed time x for repairing belong to a finite interval [ 0 , T ) , where T denotes the maximum elapsed time for repairing of the component. That is, if the component cannot be repaired within time T, the component will be regarded as non-repairable and it will be replaced at time T immediately.
Thus, if we assume that the elapsed repair time of the component x follows the general distribution with transition rate b i , j ( x ) , i , j = n + 1 , , n + m and assume that the elapsed time x is in [ 0 , T ) , then we should redefine the transition rate and denote it as β i , j ( x ) , i , j = n + 1 , , n + m , which satisfying
0 x β i , j ( s ) d s = 0 x b i , j ( s ) d s , x < T , x = T .
Then, the corresponding first-order linear IDEs are transformed as
d P i ( t ) d t = λ i , i P i ( t ) + j = 1 , j i n λ j , i P j ( t ) + j = n + 1 n + m 0 T β j , i ( x ) p j ( t , x ) d x , i = 1 , 2 , , n : p i ( t , x ) x + p i ( t , x ) t = ( λ i , i + β i , i ( x ) ) p i ( t , x ) + j = n + 1 , j i n + m λ j , i p j ( t , x ) , i = n + 1 , , n + m
with boundary conditions
p i ( t , 0 ) = j = 1 n λ j , i P j ( t ) + j = n + 1 , j i n + m 0 T β j , i ( x ) p j ( t , x ) d x , i = n + 1 , , n + m .
To investigate the well-posedness and stability of the system, the first-order linear IDEs is usually transformed into an abstract Cauchy problem by introducing a suitable state space, system operators and their domains. Readers can refer to Gupur [15] for viewing the functional analysis methods on studying the reliability models. Therefore, we define the system space of (4) and (5) as X = R n × ( L 1 [ 0 , T ] ) m with norm · X ,
P X = i = 1 n | P i | + i = n + 1 m p i ( x ) L 1 [ 0 , T ] , P = ( P 1 , P 2 , , P n , p n + 1 ( x ) , , p n + m ( x ) ) X .
Then, the system (4) and (5) can be translated into an abstract Cauchy problem and written as
d d t P ( t , · ) = A P ( t , · ) , P ( 0 , · ) = P 0 ,
where operator A is a differential operator with domain
D ( A ) = ( p 1 , p 2 , , p n , p n + 1 ( x ) , p n + 2 ( x ) , , p n + m ( x ) ) τ X | p i ( x ) are absolutely continuous functions , d p i ( x ) d x L 1 [ 0 , T ] , and   ( p n + 1 ( 0 ) , p n + 2 ( 0 ) , , p n + m ( 0 ) ) = B P , i = n + 1 , n + 2 , , n + m . .
and B is a boundary operator from X to R m . By using C 0 -semigroup theory [16], it can be verified that system (6) is well posed, which ensures the existence of the time-dependent solution. The interested reader may for instance consult [17,18].
The system (6) usually are difficult to be solved analytically, which yields that much work focus on investigating the steady state solution [19,20,21] based on the qualitative analysis. Therefore, establishing an numerical computational framework for system (6) becomes an important issue. More and more methods were arosed to study the numerical solution of the IDEs, such as the successive approximations, spectral Galerkin method and collocation methods, which can be found in [22,23,24] and the references therein. Recently, Dzhumabaev [25,26,27] devoted to the linear boundary value problem for the Fredholm IDE. These paper not only established the necessary and sufficient conditions for solvability of the IDEs, but also provided algorithms with good computational accuracy in application. To solve Fredholm integral equations on the half-line, Rahmoune [28] proposed a spectral collocation method based on the scaled Laguerre functions. The proposed method provided a good efficiency for smooth solutions decaying at infinity. To study the Volterra IDEs, Al-Ahmad et al. [29], based on the differential transform method (DTM), proposed a modified differential transform method, which can greatly improve the convergence rate of DTM’s truncated series solution with true analytic solution. Xu et al. [30] proposed a Half-Sweep Successive Over-Relaxation method to investigate the first-order linear Fredholm IDEs. Moreover, they verified that the proposed method is superior to the Full-Sweep GaussSeidel methods and Full-Sweep Successive Over-Relaxation methods. For the nonlinear problems, Dawood et al. [31] proposed a Laplace Discrete Adomian Decomposition Method to solve nonlinear Volterra-Fredholm IDEs. The proposed method was simple to execute and can be effectively used to overcome the analytical approaches in solving nonlinear Fredholm IDEs. Combining with the method of cutting functions, squaring method, the Cauchy-Bunyakovsky inequality, and the integral inequality method, Iskandarov [32] studied a weakly first order nonlinear implicit Volterra IDE on the semiaxis. All in all, the methods proposed in reference [29,30,31,32] are synthetical and possess good computational accuracy for the solution of IDEs.
However the convergence of these algorithms based on polynomials or collocation methods is always worked in the framework of L 2 space, which yields that some approximation methods are not suitable in L 1 space, since it is non-reflexivity. Therefore, the finite-differences methods became the main tool for deriving the numerical solutions of system (6). Xu [33] presented a simple finite-differences scheme and established the convergence of the scheme by employing Trotter-Kato Theorem for a reliability model consisting of two machines separated by a finite storage buffer. Applications for the finite-differences methods cover researches including software rejuvenation system [10], repairable systems with preventive maintenance [34], optimal maintenance design for a Simple Reparable System [35] and so on. However, the numerical results obtained by the finite-differences methods are in general very poor, since that finite-differences is a local method, which is to compute the derivative at a given point by using information only pertaining to a small neighborhood of the point [22]. The convergence rate of the finite-differences method is also limited by the order of the method. Therefore, how to design a computational scheme for the first-order linear IDEs with computational efficiency becomes an urgent problem to be solved.
To this end, a computational scheme based on the shifted Legendre spectral collocation method as well as the corresponding approximation algorithm is proposed. The remainder of this paper is organized as follows. Section 2 designs an approximation scheme based on the basic theoretical framework for the Legendre collocation method. Section 3 demonstrates the convergence of the scheme in the sense of L 1 -norm by applying Trotter-Kato Theorem. Section 4 proposes an computational algorithm for the scheme and performs some numerical examples. Section 5 concludes the paper.

2. Approximation Scheme Design

In this section, an approximation scheme to the abstract version (6) is designed. Firstly, by decomposing the system operator of (6), we rewrite the system into an more general form for approximation scheme construction. Secondly, some important results in approximation theory are reviewed as the basic theoretical framework for scheme construction. Finally, a polynomial approximation scheme for (6) is proposed.

2.1. System Formulating

In this subsection, we rewrite the system (6) into an more general form. We decompose the system operator A according to the system space. For any P ( t , · ) = ( P 1 ( t ) , P 2 ( t , · ) ) X , where P 1 ( t ) R n and P 2 ( t , · ) ( L 1 [ 0 , T ] ) m , we define operators A i , i = 1 , 2 , 3 , 4 , and L i X , i = 1 , 2 , and rewrite the abstract version (6) as
d d t P ( t , · ) = A P ( t , · ) = A 1 A 2 L 1 ( L 1 ) 0 m × n A 3 + L 2 ( L 1 ) · P 1 ( t ) P 2 ( t , · ) , P ( 0 , · ) = P 0 .
with the boundary conditions P 2 ( t , 0 ) = A 4 P 1 ( t ) , where 0 m × n denotes the ( m × n ) - dimensional zero matrix, A 1 : R n R n and A 4 : R n R m are linear continuous operators. Operator A 2 : ( L 1 [ 0 , T ] ) n R n and operator A 3 : ( L 1 [ 0 , T ] ) m ( L 1 [ 0 , T ] ) m are denoted as the integral operator and the differential operator, respectively, which is that for any Y = y 1 ( x ) , y 2 ( x ) , , y n ( x ) τ ( L 1 [ 0 , T ] ) n ,
A 2 Y = 0 T | y 1 ( x ) | d x , 0 T | y 2 ( x ) | d x , , 0 T | y n ( x ) | d x τ , A 3 Y = d d x y 1 ( x ) , d d x y 2 ( x ) , , d d x y n ( x ) τ .
Moreover, operator L 1 X : X m X n , and L 2 X : X m X m are linear operators. Without loss of generality, we denote L 1 ( L 1 ) : = ( u k j ( x ) ) k = 1 , 2 , , m ; j = 1 , 2 , , n and L 2 ( L 1 ) : = ( v k j ( x ) ) k , j = 1 , 2 , , m , where u k j ( x ) and v k j ( x ) are corresponding to the transient rate functions of the system and defined according to the coupling relationship among the states.
We should note that system (7) is equivalent to system (6) for the situation that without integral boundary conditions. Therefore, linear continuous operator A 4 is a constant matrix. The situation for system (6) with integral boundary conditions also can be investigated by the same way. Thus, we omit it. On the other hand, we claim that the rewriting form (7) for the abstract version (6) is unique. The boundary operators A 4 is straightforward. Then, once the integral operator A 2 and the differential operator A 3 of the system are determined, it is natural to obtain the operator A i and L i ( L 1 ) , i = 1 , 2 . A specific example will be shown in Section 4.

2.2. Approximation Theory

In this subsection, we review the theory of polynomial approximations to continuous functions. We begin with the collocation method, which is a widely applied and powerful technique in the construction of numerical methods for differential or integral equations. A collocation method is based on the idea of approximating the exact solution of a given equation with a suitable approximant belonging to a finite dimensional space, usually a piecewise algebraic polynomial, which exactly satisfies the equation on a certain finite subset of the integration interval.
Legendre spectral collocation method (pseudospectral method) [24], which is to seek an approximate solution f ( N ) ( x ) to f ( x ) in the form
f ( N ) ( x ) = k = 1 N f ( x k ) · h k ( N ) ( x ) .
where { x j } j = 1 N are the Legendre collocation points, which corresponding to the zeros of the Legendre polynomials L N ( x ) , x ( 1 , 1 ) . The Legendre polynomials are important special cases of the Jacobi polynomials and they are mutually orthogonal with respect to the uniform weight function w ( x ) = 1 . Moreover, { h k ( N ) ( x ) } k = 1 N in (8) are the Lagrange basis polynomials associated with the collocation points { x j } j = 1 N and
h k ( N ) ( x ) = i = 1 , i k N x x i x k x i , 1 k N ,
with the alternate expression
h k ( N ) ( x ) = L N ( x ) x L N ( x k ) ( x x k ) , 1 k N .
Furthermore, we introduce the shifted Legendre polynomial L N , T ( x ) to generalize the Legendra spectral collocation method to any bounded interval for application. That is
L N , T ( x ) = L N ( 2 x T 1 ) , x ( 0 , T ) ,
Readers can review the other properties and applications of Legendre polynomials as well as its shifted form in reference [24,36]. As a result, the polynomials approximation based on the shifted Legendra collocation method to function f ( x ) in bounded interval ( 0 , T ) can be defined as
f ( N ) ( x ) = k = 1 N f ( x ¯ k ) · h ¯ k ( N ) ( x ) , x ( 0 , T ) .
where x ¯ j = T 2 ( x j + 1 ) , j = 1 , 2 , , N call the shifted Legendre collocation points and { h ¯ k ( N ) ( x ) } k = 1 N are the Lagrange basis polynomials associated with { x ¯ j } j = 1 N .
However, in many practical applications, we also should consider the integral of the approximate solution, which leads to the consideration between orthogonal polynomials and Gauss-type quadrature, which is to seek the best numerical approximation of an integral by selecting optimal nodes at which the integrand is evaluated. However, it can be verified that all the nodes of the shifted Legendre polynomials lie in the interior of the interval ( 0 , T ) . This makes it difficult to impose boundary conditions when dealing with the problems as in system (6). To this end, we consider the shifted Gauss-Radau quadratures which includes the endpoints 0 as node. We introduce the following important Gauss-Radau quadrature Theorem, which established in [24].
Theorem 1
([24]). Let x 0 = a and { x j } j = 1 N be the zero of q N defined as
q N ( x ) = p N + 1 + α N p N ( x ) x a with α N = p N + 1 ( a ) p N ( a )
where p N is the N order orthogonal polynomial. { q N : N 0 } defines a sequence of polynomials orthogonal with respect to the weight function ω ¯ ( x ) : = ω ( x ) ( x a ) , and the leading coefficient of q N is k N + 1 . Then there exists a unique set of quadrature weights { ω j } j = 0 N , which defined as
ω j = a b h j ( x ) ω ( x ) d x , 0 j N ,
where { h j ( x ) } j = 0 N are the Lagrange basis polynomials associated with { x j } j = 1 N , such that
a b p ( x ) ω ( x ) d x = j = 0 N p ( x j ) ω j , p P 2 N .
Moreover, the quadrature weights are all positive and can be expressed as
ω 0 = 1 q N ( a ) a b q N ( x ) ω ( x ) d x , ω j = 1 x j a k N + 1 k N q N 1 ( x ) ω ¯ 2 q N 1 ( x j ) q N ( x j ) , 1 j N .
By using Theorem 1, we propose the following important theorm.
Theorem 2.
Let { x j , ω j } j = 0 N be the set of Legendre-Gauss Radau quadrature nodes and weights. Then, we have { x j } j = 0 N are the zeros of L N ( x ) + L N + 1 ( x ) and
ω j = 1 ( N + 1 ) 2 1 x j [ L N ( x j ) ] 2 , 0 j N .
Moreover, let { x ^ j , ω ^ j } j = 0 N be the set of shifted Legendre-Gauss Radau quadrature nodes and weights. Then, we have
x ^ j = T 2 ( x j + 1 ) , ω ^ j = T 2 ω j , 0 j N .
with the above nodes and weights, we have
0 T p ( x ) d x = j = 0 N p ( x ^ j ) · ω ^ j , p ( x ) P 2 N 1 .
where P 2 N 1 denotes the space of polynomial with degree at most 2 N 1 .
Proof. 
According to Theorem 1, the Legendre-Gauss Radau quadrature nodes and weights { x j , ω j } j = 0 N can be derived with the properties of Legendre polynomials. Moreover, the corresponding shifted nodes and weights { x ^ j , ω ^ j } j = 0 N in Equation (12) has been proved in [36]. Therefore, Theorem 2 is straightforward. □
However, Theorem 2 is strictly true only for p ( x ) P 2 N 1 , which leads to the consideration of the convergence analysis for the case of general continuous functions. We define the interpolation operators I N : C 0 [ 0 , T ] P N as
I N f ( x ) = f ^ ( N ) ( x ) : = k = 0 N f ( x ^ k ) · h ^ k ( N ) ( x ) .
where h ^ k ( N ) ( x ) denotes the Lagrange basis polynomials associated with the shifted Legendre-Gauss Radau quadrature nodes { x ^ j } j = 0 N . Then, we consider the Legendre-Gauss Radau quadrature formulas for the continuous functions with increasing number of nodes. We proposed the following theorem which as a corollary of a convergence result from [22].
Theorem 3.
Let f ( x ) is Riemann integrable in [ 0 , T ] . Then, we have
lim N j = 0 N f ( x ^ j ) · ω ^ j = lim N 0 T I N f ( x ) d x = 0 T f ( x ) d x .
where operator I N is defined in (14), { x ^ j } j = 0 N and { ω ^ j } j = 0 N are the shifted Legendre-Gauss Radau quadrature nodes and weights respectively.
Proof. 
According to [22], we can obtain the convergence theorem of the Jacobi Gaussian formulas to the function f ( x ) which satisfies that f ( x ) · w j b ( x ) is Riemann integrable in ( 1 , 1 ) , where w j b ( x ) is the Jacobi weights function. Since Legendre polynomials are the spacial cases of Jacobi polynomials with w j b ( x ) = 1 , the convergence theorem for the Legendre Gaussian formulas can be derived, which yields the establishment for the shifted Legendre-Gauss case. □
Hereafter, we consider the differentiation of Equation (14). Differentiating Equation (14) m times leads to
d m d x m f ^ ( N ) ( x ^ k ) = j = 0 N d ^ k j ( m ) f ^ ( N ) ( x ^ j ) ,
The matrix D ^ ( m ) = ( d ^ k j ( m ) ) k , j = 1 , 2 , , N is called the differentiation matrix of order m relative to { x ^ j } j = 1 N with the entries
d ^ k j ( m ) = d m d x m h j ( N ) ( x ^ k ) = N ( N + 2 ) 2 T , k = j = 0 , 2 T · x ^ k 1 x ^ k 2 + ( N + 1 ) L N ( x ^ k ) ( 1 x ^ k 2 ) Q N ( x ^ k ) , 1 k = j N , Q N ( x ^ k ) Q N ( x ^ j ) · 1 x ^ k x ^ j , k j ,
where Q N ( x ) = L N + 1 ( x ) + L N . Moreover, if we denote by f ( m ) the vector whose components are the values of d m d x m f ^ ( N ) ( x ^ j ) , j = 1 , 2 , , N , it follows from (15) and (16) that
f ( m ) = D ^ ( m ) f ( 0 ) .
Equation (17) implies that, once the differentiation matrices D ^ ( m ) are precomputed, the derivative values can be evaluated. Therefore, it is very convenient for solving problems with variable coefficients and nonlinear problems in the physical space with such collocation method.

2.3. Approximation Scheme Design

In this subsection, we construct an approximation scheme for system (7) with the shifted Legendre collocation method which shown in Section 2.2.
Let { x ^ j , ω ^ j } j = 0 N be the set of shifted Legendre-Gauss Radau quadrature nodes and weights defined in Equation (12). Let R n = P 1 = ( p 1 , p 2 , , p n ) τ R n : P 1 1 = i = 1 n | p i | and R N = P 2 = ( p 1 , p 2 , , p N ) τ R N : P 2 R N = k = 1 N ω ^ k · p k . Then, we define the approximating space as
X ( N ) = P ( N ) = ( P 1 , P 2 ( N ) ) R n × ( R N ) m : P ( N ) N = P 1 1 + P 2 ( N ) + ω ^ 0 · A 4 P 1
where operator A 4 is the boundary operator in system (7). Consequently, based on the shifted Legendre collocation method, the approximation system for system (7) can be constructed as
d d t P ( N ) ( t ) = A ( N ) P ( N ) ( t ) = A 1 + ω ^ 0 · L 1 ( R ) A 4 Ψ 1 ( N ) L 1 ( R N ) D ^ 0 A 4 D ^ 1 I N + L 2 ( R N ) · P 1 ( t ) P 2 ( N ) ( t ) P ( N ) ( 0 ) = P 0 ( N ) .
where ⊗ is the Kronecker product. The Kronecker product of matrices A = ( a i j ) R m × n and B = ( b i j ) R p × q is an m p × n q matrix defined as A B = ( ( a i j · B ) ) . I N is a n-dimensional identity matrix.
Operator A ( N ) is a linear continuous operator form X ( N ) to X ( N ) . It can be verified that A ( N ) is the infinitesimal generator of a C 0 semigroup of contraction T N ( t ) on X ( N ) . Operator A i , i = 1 , 4 are defined as in (7). Operator D ^ 0 = ( d ^ k j ) j = 0 , k = 1 , 2 , , N and D ^ 1 = ( d ^ k j ) k , j = 1 , 2 , , N are defined according to the differentiation matrix in (17). Moreover, operator L 1 ( R ) = ( u k j ( 0 ) ) k = 1 , 2 , , n , j = 1 , 2 , , m . Operator Ψ 1 ( N ) : ( R N ) n R n and L 1 ( R N ) : ( R N ) m ( R N ) n are defined to approximate the intergal items and
Ψ 1 ( N ) = diag ( ω , ω , , ω ) , ω = ( ω 1 , ω 2 , , ω N ) , L 1 ( R N ) = ( u k j ) k = 1 , 2 , , m ; j = 1 , 2 , , n , u k j = diag ( u k j ( x 1 ) , u k j ( x 2 ) , , u k j ( x N ) )
Operator L 2 ( R N ) : ( R N ) m ( R N ) m is defined as
L 2 ( R N ) = ( v k j ) k , j = 1 , 2 , , m , v k j = diag ( v k j ( x 1 ) , v k j ( x 2 ) , , v k j ( x N ) )

3. Approximation Scheme Convergence

In this section, the convergence for the scheme is investigated. Firstly, we review the Trotter-Kato theorem, which is useful for studying convergence of numerical approximations of solutions to partial differential equations. Secondly, by employing Trotter-Kato theorem, we demonstrate the convergence for the abstract version (18) to (7).

3.1. The Trotter-Kato Theorem

Let X and X N be Banach space with norm · , · N , N = 1 , 2 , , respectively, and assume that T ( t ) is a C 0 -semigroup on X . A G ( M , δ , X ) , M 1 , δ R , means that A is the infinitesimal generator of a C 0 -semigroup T ( t ) , t 0 , satisfying T M e δ t . Also assume that, for each n = 1 , 2 , , there exists bounded linear operators P N : X X N and E N : X N X satisfying.
(A1)
P N M 1 , E N M 2 , where M 1 and M 2 are independent of N;
(A2)
E N P N x x 0 as N , for all x X ;
(A3)
P N E N = I N , where I N is the identity operator on X N .
Theorem 4.
(Trotter-Kato). Assume that (A1) and (A3) are satisfied. Let A resp. A N be in G ( M , δ , X ) resp. in G ( M , δ , X N ) and let T ( t ) and T N ( t ) be the semigroup generated by A and A N on X and X N respectively. Then the following statements are equivalent:
(a)
There exists a λ 0 ρ ( A ) n = 1 ρ ( A N ) such that, for all x X ,
E N ( λ 0 I N A N ) 1 P N x ( λ 0 I A ) 1 x 0 as N ,
(b)
for every x X and t 0 ,
E N T N ( t ) P N x T ( t ) x 0 as N ,
uniformly on bounded t-intervals.
Theorem 4 is the functional analysis form of the Lax equivalent theorem [37]. In order to apply the Theorem 4 to a special problem, it is usually difficult to verify the resolvent convergence. To overcome this difficulty, reference [38] pointed out that condition (a) can be verified by the following theorem.
Theorem 5.
Let the assumptions of Theorem be satisfied. Then statement (a) of the Theroem 4 is equivalent to (A2) and the following two statements:
(C1)
There exists a subset D dom A such that D ¯ = X and ( λ 0 I A ) D ¯ = X for a λ 0 > δ ;
(C2)
For all x D there exists a sequence ( x ¯ N ) N N with ( x ¯ N ) N N dom A N such that
lim N E N x ¯ N = x , and lim N E N A N x ¯ N = A x
Condition (C1) is the statement that subset D is the core of the operator A. Thus, condition (C2) is the statement that A N converge pointwise on a core of A, see [16].

3.2. Convergence for the Scheme

In this subsection, we employ the Trotter-Kato theorem to obtain the convergence. Let X ( N ) be defined as in Section 2.3. Then, for any
y = ( p 1 , p 2 , , p n , ( p n + 1 ( x ^ 1 ) , , p n + 1 ( x ^ N ) , , p n + m ( x ^ 1 ) , , p n + m ( x ^ N ) ) τ X ( N ) , ϕ = ( ( p 1 , p 2 , , p n ) , ( p n + 1 ( x ) , , p n + m ( x ) ) ) τ X ,
let P N and E N be defined as
E N ( y ) = ( p 1 , p 2 , , p n , P n + 1 ( x ) , P n + 2 ( x ) , P n + m ( x ) ) τ ,
P N ( ϕ ) = ( p 1 , p 2 , , p n , p n + 1 ( x ^ 1 ) , , p n + 1 ( x ^ N ) , , p n + m ( x ^ 1 ) , , p n + m ( x ^ N ) ) τ ,
where
P n + j ( x ) = h ^ 0 ( N ) ( x ) · p ¯ j + k = 1 N h ^ k ( N ) ( x ) · p n + j ( x ^ k ) j = 1 , 2 , , m
p ¯ j is the jth items of A 4 ( p 1 , p 2 , , p n ) τ , and h ^ k ( N ) ( x ) , k = 1 , 2 , , N are the Lagrange basis polynomials associated with the shifted Legendre-Gauss Radau quadrature nodes { x ^ j } j = 0 N .
To derived the convergence for the scheme by using the Trotter-Kato theorem, we only need to verify that E N and P N satisfy assumption (A1)–(A3), and prove that both conditions (C1) and (C2) in Theorem 5 hold. To this end, we prove the following two propositions.
Proposition 1.
Let E N and P N be defined as (19) and (20), then,
(A1)
P N M 1 , E N M 2 , where M 1 and M 2 are independent of N;
(A2)
E N P N x x 0 as N , for all x X ;
(A3)
P N E N = I N , where I N is the n-dimensional identity operator on X ( N ) .
Proof. 
For any y = ( p 1 , p 2 , , p n , p n + 1 ( x ^ 1 ) , , p n + 1 ( x ^ N ) , , ( p n + m ( x ^ 1 ) , , p n + m ( x ^ N ) ) ) τ X ( N ) ,
E N y X = i = 1 n | p i | + j = 1 m 0 T | h ^ 0 ( N ) ( x ) · p ¯ j + k = 1 N h ^ k ( N ) ( x ) · p n + j ( x ^ k ) | d x i = 1 n | p i | + j = 1 m 0 T | h ^ 0 ( N ) ( x ) | d x · p ¯ j + k = 1 N 0 T | h ^ k ( N ) ( x ) | d x · p n + j ( x ^ k ) i = 1 n | p i | + j = 1 m ω ^ 0 · p ¯ j + k = 1 N ω ^ k · p n + j ( x ^ k ) = y N .
which yields E N 1 . Then, for any ϕ = ( p 1 , p 2 , , p n ) , ( p n + 1 ( x ) , , p n + m ( x ) ) τ X , since D ( A ) is densed in X , there exists a ϕ = ( ( p 1 , p 2 , , p n ) , ( p n + 1 ( x ) , , p n + m ( x ) ) ) τ D ( A ) such that ϕ ϕ X < ε . Moreover
P N ( ϕ ) N P N ( ϕ ϕ ) N + P N ( ϕ ) N = i = 1 n | p i p i | + j = 1 m k = 1 N ω ^ k · | p n + j ( x ^ k ) p n + j ( x ^ k ) | + j = 1 m ω ^ 0 · | p ¯ p n + j ( 0 ) | + i = 1 n | p i | + j = 1 m [ k = 0 N ω ^ k · p n + j ( x ^ k ) ] M ε + i = 1 n | p i | + k = 1 N 0 T I N p n + j ( x ) d x ϕ X
where x 0 = 0 , p ¯ j is the jth items of p ¯ = A 4 ( p 1 , p 2 , , p n ) τ and M = max { ω ^ k , k = 0 , 1 , , N } . This implies P N 1 . Therefore, (A1) is satisfied. After that, for any ϕ = ( ( p 1 , p 2 , , p n ) , ( p n + 1 ( x ) , , p n + m ( x ) ) ) τ X , if ϕ D ( A ) , then
E N P N ϕ ϕ X = j = 1 m 0 T | k = 0 N h ^ k ( N ) ( x ) · p n + j ( x ^ k ) p n + j ( x ) | d x = j = 1 m 0 T | I N p n + j ( x ) p n + j ( x ) | d x 0 as N .
Otherwise, there is a ϕ = ( p 1 , p 2 , , p n ) , ( p n + 1 ( x ) , , p n + m ( x ) ) τ D ( A ) such that ϕ ϕ X < ε , and
E N P N ϕ ϕ X E N P N ϕ E N P N ϕ X + E N P N ϕ ϕ X + ϕ ϕ X ( E N P N + 1 ) · ϕ ϕ X + E N P N ϕ ϕ X 0 as N ,
Thus, (A2) is satisfied. P N E N = I N is obvious. Therefore, (A3) follows. □
Secondly, we will prove that both (C1) and (C2) hold. To this end, we choose D = D ( A ) , which establishes condition (C1) in Theorem 5 with δ = 0 , the proof is rather straightforward and will therefore be omitted. On the other hand, we prove that condition (C2) also holds.
Proposition 2.
For all ϕ = ( ϕ 1 , ϕ 2 ) = ( ( p 1 , , p n ) , ( p n + 1 ( x ) , , p n + m ( x ) ) ) τ D ( A ) X , there exists a sequence ϕ ( N ) X ( N ) , which defined as ϕ ( N ) = ( ϕ 1 , ϕ 2 ( N ) ) =   ( ( p 1 , , p n ) , ( p n + 1 ( x ^ 1 ) , , p n + 1 ( x ^ N ) ) , , ( p n + m ( x ^ 1 ) , , p n + m ( x ^ N ) ) ) τ , such that
lim N E N ϕ ( N ) = ϕ , and lim N E N A ( N ) ϕ ( N ) = A ϕ
Proof. 
First, we consider E N ϕ ( N ) ϕ X that is
E N ϕ ( N ) ϕ X = j = 1 m 0 T | h ^ 0 ( N ) ( x ) · p ¯ j + k = 1 N h ^ k ( N ) ( x ) · p n + j ( x ^ k ) p n + j ( x ) | d x = j = 1 m 0 T | I N p n + j ( x ) p n + j ( x ) | d x 0 as N ,
Then, we investigate E N A ( N ) ϕ ( N ) A ϕ X , according to (7) and (18),
A ϕ = A 1 · ϕ 1 + A 2 L 1 ( L 1 ) ϕ 2 ( A 3 + L 2 ( L 1 ) ) ϕ 2 : = Φ 1 Φ 2
A ( N ) ϕ ( N ) = ( A 1 + ω ^ 0 · L 1 ( R ) A 4 ) · ϕ 1 + Ψ 1 ( N ) L 1 ( R N ) · ϕ 2 ( N ) D ¯ 0 A 4 · ϕ 1 + ( D ¯ I N + L 2 ( R N ) ) · ϕ 2 ( N ) : = Φ ¯ 1 Φ ¯ 2
Moreover,
Φ ¯ 1 = A 1 · ϕ 1 + ω ^ 0 · L 1 ( R ) A 4 · ϕ 1 + Ψ 1 ( N ) L 1 ( R N ) · ϕ 2 ( N ) = A 1 · ϕ 1 + ω ^ 0 · L 1 ( R ) · ϕ 2 ( N ) ( 0 ) + Ψ 1 ( N ) L 1 ( R N ) · ϕ 2 ( N ) = A 1 · ϕ 1 + A 2 I N ( n ) L 1 ( L 1 ) ϕ 2 . Φ ¯ 2 = D ¯ 0 A 4 · ϕ 1 D ¯ I N · ϕ 2 ( N ) + L 2 ( R N ) · ϕ 2 ( N ) = D ¯ 0 ϕ 2 ( 0 ) D ¯ I N · ϕ 2 ( N ) + L 2 ( R N ) · ϕ 2 ( N ) = Π ( N ) A 3 I N ( m ) ϕ 2 + L 2 ( R N ) · ϕ 2 ( N ) .
where ϕ 2 ( 0 ) = p n + 1 ( 0 ) , p n + 2 ( 0 ) , , p n + m ( 0 ) τ , operator I N ( n ) : ( L 1 [ 0 , T ] ) n ( P N ) n and operator Π ( N ) : ( P N ) n ( R N + 1 ) n , which satisfy that
I N ( n ) f = ( I N f 1 ( x ) , I N f 2 ( x ) , , I N f n ( x ) ) τ , f = f 1 ( x ) , , f n ( x ) τ ( L 1 [ 0 , T ] ) n . Π ( N ) g = ( g 1 ( x ^ 1 ) , , g 1 ( x ^ N ) ) , , ( g n ( x ^ 1 ) , , g n ( x ^ N ) ) τ , g = g 1 ( x ) , , g n ( x ) τ ( P N ) n .
Therefore,
E N A ( N ) ϕ ¯ ( N ) = E N ( Φ ¯ 1 , Φ ¯ 2 ( N ) ) τ = E N A 1 · ϕ 1 + A 2 I N ( n ) L 1 ( L 1 ) ϕ 2 Π ( N ) A 3 I N ( m ) ϕ 2 + L 2 ( R N ) · ϕ 2 ( N ) = A 1 · ϕ 1 + A 2 I N ( n ) L 1 ( L 1 ) ϕ 2 A 3 I N ( m ) ϕ 2 + I N ( m ) L 2 ( L 1 ) · ϕ 2 + R m
where R m is a m-dimensional column vector, and
R m = h ^ 0 ( N ) ( x ) · A 4 Φ ¯ 1 A 3 I N ( m ) ϕ 2 | x = 0 L 2 ( R ) · ϕ 2 ( 0 ) ,
with L 2 ( R ) = ( v k j ( 0 ) ) k , j = 1 , 2 , , n . It can be verified that
R m = h ^ 0 ( N ) ( x ) · A 4 A 1 · ϕ 1 + A 2 I N ( n ) L 1 ( L 1 ) ϕ 2 A 3 I N ( m ) ϕ 2 | x = 0 L 2 ( R ) · ϕ 2 ( 0 ) = h ^ 0 ( N ) ( x ) · A 4 A 1 · ϕ 1 + A 2 L 1 ( L 1 ) ϕ 2 + A 2 I N ( n ) L 1 ( L 1 ) ϕ 2 A 2 L 1 ( L 1 ) ϕ 2 A 3 I N ( m ) ϕ 2 | x = 0 L 2 ( R ) · ϕ 2 ( 0 ) = h ^ 0 ( N ) ( x ) · A 4 A 2 I N ( n ) L 1 ( L 1 ) ϕ 2 A 2 L 1 ( L 1 ) ϕ 2 + A 4 Φ 1 A 3 I N ( m ) ϕ 2 | x = 0 L 2 ( R ) · ϕ 2 ( 0 ) = h ^ 0 ( N ) ( x ) · A 4 A 2 I N ( n ) L 1 ( L 1 ) ϕ 2 A 2 L 1 ( L 1 ) ϕ 2 + A 3 ϕ 2 | x = 0 A 3 I N ( m ) ϕ 2 | x = 0 + L 2 ( L 1 ) ϕ 2 L 2 ( R ) · ϕ 2 ( 0 ) .
Thus, it is obvious that R m 0 as N . Consequently, we have
E N A ( N ) ϕ ¯ ( N ) A ϕ X = A 2 I N ( n ) L 1 ( L 1 ) ϕ 2 A 2 L 1 ( L 1 ) ϕ 2 1 + A 3 I N ( m ) ϕ 2 + I N ( m ) L 2 ( L 1 ) · ϕ 2 A 3 ϕ 2 L 2 ( L 1 ) ϕ 2 + R m i = 1 n 0 T | I N ( L 1 ( L 1 ) ϕ 2 ) i ( L 1 ( L 1 ) ϕ 2 ) i | d x + i = 1 m 0 T | A 3 I N p n + i ( x ) p n + i ( x ) | d x + i = 1 m 0 T | I N ( L 2 ( L 1 ) ϕ 2 ) i ( L 2 ( L 1 ) ϕ 2 ) i | d x + R m 0 as N .
where ( L 1 ( L 1 ) ϕ 2 ) i and ( L 2 ( L 1 ) ϕ 2 ) j denote the ith items of the n-dimensional column vector L 1 ( L 1 ) ϕ 2 and the jth items of the m-dimensional column vector L 2 ( L 1 ) ϕ 2 , respectively. □
Since the operator A generates a C 0 -semigroup T ( t ) of contraction in X , and operator A ( N ) generates a C 0 -semigroup T N ( t ) of contraction in X ( N ) , we can summarize the convergence theorem according to Propositions 1 and 2 and Theorem 5 that
Theorem 6.
Let E N and P N be defined as (19) and (20), A and A ( N ) be the system operators defined in system (7) and system (18). Let T ( t ) and T N ( t ) be the semigroup generated by A and A ( N ) on X and X ( N ) , respectively, then for any x X and t [ 0 , T ] , E N T N ( t ) P N x T ( t ) x X 0 as N uniformly on bounded t-intervals.
As a result, we have the following corollary.
Corollary 1.
Since the dynamical solution of (7) can be expressed as p * ( t ) = T ( t ) P 0 , here P 0 is the initial value of the system. Then, we have following relationship,
E N T N ( t ) P N P 0 T ( t ) P 0 X 0 as N ,
where T N ( t ) is the semigroup with generator A ( N ) , and T N ( t ) = e A ( N ) t .
Proof. 
Since A can generate C 0 -semigroup T N ( t ) , the dynamical solution could be expressed as p * ( t ) = T ( t ) P 0 , where P 0 is the initial value of the (7). Moreover, A ( N ) is bounded operator, which generates a uniformly continuous semigroup. Therefore, the conditions of Theorem 6 are satisfied, and
E N T N ( t ) P N P 0 T ( t ) P 0 X 0 as N ,
holds. Corollary 1 provides us with a useful tool to calculate the dynamical solution and instantaneous indexed of the model. □

4. Scheme Application

In this section, we summarize the usage for the approximation scheme in Section 2.3 into an algorithm and present some numerical experiments as applications to the algorithm.

4.1. Approximation Algorithm

In this subsection, we propose an approximation algorithm, which can be an useful tool in calculation for the dynamical solution or instantaneous indexes to the system. According to aforementioned analysis, an algorithm for the system which can be described as an Cauchy problem as (6) is proposed as follows.
  • Step 1: Inputting system parameters. Input the system parameters for the system operator A and boundary operator B . Input T as the upper boundary value of the transient rate function. Input N as the number of discrete nodes to approximation.
  • Step 2: Calculating polynomial parameters. Calculating the shifted Legendre-Gauss Radau quadrature nodes { x ^ j } j = 0 N weights { ω ^ j } j = 0 N according to Theorem 2 and Equation (12). Calculating the differention matrix D ( 1 ) with Equation (16).
  • Step 3: Decompositing the system operator. By definded the integral operator A 2 , differential operator A 3 and boundary operator A 4 , we decompose the system operator A into A i , i = 1 , 2 , 3 , and L i ( L 1 ) , i = 1 , 2 , which are defined as in Section 2.3, and
    A = A 1 A 2 L 1 ( L 1 ) 0 m × n A 3 + L 2 ( L 1 )
  • Step 4: Constructing an approximating operator. Based on Step 1, Step 2 and Step 3, we construct an approximating operator A ( N ) as
    A ( N ) = A 1 + ω ^ 0 · L 1 ( R ) A 4 Ψ 1 ( N ) L 1 ( R N ) D ^ 0 A 4 D ^ 1 I N + L 2 ( R N )
    where L 1 R , L i R N , i = 1 , 2 and Ψ 1 ( N ) are defined in Section 2.3, ⊗ is the Kronecker product and I N is the n-dimensional identity matrix.
  • Step 5: Approximating the dynamic solution. Once the approximating operator A ( N ) is obtained, by giving the initial conditions P 0 of (6), then the dynamic solution of the sytem can be approximated by calculated E N P ( N ) ( t ) = E N e A ( N ) t P N P 0 , where E N and P N are defined in (19) and (20).
By using the algorithm, we can obtain the approximating solution for the system in bounded interval [ 0 , T ] . Moreover, the convergence Theorem 6 and Corollary 1 also indicate that the approximating solution can approximate to the transient solution with arbitrary precision as N .

4.2. Numerical Example

In this subsection, two numerical examples will be presented. Firstly, the software rejuvenation system in [10] is taken as the first example to show the usage for the algorithm. Then, we take a two-state system as the second example, in which the approximation solutions derived by the algorithm and the finite-differences method are considered.
Example 1.
According to [10], we take an software rejuvenation system as an example and investigate the dynamic solution under the algorithm. The software rejuvenation system can be formulated as the following coupled integro differential equations
d P 0 ( t ) d t = ( α 0 + β ) P 0 ( t ) + 0 T μ ( x ) p 2 ( t , x ) d x d P 1 ( t ) d t = β P 0 ( t ) α 1 P 1 ( t ) p 2 ( t , x ) x + p 2 ( t , x ) t = μ ( x ) p 2 ( t , x )
with the boundary condition
p 2 ( t , 0 ) = α 0 P 0 ( t ) + α 1 P 1 ( t ) ,
where α 0 , α 1 and β are constants. μ ( x ) stands for the time-dependent repair rate with an elapsed repair time x [ 0 , T ] and satisfies condition (3). P 0 ( t ) , P 1 ( t ) and p 2 ( t , x ) denotes the probability of the system in the normal operation state, failure probable state and the failure state with elapsed repair time of x at time t, respectively.
By defining the state space X = P = ( P 0 , P 1 , p 2 ( x ) ) τ R 2 × L 1 [ 0 , T ] :   P X = | P 0 | + | P 1 | + p 2 ( x ) L 1 [ 0 , T ] < , and the system operators as follows,
A = α 0 + β 0 0 0 α 1 0 0 0 d d x μ ( x ) , B P 0 P 1 p 2 ( x ) = 0 T μ ( x ) p 2 ( x ) d x β P 0 0 ,
D ( A + B ) = ( P 0 , P 1 , p 2 ( x ) ) T X | p 2 ( x ) is absolutely continuous function , d p 2 ( x ) d x L 1 [ 0 , T ] , p 2 ( 0 ) = α 0 P 0 + α 1 P 1 ,
The software rejuvenation system (21) and (22) can be described as an abstract Cauchy problem in X:
d P ( t ) d t = ( A + B ) P ( t ) P ( 0 ) = P 0
It can be verified that system (23) is well-posed and has a time-dependent solution. Then, we construct an approximation scheme for the abstract version (23). First, we decomposite the system operator A + B and define the boundary operator A 4 for system abstraction. We have
A + B = A 1 A 2 L 1 ( L 1 ) 0 m × n A 3 + L 2 ( L 1 ) , A 4 = α 0 0 0 α 1 ,
where
A 1 = α 0 + β 0 β α 1 , L 1 ( L 1 ) = μ ( x ) 0 , L 2 ( L 1 ) = μ ( x ) ,
and A 2 and A 3 denotes the integral operator and differential operator, respectively. Then, the approximating equation can be constructed as
d P 0 ( N ) ( t ) d t = ( α 0 + β ) P 0 ( N ) ( t ) + μ ( 0 ) · ( α 0 P 0 ( N ) ( t ) + α 1 P 1 ( N ) ( t ) ) · ω ^ 0 + i = 1 N μ ( x i ) p 2 ( N ) ( t , x ^ i ) · ω ^ i d P 1 ( N ) ( t ) d t = β P 0 ( N ) ( t ) α 1 P 1 ( N ) ( t ) d p 2 ( N ) ( t , x ^ 1 ) d t = d ^ 10 ( 1 ) · ( α 0 P 0 ( N ) ( t ) + α 1 P 1 ( N ) ( t ) ) i = 1 N d ^ 1 i ( 1 ) p 2 ( N ) ( t , x ^ i ) μ ( x ^ 1 ) p 2 ( N ) ( t , x ^ 1 ) : d p 2 ( N ) ( t , x ^ N ) d t = d ^ N 0 ( 1 ) · ( α 0 P 0 ( N ) ( t ) + α 1 P 1 ( N ) ( t ) ) i = 1 N d ^ N i ( 1 ) p 2 ( N ) ( t , x ^ N ) μ ( x 1 ) p 2 ( N ) ( t , x N )
where { x ^ j , ω ^ j } j = 0 N be the set of shifted Legendre-Gauss Radau quadrature nodes and weights, ( d k j ( 1 ) ) k , j = 1 , 2 , , N are defined in Equation (16). Moreover, we define the approximating space X ¯ ( N ) = R N + 2 with the norm, for any P ( N ) = P 0 ( N ) , P 1 ( N ) , p 21 ( N ) , , p 2 N ( N ) R N + 2
P ( N ) N = i = 0 1 | P i ( N ) | + ω ^ 0 · ( α 0 P 0 ( N ) + α 1 P 1 ( N ) ) + i = 1 N ω ^ i · p 2 i ( N ) ,
and the approximating operator as
A ( N ) = A 00 μ ( x ^ 0 ) α 1 · ω ^ 0 μ ( x ^ 1 ) · ω ^ 1 μ ( x ^ 2 ) · ω ^ 2 μ ( x ^ N ) · ω ^ N β α 1 0 0 0 d ^ 10 ( 1 ) α 0 d ^ 10 ( 1 ) α 1 ( d ^ 11 ( 1 ) + μ ( x ^ 1 ) ) d ^ 12 ( 1 ) d ^ 1 N ( 1 ) d ^ 20 ( 1 ) α 0 d ^ 20 ( 1 ) α 1 d ^ 21 ( 1 ) ( d ^ 22 ( 1 ) + μ ( x ^ 2 ) ) d ^ 2 N ( 1 ) : : : : : d ^ N 0 ( 1 ) α 0 d ^ N 0 ( 1 ) α 1 d ^ N 1 ( 1 ) d ^ N 2 ( 1 ) ( d ^ N N ( 1 ) + μ ( x ^ N ) ) .
where A 00 = 1 μ ( x ^ 0 ) ω ^ 0 · α 0 β . Then, by denoting P ( N ) ( t ) = ( P 0 ( N ) ( t ) , P 1 ( N ) ( t ) , p 2 ( N ) ( t , x ^ 1 ) , , p 2 ( N ) ( t , x ^ N ) ) , the approximating system for (23) can be described as the following Cauchy problem
d d t P ( N ) ( t ) = A ( N ) P ( N ) ( t ) P ( N ) ( 0 ) = P 0 ( N ) .
Therefore, by denoting the approximation solution on bounded t-intervals as P ^ ( N ) ( t ) , we can obtain that
P ^ ( N ) ( t ) = E N P ( N ) ( t ) = E N · e A ( N ) t · P 0 ( N ) .
Based on above analysis, we propose some numerical result as follows. According to [10], let α 0 = 1 120 , α 1 = 1 3 , β = 1 7 . Let the maximum repair time be T = 20 , then the hazard rate of the repair time μ ( x ) are assumed to satisfies
0 x μ ( s ) d s = 0 x 2 s d s , x < T , x = T .
By solving the Equation (26), we can derive the approximating solution of system (25) directly. Figure 1 shows dynamical behaviour of the approximation solution P 0 ( N ) ( t ) , P 1 ( N ) ( t ) and P 2 ( N ) ( t ) , which corresponding to the probability distributions of the software system in [ 0 , 20 ) with N = 15 , where
P 2 ( N ) ( t ) = ω ^ 0 · ( α 0 P 0 ( N ) ( t ) + α 1 P 1 ( N ) ( t ) ) + i = 1 N ω ^ i · p 2 ( N ) ( t , x ^ i ) .
Figure 1 indicates that the probability of the system in each state will gradually tend to stable as time increases, which corresponds to the fact that the solution of system (21) and (22) gradually tends to a steady-state solution as t tends to infinity. On the other hand, Figure 2 plots the dynamical behaviour of p 2 ( N ) ( t , x ) with N = 15 . Figure 2 shows that for any fixed t, p 2 ( N ) ( t , x ) 0 as x grows. This situation agrees with p 2 ( N ) ( t , x ) L 1 [ 0 , T ] for any given time t and also conforms to reality. As a result, Figure 1 and Figure 2 show that, by using a few nodes, the approximation solution for the software rejuvenation system (21) and (22) on bounded t-intervals can be derived.
Example 2.
In this example, we consider the approximation scheme proposed in this paper as well as the finite-differences method in numerical solution computation to the following first-order linear IDEs with two system state.
d P 0 ( t ) d t = λ P 0 ( t ) + 0 h ( x ) · p 1 ( t , x ) d x p 1 ( t , x ) x + p 1 ( t , x ) t = h ( x ) · p 1 ( t , x )
with the boundary condition
p 1 ( t , 0 ) = λ P 0 ( t ) .
System (27) and (28) can be regarded as a classical two-state repairable system, see [6]. The failure time of a component follows an exponential distribution with rate λ and the repair time follows a general distribution with repair rate h ( x ) , where x stands for the elapsed repair time. Also, System (27) and (28) can be regarded as an single server queuing system with general service rate h ( x ) , where x stands for the elapsed service time, see [5].
By assuming that the repair time (or service time) follows an Erlang distribution with h ( x ) = λ 2 x 1 + λ x , the analytical solution for system (27) and (28) can be calculated as
P 0 ( t ) = 1 3 + 2 3 e 3 2 λ t cos 3 2 λ t .
Moreover, we assume that the maximum repair time (or service time) is T. Then, by performing the approximation algorithm according to Section 4.1, we can derive the corresponding approximation operator A ¯ ( N ) as
A ¯ ( N ) = λ λ 2 x ^ 1 1 + λ x ^ 1 · ω ^ 1 λ 2 x ^ 2 1 + λ x ^ 2 · ω ^ 2 λ 2 x ^ N 1 + λ x ^ N · ω ^ N d ^ 10 ( 1 ) λ ( d ^ 11 ( 1 ) + λ 2 x ^ 1 1 + λ x ^ 1 ) d ^ 12 ( 1 ) d ^ 1 N ( 1 ) d ^ 20 ( 1 ) λ d ^ 21 ( 1 ) ( d ^ 22 ( 1 ) + λ 2 x ^ 2 1 + λ x ^ 2 ) d ^ 2 N ( 1 ) : : : : d ^ N 0 ( 1 ) λ d ^ N 1 ( 1 ) d ^ N 2 ( 1 ) ( d ^ N N ( 1 ) + λ 2 x ^ N 1 + λ x ^ N ) .
Therefore, by denoting the initial condition as P 0 , the approximation solution P ^ ( N ) ( t ) = ( P 0 ( N ) , P 1 ( N ) ) τ for the liner IDEs system (27) and (28) on bounded t-intervals can be obtain as
P ^ ( N ) ( t ) = E N · e A ¯ ( N ) t · P N P 0 .
Let λ = 0.8 , T = 30 . Figure 3 plots the behavior of solution P 0 ( N ) ( t ) with different N in interval [ 0 , 10 ] . Table 1 shows the corresponding error analysis of the approximation solution derived by (30) in different order. Figure 3 and Table 1 show that P 0 ( N ) ( t ) will gradually tends to the analytical solution P 0 ( t ) as N increases, in particular, the curves of P 0 ( 20 ) ( t ) are almost coincident with P 0 ( t ) , which agrees with the conclusion of Corollary 1.
On the other hand, we consider the approximation solution obtained by the finite-differences method. Let λ = 0.8 , T = 30 and the difference notes are denoted as N. And let P ¯ 0 ( N ) ( t ) stands for the approximation solution obtained by the finite-differences method with N difference notes.
Figure 4 plots the behavior of solution P ¯ 0 ( N ) ( t ) which derived by the finite-difference methods with different N in interval [ 0 , 10 ] . Table 2 shows the corresponding error analysis of the approximation solution. Figure 4 and Table 2 shows that the approximation solution P ¯ 0 ( N ) ( t ) will also gradually tends to the analytical solution P 0 ( t ) as N increases. When N = 1000 , the curves of P ¯ 0 ( 20 ) ( t ) are almost coincident with P 0 ( t ) .
Although the finite-differences method can also provide a approximation solution for the system, we can find that the approximation error of P 0 ( 20 ) ( t ) is obviously smaller that P ¯ 0 ( 1000 ) ( t ) , which shows that the convergence rate of the finite-difference method is obviously lower than the convergence rate of the proposed method by comparing with Table 1 and Table 2. However, by comparing with Figure 3 and Figure 4, the approximation scheme proposed in this paper not only exhibits a faster convergence than the finite-differences method, but also requires fewer nodes.

5. Conclusions

In this paper, a polynomials approximation scheme based on the shifted Legendre spectral collocation method was proposed for the numerical computation of first-order linear Integro-Differential Equations (IDEs) in L 1 space. We transformed the first-order linear IDEs into an abstract Cauchy problem to establish the basic framework of theoretical analysis. By decomposing the system operator, we rewrote the Cauchy problem into a more general form for approximation. Then, we constructed an approximation scheme based on the shifted Legendre spectral collocation method and demonstrated the convergence of the approximation scheme in the sense of L 1 -norm by employing Trotter-Kato theorem. The convergence theorem derived in this paper pointed out that for any bounded t-intervals, the system solution can be approximated with arbitrary precision. Moreover, we summarized the approximation scheme into an algorithm for application and presented two numerical examples.
The computational scheme proposed in this paper inherits the advantages of the spectral collocation method with computational efficiency and easily executive. As a result, it provides us with a useful tool to calculate the numerical solution for the first-order linear IDEs. In addition, it also can be applied in computation of uncertainty problems, such as the fuzzy integro-differential equations shown in [39,40].
In this paper, we proved the convergence of the proposed scheme in the system space of R + × L 1 [ 0 , T ] . However, if the elapsed time x belongs to [ 0 , ) , then, the system space will become R + × L 1 [ 0 , ) and the convergence of the proposed scheme will not be guaranteed any more. Therefore, it is a significant job to design an innovative approximation scheme in L 1 [ 0 , ) , which will become a new research direction in the future.

Author Contributions

Investigation, methodology, writing—original draft, Z.C. and H.X.; writing—review and editing, Z.C. and H.H. All authors have read and agreed to the published version of the manuscript.

Funding

The project is funded by the NSAF (Grant No. U1430125).

Data Availability Statement

Not applicable.

Acknowledgments

The authors are thankful to the referees for their constructive suggestions which helped to improve the presentation of this paper.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Yang, X.J. General Fractional Derivatives: Theory, Methods and Applications; Chapman and Hall/CRC: Boca Raton, FL, USA, 2019; pp. 52–85. ISBN 978-04-2928-408-3. [Google Scholar]
  2. Shiri, B. A note on using the Differential Transformation Method for the Integro-Differential Equations. Appl. Math. Comput. 2013, 219, 7306–7309. [Google Scholar] [CrossRef]
  3. Hariharan, G.; Kannan, K. Review of wavelet methods for the solution of reaction-diffusion problems in science and engineering. Appl. Math. Model. 2014, 38, 799–813. [Google Scholar] [CrossRef]
  4. Iannelli, M.; Milner, F. The Basic Approach to Age-Structured Population Dynamics, 1st ed.; Springer: Dordrecht, The Netherlands, 2017; pp. 10–250. ISBN 978-94-0241-146-1. [Google Scholar]
  5. Shortle, J.F.; Thompson, J.M.; Gross, D.; Harris, C.M. Fundamentals Of Queueing Theory, 3rd ed.; John Wiley and Sons: New York, NY, USA, 2018; pp. 25–75. ISBN 978-71-1556-998-1. [Google Scholar]
  6. Ascher, H.; Feingold, H. Repairable Systems Reliability: Modeling, Inference, Misconceptions, and Their Causes; Wiley Online Library: Hoboken, NJ, USA, 1985; pp. 119–223. ISBN 0-8247-7276-8. [Google Scholar]
  7. Sharma, G.; Rai, R.N. Reliability modeling and analysis of environmental control and life support systems of space stations: A literature survey. Acta Astronaut. 2018, 155, 238–246. [Google Scholar] [CrossRef]
  8. Cox, D. The analysis of non-Markovian stochastic processes by the inclusion of supplementary variables. Math. Proc. Camb. 1955, 51, 433–441. [Google Scholar] [CrossRef]
  9. Li, Y.; Meng, X.Y. Reliability analysis of a warm standby repairable system with priority in use. Appl. Math. Model. 2011, 35, 4295–4303. [Google Scholar] [CrossRef]
  10. Huo, H.; Xu, H.; Chen, Z. Modelling and dynamic behaviour analysis of the software rejuvenation system with periodic impulse. Math. Comp. Model. Dyn. 2021, 27, 522–542. [Google Scholar] [CrossRef]
  11. Zheng, F.; Xu, S.S.; Li, X. Numerical solution of the steady-state probability and reliability of a repairable system with three unites-ScienceDirect. Appl. Math. Comput. 2015, 263, 251–267. [Google Scholar] [CrossRef]
  12. Nazarov, A.; Melikov, A.; Pavlova, E. Analyzing an M/M/N Queueing System with Feedback by the Method of Asymptotic Analysis. Cybern. Syst. Anal. 2021, 57, 57–65. [Google Scholar] [CrossRef]
  13. Yuan, L.; Guan, J. An optimal repair–replacement policy for a cold standby system with use priority. Appl. Math. Model. 2013, 35, 1222–1230. [Google Scholar] [CrossRef]
  14. Zong, S.; Chai, G.; Zhe, G. Optimal replacement policy for a deteriorating system with increasing repair times. Appl. Math. Model. 2013, 37, 9768–9775. [Google Scholar] [CrossRef]
  15. Gupur, G. Functional Analysis Methods for Reliability Models; Springer: Basel, Switzerland, 2011; pp. 59–133. ISBN 978-3-0348-0101-0. [Google Scholar]
  16. Pazy, A. Semigroups of Linear Operator and Applications to Partial Differential Equations; Springer: New York, NY, USA, 1983; pp. 8–22. ISBN 978-01-4612-5563-8. [Google Scholar]
  17. Xu, H.; Yu, J.; Zhu, G. Asymptotic property of a reparable multi-state device. Quart. Appl. Math. 2005, 63, 779–789. [Google Scholar] [CrossRef] [Green Version]
  18. Hu, W. Differentiability and compactness of the C0-semigroup generated by the reparable system with finite repair time. J. Math. Anal. Appl. 2016, 433, 1614–1625. [Google Scholar] [CrossRef]
  19. Wang, W.L.; Xu, G.Q. Stability analysis of a complex standby system with constant waiting and different repairman criteria incorporating environmental failure. Appl. Math. Model. 2009, 33, 724–743. [Google Scholar] [CrossRef]
  20. Rhodes, C.A.; House, T. The rate of convergence to early asymptotic behaviour in age-structured epidemic models. Theor. Popul. Biol. 2013, 85, 58–62. [Google Scholar] [CrossRef] [Green Version]
  21. Fu, Z.; Zhu, G.; Chao, G. Well-posedness and stability of the repairable system with N failure modes and one standby unit. J. Math. Anal. Appl. 2011, 375, 174–184. [Google Scholar] [CrossRef] [Green Version]
  22. Funaro, D. Polynomial Approximation of Differential Equations; Springer: New York, NY, USA, 1992; ISBN 3-540-55230-8. [Google Scholar]
  23. Mastroianni, G.; Milovanovic, G. Interpolation Processes; Springer: New York, NY, USA, 2008; pp. 48–68. ISBN 978-35-4068-346-9. [Google Scholar]
  24. Jie, S.; Tao, T.; Wang, L.L. Spectral Methods: Algorithms, Analysis and Applications; Springer: New York, NY, USA, 2011; pp. 47–105. ISBN 978-3-540-71041-7. [Google Scholar]
  25. Dzhumabaev, D.S. On one approach to solve the linear boundary value problems for Fredholm integro-differential equations. J. Comput. Appl. Math. 2016, 294, 342–357. [Google Scholar] [CrossRef]
  26. Dzhumabaev, D.S. New general solutions to linear Fredholm integro-differential equations and their applications on solving the boundary value problems. J. Comput. Appl. Math. 2018, 327, 79–108. [Google Scholar] [CrossRef]
  27. Dzhumabaev, D.S. Computational methods of solving the boundary value problems for the loaded differential and Fredholm integro-differential equations. Math. Methods Appl. 2018, 41, 1439–1462. [Google Scholar] [CrossRef]
  28. Rahmoune, A. Spectral collocation method for solving Fredholm integral equations on the half-line. Appl. Math. Comput. 2013, 219, 9254–9260. [Google Scholar] [CrossRef]
  29. Al-Ahmad, S.; Sulaiman, I.B.; Nawi, M.A.A.; Mamat, M.; Ahmad, M.Z. Analytical solution of systems of Volterra integro-differential equations using modified differential transform method. J. Math. Comput. Sci. 2022, 26, 1–9. [Google Scholar] [CrossRef]
  30. Xu, M.M.; Sulaiman, J.; Ali, L.H. Half-sweep SOR iterative method using linear rational finite difference approximation for first-order Fredholm integro-differential equations. Int. J. Math. Comput. Sci. 2021, 16, 1555–1570. [Google Scholar]
  31. Dawood, L.A.; Hamoud, A.A.; Mohammed, N.M. Laplace discrete decomposition method for solving nonlinear Volterra-Fredholm integro-differential equations. J. Math. Comput. Sci. 2020, 2, 158–163. [Google Scholar] [CrossRef]
  32. Iskandarov, S. Estimate and asymptotic smallness of solutions of a weakly nonlinear implicit Volterra integro-differential equation of the first order on the semiaxis. Lobachevskii J. Math. 2021, 42, 3645–3651. [Google Scholar] [CrossRef]
  33. Xu, H.; Hu, W. Analysis and approximation of a reliable model. Appl. Math. Model. 2013, 37, 3777–3788. [Google Scholar] [CrossRef]
  34. Xu, H.; Hu, W. Modelling and analysis of repairable systems with preventive maintenance. Appl. Math. Comput. 2013, 224, 46–53. [Google Scholar] [CrossRef]
  35. Boardman, N.; Hu, W.; Mishra, R. Optimal Maintenance Design for a Simple Reparable System. In Proceedings of the 2019 IEEE 58th Conference on Decision and Control (CDC), Nice, France, 11–13 December 2019; pp. 3098–3103. [Google Scholar] [CrossRef]
  36. Guo, B.Y.; Wang, Z.Q. Legendre-Gauss collocation methods for ordinary differential equations. Adv. Comput. Math. 2009, 30, 249–280. [Google Scholar] [CrossRef]
  37. Lax, P.D.; Richtmyer, R.D. Survey of the stability of linear finite differential equations. Commun. Pure Appl. Math. 1956, 9, 267–293. [Google Scholar] [CrossRef]
  38. Ito, K.; Kappel, K. The Trotter-Kato theorem and approximation of PDEs. Math. Comput. 1998, 67, 21–44. [Google Scholar] [CrossRef] [Green Version]
  39. Issa, M.B.; Hamoud, A.; Ghadle, K. Numerical solutions of fuzzy integro-differential equations of the second kind. J. Math. Comput. Sci. 2021, 23, 67–74. [Google Scholar] [CrossRef]
  40. Ghanbari, M. A new computational method for solving the first order linear fuzzy Fredholm integro-differential equations. J. Interpolat. Approx. Sci. Comput. 2013, 13, 89–101. [Google Scholar]
Figure 1. Probability distributions of the software system.
Figure 1. Probability distributions of the software system.
Mathematics 10 04117 g001
Figure 2. Behavior of solution p 2 ( 15 ) ( t , x ) .
Figure 2. Behavior of solution p 2 ( 15 ) ( t , x ) .
Mathematics 10 04117 g002
Figure 3. Approximation solutions P 0 ( N ) ( t ) derived by (30) and the analytical solution P 0 ( t ) .
Figure 3. Approximation solutions P 0 ( N ) ( t ) derived by (30) and the analytical solution P 0 ( t ) .
Mathematics 10 04117 g003
Figure 4. Approximation solution P ¯ 0 ( N ) ( t ) obtained by the finite-differences method and the analytical solution P 0 ( t ) .
Figure 4. Approximation solution P ¯ 0 ( N ) ( t ) obtained by the finite-differences method and the analytical solution P 0 ( t ) .
Mathematics 10 04117 g004
Table 1. The error analysis of solutions derived by the proposed scheme in different order.
Table 1. The error analysis of solutions derived by the proposed scheme in different order.
Analytical SolutionApproximation Solutions P 0 ( N ) ( t ) Derived by (30)
P 0 ( t ) P 0 ( 12 ) ( t ) P 0 ( 16 ) ( t ) P 0 ( 20 ) ( t )
ValueValueErrorValueErrorValueError
t = 1 0.487835 0.489046 0.002476 0.486001 0.003772 0.486002 0.003772
t = 2 0.344467 0.330438 0.042456 0.341238 0.009463 0.347338 0.008265
t = 3 0.324477 0.309189 0.049446 0.329176 0.014272 0.322785 0.005242
t = 4 0.328218 0.320813 0.023081 0.329160 0.002860 0.327253 0.002948
t = 5 0.331766 0.327079 0.014328 0.328372 0.010334 0.333693 0.005775
t = 6 0.333070 0.325695 0.022643 0.331923 0.003456 0.332196 0.002631
t = 7 0.333353 0.322681 0.033072 0.334367 0.003031 0.333290 0.000190
t = 8 0.333366 0.321002 0.038516 0.333456 0.000269 0.333566 0.000599
t = 9 0.333346 0.320441 0.040272 0.332444 0.002713 0.333125 0.000664
t = 10 0.333336 0.320042 0.041537 0.332729 0.001825 0.333548 0.000634
Table 2. The error analysis of solutions derived by the finite-differences method in different order.
Table 2. The error analysis of solutions derived by the finite-differences method in different order.
Analytical SolutionApproximation Solutions P ¯ 0 ( N ) ( t ) by Finite-Differences Method
P 0 ( t ) P ¯ 0 ( 100 ) ( t ) P ¯ 0 ( 500 ) ( t ) P ¯ 0 ( 1000 ) ( t )
ValueValueErrorValueErrorValueError
t = 1 0.487835 0.506053 0.036001 0.491842 0.008147 0.489861 0.004136
t = 2 0.344467 0.366708 0.060649 0.349087 0.013235 0.346785 0.006685
t = 3 0.324477 0.342710 0.053201 0.328144 0.011174 0.326310 0.005616
t = 4 0.328218 0.343983 0.045829 0.331388 0.009566 0.329803 0.004806
t = 5 0.331766 0.346979 0.043844 0.334856 0.009229 0.333313 0.004642
t = 6 0.333070 0.348381 0.043946 0.336197 0.009299 0.334637 0.004680
t = 7 0.333353 0.348799 0.044283 0.336510 0.009380 0.334935 0.004721
t = 8 0.333366 0.348873 0.044449 0.336534 0.009413 0.334953 0.004738
t = 9 0.333346 0.348869 0.044494 0.336516 0.009419 0.334934 0.004741
t = 10 0.333336 0.348860 0.044497 0.336506 0.009418 0.334924 0.004740
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Chen, Z.; Xu, H.; Huo, H. Computational Scheme for the First-Order Linear Integro-Differential Equations Based on the Shifted Legendre Spectral Collocation Method. Mathematics 2022, 10, 4117. https://doi.org/10.3390/math10214117

AMA Style

Chen Z, Xu H, Huo H. Computational Scheme for the First-Order Linear Integro-Differential Equations Based on the Shifted Legendre Spectral Collocation Method. Mathematics. 2022; 10(21):4117. https://doi.org/10.3390/math10214117

Chicago/Turabian Style

Chen, Zhuoqian, Houbao Xu, and Huixia Huo. 2022. "Computational Scheme for the First-Order Linear Integro-Differential Equations Based on the Shifted Legendre Spectral Collocation Method" Mathematics 10, no. 21: 4117. https://doi.org/10.3390/math10214117

APA Style

Chen, Z., Xu, H., & Huo, H. (2022). Computational Scheme for the First-Order Linear Integro-Differential Equations Based on the Shifted Legendre Spectral Collocation Method. Mathematics, 10(21), 4117. https://doi.org/10.3390/math10214117

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop