- Email: [email protected]
Two-dimensional integral–algebraic systems: Analysis and computational methods
Journal of Computational and Applied Mathematics 236 (2011) 132–140
Contents lists available at ScienceDirect
Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam
Two-dimensional integral–algebraic systems: Analysis and computational methods M.V. Bulatov a , P.M. Lima b,∗ a
664033, Irkutsk, Lermontov St., 134, Institute of System Dynamics and Control Theory SB RAS, Russia
b
Centro de Matemática e Aplicações, Instituto Superior Técnico, Universidade Técnica de Lisboa, Av. Rovisco Pais, 1049-001 Lisboa, Portugal
article
abstract
info
Article history: Received 7 February 2011 Received in revised form 10 May 2011 MSC: 65R20 45F15
In this article we formulate sufficient conditions for the existence and uniqueness of solution to systems of two-dimensional Volterra integral equations, in which the coefficient of the main term is a singular matrix. A numerical method is introduced which can be applied to approximate the solution when the given conditions are satisfied. The convergence of this method is proved and illustrated by numerical examples. © 2011 Elsevier B.V. All rights reserved.
Keywords: Integral–algebraic systems Two-dimensional Volterra integral equations Implicit Euler method Right rectangles rule
1. Introduction In this work we are concerned with two-dimensional linear systems of Volterra integral equations (VIE), for which the determinant of the matrix in the main term is zero on the whole considered domain. Such problems may arise from the analysis of some classes of differential–algebraic systems of partial differential equations [1] and in the modeling of certain heat conduction processes. Sets of coupled two-dimensional VIE of first and second kind, including algebraic relationships, can be considered as particular cases of such systems. As an important particular case of the integral–algebraic equations (IAE) we mention the Volterra integral equations of the first kind:
∫ t∫ 0
x
K (t , x, τ , s)u(τ , s)dsdτ = f (t , x),
(1)
0
in the rectangular domain
Ω = {0 ≤ s ≤ x ≤ a, 0 ≤ τ ≤ t ≤ b}, where K (t , x, τ , s) is a given matrix, K (t , x, t , x) ̸= 0,
(det K (t , x, t , x) ̸= 0) ∀(t , x) ∈ Ω .
Such equations and systems have been object of research since a long time (see, for example [2,3] and the references therein).
∗
Corresponding author. Tel.: +351 218417071; fax: +351 218417048. E-mail addresses: [email protected] (M.V. Bulatov), [email protected] (P.M. Lima).
0377-0427/$ – see front matter © 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.cam.2011.06.001
M.V. Bulatov, P.M. Lima / Journal of Computational and Applied Mathematics 236 (2011) 132–140
133
On the other hand, some systems of coupled partial differential equations and algebraic relations can also be rewritten in the form of IAE. As an example, consider the system utx (t , x) + a(t , x)u(t , x) + b(t , x)v(t , x) = f (t , x), c (t , x)ut (t , x) + d(t , x)v(t , x) = g (t , x), where a(.), b(.), c (.), d(.) are matrices of a certain dimension, u(.), v(.) are the unknown functions, f (.), g (.) are given functions of the corresponding dimension. The qualitative behavior of such systems has a number of specific characteristics, which must be taken into account in the construction of efficient numerical methods. In particular, such systems may have more than one solution, as well as they may have no solution at all. Here we call a solution a sufficiently smooth vector function (with respect to all its arguments), which once replaced in the system of equations transforms it into an equality. Even when the considered problem has a unique solution some numerical methods may lead to unstable algorithms. As far as we know, the subject of two-dimensional integro-algebraic systems of equations is new in the literature and there are no available results that we could use as a starting point. Therefore, the main purpose of this work is to make a first step towards the study of this class of equations. In Section 2, after formulating the problem and introducing some definitions and preliminary results, we prove the main theorem on existence and uniqueness of solution. In Section 3 we describe a computational method for the numerical solution of the considered class of problems and prove its convergence. Numerical results illustrating the performance of this algorithm are presented in Section 4, as well as the main conclusions of this work. 2. Existence and uniqueness of solution Consider the system A(t , x)u(t , x) +
∫ t∫ 0
x
K (t , x, τ , s)u(τ , s)dsdτ = f (t , x)
(2)
0
in the rectangular domain
Ω = {0 ≤ s ≤ x ≤ a, 0 ≤ τ ≤ t ≤ b}, where A(t , x) and K (t , x, τ , s) are (n × n) matrices, f (t , x) is a given vector function and u(t , x) is the unknown. In this paper we assume that system (1) satisfies det A(t , x) = 0,
∀(t , x) ∈ Ω .
(3)
We shall denominate such systems as two-dimensional integral–algebraic equations (TIAE). Remark that in the case n = 1 an equation of the form (1) satisfying condition (2) is a first kind Volterra integral equation. An equation of the form A(t )u(t ) +
t
∫
K (t , τ )u(τ )dτ = f (t ),
0 ≤ τ ≤ t ≤ 1,
0
satisfying the condition det A(t ) ≡ 0 is usually called an integral–algebraic equation. General information about this type of equations can be found in [4–8]. Let us describe some of the main properties of the considered class of problems.
• A TIAE of the form (1) may have more than one solution, or no solution at all. • We say that the point (t0 , x0 ) is nonsingular if for each n-dimensional vector u0 there is a unique vector function u(t , x), such that u is a solution of the given equation and u(t0 , x0 ) = u0 ; otherwise, (t0 , x0 ) is said to be singular. Under this definition, the points where the rank of matrix A(t , x) changes may be either singular or nonsingular. We will now give some examples that illustrate these properties. Example 1. We start by considering the following system:
ϕ(t , x) 0
∫ t∫ x v(t , x) −1 + w(t , x) 0 0 0
0 0
v(τ , s) 0 dsdτ = . w(τ , s) f (t , x)
0 1
In order to assure that this system is solvable we require that f (0, x) = f (t , 0) = 0. In this case we have w(t , x) = f ′′ tx . If we define ϕ(t , x) = tx2 , we see that all the functions of the form v(t , x) = α(β t + γ x), where α, β and γ are arbitrary real numbers, are solutions of the considered equation and satisfy v(0, 0) = 0, which means that the origin is a singular point. Example 2. Consider now the system
0 0
∫ t∫ x ϕ(t , x) v(t , x) −1 + 0 w(t , x) 0 0 0
v(τ , s) 0 dsdτ = , w(τ , s) f (t , x)
0 1
134
M.V. Bulatov, P.M. Lima / Journal of Computational and Applied Mathematics 236 (2011) 132–140
where f (0, x) = f (t , 0) = 0, which is very close to the first example; this system has as its unique solution w(t , x) =
f ′′ tx , v(t , x) = ∂∂t ∂ x [ϕ(t , x)f ′′ tx ]. In this case, if ϕ(t , x) satisfies the appropriate smoothness conditions the points where this function vanishes are nonsingular. 2
Example 3. Consider the homogeneous system
1 0
∫ t∫ x 0 v(t , x) + 1 w(t , x) 0 0
0 0
1
v(τ , s) 0 dsdτ = . w(τ , s) 0
d(t − τ )(x − s)
Differentiating both sides of this system with respect to t and x, we obtain
v(t , x) + d
x
∫ t∫ 0
w(τ , s)dsdτ = 0, 0
which in the case d = 1 coincides with the first equation of the considered system. Hence in the case d = 1 this system has a family of solutions of the form v(t , x) = ψ(t , x), w(t , x) = −ψ ′′ tx (t , x), where ψ(t , x) is any function with continuous second order mixed derivative. Some other properties of the systems of this class are: (a) they may not have any solutions (b) small perturbations (in the sense of the maximum norm) in the right-hand side may lead to completely different solutions (see [8]). Example 4. We will now analyze the system
1 0
∫ t∫ x v(t , x) 1 + wδ (t , x) 0 0 0
0 0
v(τ , s) 0 dsdτ = . w(τ , s) δ(t , x)
0 1
If δ(0, x) ̸= 0 or δ(t , 0) ̸= 0, the considered problem has no continuous solution. If δ(t , x) ≡ 0, the unique solution is the trivial one. But if δ(t , x) = ε sin t +ε x we conclude that wδ (t , x) = − 1ε sin t +ε x . We will now present some auxiliary results that are necessary for the theorem we want to prove. Theorem 1 ([9]). The system of second kind integral equations u(t , x) +
t
∫
K1 (t , x, τ , x)u(τ , x)dτ + 0
x
∫
K2 (t , x, t , s)u(t , s)ds + 0
∫ t∫ 0
x
K3 (t , x, τ , s)u(τ , s)dsdτ = f (t , x),
(4)
0
where K1 (·), K2 (·), K3 (·) are (n × n) matrices with continuous elements, f (t , x) is an n-dimensional vector function with continuous elements, has a unique solution. Lemma 1 (See, for Example [2]). Let the sequence ψi+1,j+1 be defined for i = 0, 1, . . . , N − 1 and j = 0, 1, . . . , M − 1, h, q > 0, α, β, γ , l ≥ 0 and
ψi+1,j+1 ≤ l + α h
i −
ψk,j+1 + β q
k=1
j −
ψi+1,p + γ hq
p=1
j i − −
ψk,p .
k=1 p=1
Then the following inequality holds true
ψi+1,j+1 ≤ l exp(α ti + β xj )J0 2I (αβ + γ )ti xj , where ti = ih, xj = jq, h = b/N , q = a/M , J0 (·) is the zero order Bessel function and I stands for the unity of imaginary numbers. Example 5. This example shows that the regularity of the matrix pencil λA(t , x) + K (t , x) does not assure the existence and uniqueness of solution of a IAE, even in the one-dimensional case. Consider the system of one-dimensional equations of the form
1 0
∫ t v(t ) 0 + w(t ) 1 0
t 0
r
2t − τ
v(τ ) 0 dτ = , w(τ ) 0
t ∈ [0, 1],
where r is a scalar. It is easy to verify that for r ̸= 0 the matrix pencil λA(t ) + K (t , t ) is regular, though for r = 2 the given t equation has an infinite set of solutions of the form v(t ) = −t w(t ) − 2 0 w(τ )dτ , where w(t ) is an arbitrary function. Hence, even in the one-dimensional case the regularity (or singularity) of a given matrix pencil λA(t ) + K (t , t ) does not give any information on the existence of a unique solution. Now we will define a condition, independent of the regularity of a matrix pencil, which can be used in the existence and uniqueness theory for TIAE.
M.V. Bulatov, P.M. Lima / Journal of Computational and Applied Mathematics 236 (2011) 132–140
135
Definition 1 (See, for Example, [10]). We say that the matrix pencil λA(t , x) + K (t , x) satisfies the rank-power criterion (or has index one, or has a simple structure) in the domain Ω , if rank A(t , x) = deg det(λA(t , x) + K (t , x)) = k = const, or rank C (t , x) = rank C 2 (t , x) = k = const
∀(t , x) ∈ Ω ,
where C (t , x) = (λA(t , x) + K (t , x)) A(t , x) and deg(·) denotes the degree of a function as a polynomial in λ. If the rank of matrix A(t , x) is constant and equal to k on the domain Ω , it is known [10] that there exists a matrix P (t , x), nonsingular for any (t , x) ∈ Ω , whose elements have the same degree of smoothness as the elements of A(t , x), and which satisfies −1
P (t , x)A(t , x) =
B(t , x) , 0
(5)
where B(t , x) has dimension (k × n), 0 is a null matrix of dimension ((n − k) × n) and rank B(t , x) = k = const∀(t , x) ∈ Ω . Let us denote P (t , x)K (t , x) =
S ( t , x) , Q (t , x)
(6)
where K (t , x) is the same matrix as in Definition 1, P (t , x) is defined by (4), S (t , x) and Q (t , x) are matrices with dimension (k × n) and (n − k × n), respectively. Lemma 2. Let A(t , x) have the form moreover, let K (t , x) =
S (t , x) Q (t , x)
B(t , x) 0
, where B(t , x) is a (k × n) matrix and rank B(t , x) = k = const∀(t , x) ∈ Ω ;
, where S (t , x) and Q (t , x) are (k × n) and (n − k × n) matrices, respectively and the pencil λA(t , x) + K (t , x) satisfies the rank-power criterion. Then B(t , x) det ̸= 0 ∀(t , x) ∈ Ω . Q ( t , x) The proof of this lemma follows from [11]. Let us now formulate the Theorem about existence and uniqueness of solution of a system of the form (1) satisfying condition (2). Theorem 2. Assume that problem (1) satisfies the following conditions: 1. The elements of the matrices A(t , x), K (t , x, τ , s) and f (t , x) are continuous functions, as well as their first order partial derivatives and second order mixed derivatives with respect to t and x. 2. rank A(t , 0) = rank(A(t , 0)|f (t , 0)), rank A(0, x) = rank(A(0, x)|f (0, x)); 3. The matrix pencil λA(t , x) + K (t , x, t , x) satisfies the rank-power criterion in Ω . Then the considered system has a unique solution in Ω . Proof. Substituting into (1) x = 0 or t = 0, we obtain A(t , 0)u(t , 0) = f (t , 0) or A(0, x)u(0, x) = f (0, x), respectively. The solvability of the corresponding systems is guaranteed by condition 2. Multiplying both sides of Eq. (1) by P (t , x) (the matrix considered in (4)) and writing the result in the block form, we obtain
B(t , x) u(t , x) + 0
S (t , x, τ , s) Q (t , x, τ , s)
∫ t∫ x 0
S (t , x, τ , s) u(τ , s)dsdτ = Q (t , x, τ , s)
0
ϕ(t , x) , ψ(t , x)
(7)
= P (t , x)K (t , x, τ , s), (ϕ T (t , x), ψ T (t , x))T = P (t , x)f (t , x). Let us remind that B(·), S (·) and Q (·) are, respectively, (k × n), (k × n) and ((n − k) × n) matrices. Differentiating the second block line of (6) with respect to t and x, where
which is possible according to condition 1, we obtain Q (t , x, t , x)u(t , x) +
x
∫
Q x (t , x, t , s)u(t , s)ds + ′
0
∫ t∫
t
Q ′ t (t , x, τ , x)u(τ , x)dτ
0
x
Q ′′ tx (t , x, τ , s)u(τ , s)dsdτ = ψ ′′ tx (t , x).
+ 0
∫
0
Combining this last equality with the first line of (6) we obtain
B(t , x) u(t , x) + Q (t , x)
∫ t∫ x + 0
0
∫ x
0 ds + Q x (t , x, t , s)u(t , s)
∫ t
0
S (t , x, τ , s) u(τ , s)dsdτ = Q ′′ tx (t , x, τ , s)
0
0
ϕ(t , x) . ′′ ψ tx (t , x)
Q t (t , x, τ , x)u(τ , x) ′
′
ds
136
M.V. Bulatov, P.M. Lima / Journal of Computational and Applied Mathematics 236 (2011) 132–140
From the second condition of the theorem, taking into consideration Lemma 2, we conclude that the matrix nonsingular for (t , x) ∈ Ω . Finally, multiplying both sides of (7) by
B(t , x) Q (t , x)
−1
B(t , x) Q (t , x)
is
, we obtain a system of the form (3); this
system satisfies the conditions of Theorem 1 and hence has a unique solution. This concludes the proof of the theorem.
3. Numerical algorithm We will now describe and analyze a numerical algorithm for the solution of problem (1), assuming that this system satisfies the conditions of Theorem 2. This method is based on the right rectangles cubature formula (implicit Euler Method) and is analogous to the method introduced in [2,3] for the solution of two-dimensional first kind Volterra integral equations (which are a particular case of TIAE, as remarked above). We have decided to start with such an elementary method because very little is known about this class of equations. Moreover, it is very difficult to analyze the stability and convergence of higher order methods when applied to them. Numerical evidence suggests that they are often unstable. Another advantage of the considered method (right rectangles rule) is that it has a regularization effect with respect to rounding-off errors (if the stepsize is correctly chosen). For simplicity, we will illustrate this in the one-dimensional case. Consider for example the problem of computing u from its primitive f : t
∫
u(τ )dτ = f (t ),
t ∈ [0, 1], f (0) = 0, u(t ) = f ′ (t ).
(8)
0
Assume that the right-hand side function f is replaced by g, such that |g (t ) − f (t )| = δ(t ) ≤ δ . Suppose that t
∫
v(τ )dτ = g (t ),
t ∈ [0, 1].
(9)
0
Note that the considered problem is ill-conditioned, since the solution of (9) may be very different from the one of (8), even for small δ (for example, if f (t ) = 0 and g (t ) = δ sin(t /δ 2 )). However, if the perturbed Eq. (9) is solved by the right rectangles formula (analogous to the method, considered in the present article), omitting the algebraic manipulation, we obtain:
|vi − u(ti )| ≤ |(f (ti ) − f (ti−1 ))/h| + 2δ/h ≤ Lh + 2δ/h, where L < ∞. From the last inequality we conclude that the optimal precision (with minimal effect of the input data error) is obtained when h is of the order of δ 1/2 . The fact that the discretization procedure for integral equations of the first kind has a regularization effect, when there is a certain relationship between δ and h, was remarked for the first time in [12]. In the two-dimensional case, for equations that satisfy k(t , x, t , x) ̸= 0, with (t , x) ∈ Ω , a similar result can be proved; in this case, the optimal order for h and q is δ 1/3 [2]. Returning to the construction of a numerical method for TIAE, let us define in the domain Ω the mesh ti = ih, i = 1, 2, . . . , N , h = b/N , xj = jq, j = 1, 2, . . . , M , q = a/M and denote Aij = A(ti , xj ), fij = f (ti , xj ), Kijlm = K (ti , xj , tl , xm ), l ≤ i, m ≤ j; let uij be the approximate value of u(ti , xj ). The numerical method for Eq. (1), based on the right rectangles rule, has the form Aij uij + hq
j i − −
Kijlm ulm = fij .
(10)
l=1 m=1
The main result about the convergence of the described method is given by the following theorem. Theorem 3. Let the following conditions be satisfied: 1. The elements of matrices A(t , x), K (t , x, τ , s) and the vector function f (t , x) are continuous functions, as well as their partial and mixed derivatives, up to the third order; 2. The second and third conditions of Theorem 2 are satisfied. Then the following estimate is true: max ‖uij − u(ti , xj )‖ = O(h + q). i,j
Proof. Replacing in (10) ui,j by the exact value u(ti , xj ), we obtain Aij u(ti , xj ) + hq
j i − − l=1 m=1
Kijlm u(tl , xm ) = fij + Rij ,
(11)
M.V. Bulatov, P.M. Lima / Journal of Computational and Applied Mathematics 236 (2011) 132–140
137
where Rij is the error of the cubature formula of the right rectangles: Rij =
j ∫ i − −
tl
∫
(K (ti , xj , τ , s)u(τ , s) − Kijlm ulm )dsdτ .
(12)
xm−1
tl−1
l=1 m=1
xm
Denoting εij = u(ti , xj ) − uij and subtracting (10) from (11), we conclude that the error εij satisfies the following equation j i − −
Aij εij + hq
Kijlm εlm = Rij .
(13)
l=1 m=1
If we multiply both sides of (13) by Pij = P (ti , xj ) (the same matrix as in (4)), we obtain
Bij 0
εij + hq
j i − − S l=1 m=1
ijlm
Qijlm
εlm = Pij Rij .
(14)
Let us rewrite these equalities in the form of two systems Bij εij + hq
j i − −
Sijlm εlm = ρij ,
(15)
l=1 m=1
hq
j i − −
Qijlm εlm = δij ,
(16)
l=1 m=1
where (ρijT , δijT )T = Pij Rij . Let us introduce the difference operators ∆h and ∆q , which are defined by the equalities
∆h g (ti , xj ) = ∆q g (ti , xj ) =
1 h 1 q
(g (ti , xj ) − g (ti−1 , xj )),
i = 2, 3, . . . , N ,
(g (ti , xj ) − g (ti , xj−1 )),
j = 2, 3, . . . , M ,
for any given vector function g (t , x). These are the simplest difference analogs of the partial derivatives with respect to t and to x, respectively. Applying successively these two operators to both sides of (16), we obtain Qijij εij +
j −1 i−1 − − (Qijlj − Qi−1jlj )εlj + (Qijim − Qij−1im )εim l =1
m=1
i −1 j −1
+
−−
(Qijlm − Qi−1jlm − Qij−1lm + Qi−1j−1lm )εlm =
l=1 m=1
1 h·q
(δij − δij−1 − δi−1j + δi−1j−1 ).
Combining this system with (15) we obtain
Bij + hqSijij Qijij
+
j −1 i −1 − − l=1 m=1
εij +
i −1 − l =1
Qijlj
0 − Qi−1jlj
εlj +
j −1 − m=1
hqSijlm Qijlm − Qi−1jlm − Qij−1lm + Qi−1j−1lm
Qijim
0 − Qij−1im
εlm = 1 hq
εim ρij
, (δij − δij−1 − δi−1j + δi−1j−1 )
(17)
where 0 denotes a (k × n) null matrix. From the third condition of Theorem 2 and Lemma 2 it follows that for sufficiently small h and q the matrix
Bij + hqSijij Qijij
is nonsingular, that is,
−1 Bij + hqSijij ≤ K1 < ∞. Qijij
(18)
Taking into consideration the first condition of the theorem, we conclude that the differences on the left-hand side of (17) satisfy
‖Qijlj − Qi−1jlj ‖ ≤ hK2 ,
K2 < ∞,
(19)
138
M.V. Bulatov, P.M. Lima / Journal of Computational and Applied Mathematics 236 (2011) 132–140
‖Qijim − Qij−1im ‖ ≤ qK3 ,
K3 < ∞,
(20)
‖Qijlm − Qi−1jlm − Qij−1lm − Qi−1j−1lm ‖ ≤ (h + q)K4 ,
K4 < ∞.
(21)
Moreover, from the error estimate for the right rectangles cubature formula it follows that
‖ρij ‖ ≤ (h + q)K5 ,
K5 < ∞.
(22)
According to formula (12) and the first condition of the theorem, after some cumbersome algebraic manipulations (which we omit here) we conclude that 1 hq
‖δij − δij−1 − δi−1j + δi−1j−1 ‖ ≤ (h + q)K6 ,
K6 < ∞.
(23)
Finally, introducing in (17) the estimates (18)–(23), we obtain
‖εij ‖ ≤
K6 K1
(h + q) +
j −1 i−1 i−1 j−1 K2 − K3 − K4 − − ‖εim ‖ + hq ‖εlm ‖. h ‖εlj ‖ + q K1 l = 1 K1 m=1 K1 l=1 m=1
By Lemma 1, from this last inequality we conclude that
‖εij ‖ = O(h + q). This concludes the proof of the theorem.
4. Numerical examples Example 1. Consider the integro-algebraic system A(t , x)u(t , x) +
∫ t∫ 0
x
K (t , x, τ , s)u(τ , s)dsdτ = f (t , x),
(24)
0
where
[
1 A= 0
[ K =
0 1
]
0 , 0
]
1 . 1
The right-hand side function is given by: f 1 ( t , x) = 1 +
x2 t 2 4
,
f2 (t , x) = xt +
x2 t 2 4
.
Let us check that this system satisfies the conditions of Theorem 2. 1. A, K and f are smooth functions of t and x. 2. rank A(t , 0) = rank(A(t , 0)|f (t , 0)) = 1; rank(A(0, x)) = rank(A(0, x)|f (0, x)) = 1. 3. We have
λA + K =
[
λ 1
]
1 1
and therefore det(λA(t , x) + K (t , x, t , x)) = λ − 1. Hence the degree of det(λA + K ) is 1, for t , x > 0, the same as rank A. Since the conditions of Theorem 2 are satisfied, the considered system has a unique continuous solution on any domain
Ω = {(t , x) : 0 ≤ t ≤ a, 0 ≤ x ≤ b}. By direct substitution we can check that this solution is u(t , x) = (v(t , x), w(t , x))T = (1, tx)T . We shall consider in D a uniform mesh, with stepsize h on t and x: ti = ih,
i = 1, . . . , N ,
xj = jh,
j = 1, . . . , M ,
where h = a/N = b/M. With these notations we have applied the numerical method (10). The numerical results are displayed in Table 1. We have considered a = b = 1 and used five different meshes with h = q = 1/4, 1/8, 1/16, 1/32, 1/64. The error norm, given by
‖ϵh,q ‖∞ =
max
i=1,...,N ,j=1,...,M
‖ui,j − u(ti , xj )‖∞ ,
(25)
M.V. Bulatov, P.M. Lima / Journal of Computational and Applied Mathematics 236 (2011) 132–140
139
Table 1 Numerical results for Example 1. h=q
1/4
1/8
1/16
1/32
1/64
ϵhq
0.234375 0.953
0.121094 0.977
0.0615234 0.986
0.0310059 0.994
0.0155064
p
Table 2 Numerical results for Example 2. h=q
1/8
1/16
1/32
1/64
1/128
ϵhq
0.0314114 0.864
0.0172628 0.937
0.00901509 0.970
0.00460285 0.985
0.00232518
p
is given for each different mesh. We also display the estimates of the convergence order: p = log2
‖ϵh,q ‖∞ . ‖ϵh/2,q/2 ‖∞
(26)
The numerical results are in good agreement with the theoretical convergence order, given by Theorem 3. Example 2. Now let A(t , x) =
[
1 ψ(t , x)
K (t , x, τ , s) =
] ψ(t , x) 2 , ψ(t , x)
t +x−τ +s 1
[
(27)
]
0 , 5 exp(t − x + τ + s)
(28)
ψ(t , x) = t + x. In this case the right-hand side function is given by f ( t , x) =
f1 (t , x) 1 + t 2 − x + t 2 x/2 + t 4 x/12 + 3tx2 /2 − t 2 x2 /4 + t 3 x2 /2 − 5tx3 /6 + e−t −x (t + x), = f2 (t , x) tx + 5 exp(t − x)tx + t 3 x/3 − tx2 /2 + (1 + t 2 − x)(t + x) + exp(−t − x)(t + x)2
[
]
[
]
and the exact solution is u(t , x) = (v(t , x), w(t , x))T = (1 + t 2 − x, exp(−t − x))T ,
(29)
as can be verified, by direct substitution. Let us check that this system satisfies the conditions of Theorem 2. 1. A, K and f are smooth functions of t and x. 2. rank A(t , 0) = rank(A(t , 0)|f (t , 0)) = 1; rank(A(0, x)) = rank(A(0, x)|f (0, x)) = 1. 3. We have
] λ(t + x) λ(t + x)2 + 5 exp(2t ) and therefore det(λA(t , x) + K (t , x, t , x)) = λ 5 exp(2t ) + 2x(t + x)2 − (t + x) + 10x exp(2t ). Hence the degree of det(λA(t , x) + K (t , x, t , x)) is 1, for t , x > 0, the same as rank A. λA(t , x) + K (t , x, t , x) =
[
λ + 2x 1 + λ(t + x)
Since the conditions of Theorem 2 are satisfied, we may assure that (29) is the unique solution of the considered system on any domain Ω = {(t , x) : 0 ≤ t ≤ a, 0 ≤ x ≤ b}. The numerical results for this example are displayed in Table 2. As before, we have considered a = b = 1. We have used five different meshes, this time with h = q = 1/8, 1/16, 1/32, 1/64, 1/128. The error norm ϵhq is given by (25) and the empirical estimate of the convergence order p is given by (26). As in the first example, the numerical results confirm the theoretical convergence order. Example 3. Now let A(t , x) =
[
1 ψ(t , x)
K (t , x, τ , s) =
] ψ(t , x) 2 , ψ(t , x)
t +x+τ −s τ + 2s
[
ψ(t , x) = t − x.
(30)
]
0 , 5 exp(t − x − τ + s)
(31)
140
M.V. Bulatov, P.M. Lima / Journal of Computational and Applied Mathematics 236 (2011) 132–140 Table 3 Numerical results for Example 3. h=q
1/8
1/16
1/32
1/64
1/128
ϵhq
0.626502 0.752
0.0371996 0.881
0.020194 0.942
0.0105117 0.971
0.00536162
p
In this case the right-hand side function is given by f (t , x) =
=
f 1 ( t , x) f 2 ( t , x)
[
]
1 + t + e−t −x (t − x) + 3t 2 x/2 + 5t 3 x/6 − x2 + tx2 /2 + t 2 x2 /4 − t 2 x3 /2 − tx4 /12, tx exp(−t − x)(t − x)2 + (t − x)(1 + t − x2 ) + t (3 + 2t ) + 3(2 + t )x − tx2 − 3x3 + 5 exp(−x)x sinh(t ) 6
and the exact solution is u(t , x) = (v(t , x), w(t , x))T = (1 + t − x2 , exp(−t − x))T ,
(32)
as can be verified, by direct substitution. Let us check that this system satisfies the conditions of Theorem 2. 1. A, K and f are smooth functions of t and x. 2. rank A(t , 0) = rank(A(t , 0)|f (t , 0)) = 1; rank(A(0, x)) = rank(A(0, x)|f (0, x)) = 1. 3. We have
λA(t , x) + K (t , x, t , x) =
[
λ + 2t λ(t − x) + t + 2x
λ(t − x) λ(t − x)2 + 5
]
and therefore det(λA(t , x) + K (t , x, t , x)) = λ(5 + 2t 3 + 2x2 + tx(−1 + 2x) − t 2 (1 + 4x)) + 10t. Hence the degree of det(λA(t , x) + K (t , x, t , x)) is 1, for t , x > 0, the same as rank A. Since the conditions of Theorem 2 are satisfied, we may assure that (32) is the unique solution of the considered system on any domain Ω = {(t , x) : 0 ≤ t ≤ a, 0 ≤ x ≤ b}. The numerical results for this example are displayed in Table 3. As in the previous example, we have considered a = b = 1 and used five different meshes, with h = q = 1/8, 1/16, 1/32, 1/64, 1/128. The error norm ϵhq is given again by (25) and the estimated convergence order p is defined by (26). Once again the numerical results are in agreement with Theorem 3. The numerical method (10) was coded as a program in Mathematica. Floating point was used with the default precision (16 digits). All the computations were carried out in a personal computer with a Pentium 2.3 GHz processor. The obtained results confirm that the computational method (10) provides a simple and stable algorithm to approximate the solution of system (1). Its low convergence order is in some sense compensated by the simplicity of computations. We cannot provide here a comparison of methods, since (as far as we know) there were no previous attempts to solve such problems numerically. We believe that certain higher order methods can be applicable. In particular, it seems reasonable to apply a method based on the trapezoidal rule. However, the stability analysis is much more difficult (as remarked above) for higher order methods. This subject deserves further analysis and we leave it as future work. Acknowledgments M. Bulatov acknowledges financial support from RFBR, grant 10-01-00571a; P. Lima acknowledges support from FCT, grant PTDC/MAT/101867/2008. References [1] G.V. Demidenko, S.V. Uspenskii, Partial Differential Equations and Systems not Solvable with Respect to the Highest-Order Derivative, Taylor & Francis, 2005. [2] Ten Men Yan, Approximate solution of linear volterra integral equations of the first kind, Ph.D. Dissertation, Irkutsk, 1985 (in Russian). [3] S. McKee, Tao Tang, T. Diogo, An Euler-type method for two-dimensional Volterra integral equations of the first kind, IMA J. Numer. Anal. 20 (2000) 423–440. [4] H. Brunner, Collocation Methods for Volterra Integral and Related Functional Equations, University Press, Cambridge, 2004. [5] M. Hadizadeh, F. Ghoreishi, S. Pishbin, Jacobi spectral solution for integral algebraic equations of index-2, Appl. Numer. Math. 61 (2011) 131–148. [6] J.-P. Kauthen, The numerical solution of integral-algebraic equations of index-1 by pollinomial spline collocation methods, Math. Comp. 70 (2001) 1503–1514. [7] M.V. Bulatov, V.F. Chistyakov, The properties of differential–algebraic systems and their integral analogs, September, Memorial University of Newfoundland, 1997, Preprint. [8] M.V. Bulatov, Regularization of singular systems of Volterra integral equations, Comput. Math. Math. Phys. 42 (2002) 315–320. [9] V.I. Smirnov, A Course of Higher Mathematics, vol.4, Pergamon Press, Oxford, 1964. [10] V.F. Chistyakov, Differential-Algebraic Operators with Finite-Dimensional Kernel, Nauka, Novosibirsk, 1996 (in Russian). [11] M.V. Bulatov, Transformation of differential–algebraic systems of equations, Comput. Math. Math. Phys. 34 (1994) 301–311. [12] A.S. Apartsyn, A.B. Bakushinsky, Approximate solution of Volterra integral equations of the first kind by the quadratures method, Differential and Integral Equations, V. 1, Irkustk, 1972, pp. 248–258 (in Russian).
Contents lists available at ScienceDirect
Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam
Two-dimensional integral–algebraic systems: Analysis and computational methods M.V. Bulatov a , P.M. Lima b,∗ a
664033, Irkutsk, Lermontov St., 134, Institute of System Dynamics and Control Theory SB RAS, Russia
b
Centro de Matemática e Aplicações, Instituto Superior Técnico, Universidade Técnica de Lisboa, Av. Rovisco Pais, 1049-001 Lisboa, Portugal
article
abstract
info
Article history: Received 7 February 2011 Received in revised form 10 May 2011 MSC: 65R20 45F15
In this article we formulate sufficient conditions for the existence and uniqueness of solution to systems of two-dimensional Volterra integral equations, in which the coefficient of the main term is a singular matrix. A numerical method is introduced which can be applied to approximate the solution when the given conditions are satisfied. The convergence of this method is proved and illustrated by numerical examples. © 2011 Elsevier B.V. All rights reserved.
Keywords: Integral–algebraic systems Two-dimensional Volterra integral equations Implicit Euler method Right rectangles rule
1. Introduction In this work we are concerned with two-dimensional linear systems of Volterra integral equations (VIE), for which the determinant of the matrix in the main term is zero on the whole considered domain. Such problems may arise from the analysis of some classes of differential–algebraic systems of partial differential equations [1] and in the modeling of certain heat conduction processes. Sets of coupled two-dimensional VIE of first and second kind, including algebraic relationships, can be considered as particular cases of such systems. As an important particular case of the integral–algebraic equations (IAE) we mention the Volterra integral equations of the first kind:
∫ t∫ 0
x
K (t , x, τ , s)u(τ , s)dsdτ = f (t , x),
(1)
0
in the rectangular domain
Ω = {0 ≤ s ≤ x ≤ a, 0 ≤ τ ≤ t ≤ b}, where K (t , x, τ , s) is a given matrix, K (t , x, t , x) ̸= 0,
(det K (t , x, t , x) ̸= 0) ∀(t , x) ∈ Ω .
Such equations and systems have been object of research since a long time (see, for example [2,3] and the references therein).
∗
Corresponding author. Tel.: +351 218417071; fax: +351 218417048. E-mail addresses: [email protected] (M.V. Bulatov), [email protected] (P.M. Lima).
0377-0427/$ – see front matter © 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.cam.2011.06.001
M.V. Bulatov, P.M. Lima / Journal of Computational and Applied Mathematics 236 (2011) 132–140
133
On the other hand, some systems of coupled partial differential equations and algebraic relations can also be rewritten in the form of IAE. As an example, consider the system utx (t , x) + a(t , x)u(t , x) + b(t , x)v(t , x) = f (t , x), c (t , x)ut (t , x) + d(t , x)v(t , x) = g (t , x), where a(.), b(.), c (.), d(.) are matrices of a certain dimension, u(.), v(.) are the unknown functions, f (.), g (.) are given functions of the corresponding dimension. The qualitative behavior of such systems has a number of specific characteristics, which must be taken into account in the construction of efficient numerical methods. In particular, such systems may have more than one solution, as well as they may have no solution at all. Here we call a solution a sufficiently smooth vector function (with respect to all its arguments), which once replaced in the system of equations transforms it into an equality. Even when the considered problem has a unique solution some numerical methods may lead to unstable algorithms. As far as we know, the subject of two-dimensional integro-algebraic systems of equations is new in the literature and there are no available results that we could use as a starting point. Therefore, the main purpose of this work is to make a first step towards the study of this class of equations. In Section 2, after formulating the problem and introducing some definitions and preliminary results, we prove the main theorem on existence and uniqueness of solution. In Section 3 we describe a computational method for the numerical solution of the considered class of problems and prove its convergence. Numerical results illustrating the performance of this algorithm are presented in Section 4, as well as the main conclusions of this work. 2. Existence and uniqueness of solution Consider the system A(t , x)u(t , x) +
∫ t∫ 0
x
K (t , x, τ , s)u(τ , s)dsdτ = f (t , x)
(2)
0
in the rectangular domain
Ω = {0 ≤ s ≤ x ≤ a, 0 ≤ τ ≤ t ≤ b}, where A(t , x) and K (t , x, τ , s) are (n × n) matrices, f (t , x) is a given vector function and u(t , x) is the unknown. In this paper we assume that system (1) satisfies det A(t , x) = 0,
∀(t , x) ∈ Ω .
(3)
We shall denominate such systems as two-dimensional integral–algebraic equations (TIAE). Remark that in the case n = 1 an equation of the form (1) satisfying condition (2) is a first kind Volterra integral equation. An equation of the form A(t )u(t ) +
t
∫
K (t , τ )u(τ )dτ = f (t ),
0 ≤ τ ≤ t ≤ 1,
0
satisfying the condition det A(t ) ≡ 0 is usually called an integral–algebraic equation. General information about this type of equations can be found in [4–8]. Let us describe some of the main properties of the considered class of problems.
• A TIAE of the form (1) may have more than one solution, or no solution at all. • We say that the point (t0 , x0 ) is nonsingular if for each n-dimensional vector u0 there is a unique vector function u(t , x), such that u is a solution of the given equation and u(t0 , x0 ) = u0 ; otherwise, (t0 , x0 ) is said to be singular. Under this definition, the points where the rank of matrix A(t , x) changes may be either singular or nonsingular. We will now give some examples that illustrate these properties. Example 1. We start by considering the following system:
ϕ(t , x) 0
∫ t∫ x v(t , x) −1 + w(t , x) 0 0 0
0 0
v(τ , s) 0 dsdτ = . w(τ , s) f (t , x)
0 1
In order to assure that this system is solvable we require that f (0, x) = f (t , 0) = 0. In this case we have w(t , x) = f ′′ tx . If we define ϕ(t , x) = tx2 , we see that all the functions of the form v(t , x) = α(β t + γ x), where α, β and γ are arbitrary real numbers, are solutions of the considered equation and satisfy v(0, 0) = 0, which means that the origin is a singular point. Example 2. Consider now the system
0 0
∫ t∫ x ϕ(t , x) v(t , x) −1 + 0 w(t , x) 0 0 0
v(τ , s) 0 dsdτ = , w(τ , s) f (t , x)
0 1
134
M.V. Bulatov, P.M. Lima / Journal of Computational and Applied Mathematics 236 (2011) 132–140
where f (0, x) = f (t , 0) = 0, which is very close to the first example; this system has as its unique solution w(t , x) =
f ′′ tx , v(t , x) = ∂∂t ∂ x [ϕ(t , x)f ′′ tx ]. In this case, if ϕ(t , x) satisfies the appropriate smoothness conditions the points where this function vanishes are nonsingular. 2
Example 3. Consider the homogeneous system
1 0
∫ t∫ x 0 v(t , x) + 1 w(t , x) 0 0
0 0
1
v(τ , s) 0 dsdτ = . w(τ , s) 0
d(t − τ )(x − s)
Differentiating both sides of this system with respect to t and x, we obtain
v(t , x) + d
x
∫ t∫ 0
w(τ , s)dsdτ = 0, 0
which in the case d = 1 coincides with the first equation of the considered system. Hence in the case d = 1 this system has a family of solutions of the form v(t , x) = ψ(t , x), w(t , x) = −ψ ′′ tx (t , x), where ψ(t , x) is any function with continuous second order mixed derivative. Some other properties of the systems of this class are: (a) they may not have any solutions (b) small perturbations (in the sense of the maximum norm) in the right-hand side may lead to completely different solutions (see [8]). Example 4. We will now analyze the system
1 0
∫ t∫ x v(t , x) 1 + wδ (t , x) 0 0 0
0 0
v(τ , s) 0 dsdτ = . w(τ , s) δ(t , x)
0 1
If δ(0, x) ̸= 0 or δ(t , 0) ̸= 0, the considered problem has no continuous solution. If δ(t , x) ≡ 0, the unique solution is the trivial one. But if δ(t , x) = ε sin t +ε x we conclude that wδ (t , x) = − 1ε sin t +ε x . We will now present some auxiliary results that are necessary for the theorem we want to prove. Theorem 1 ([9]). The system of second kind integral equations u(t , x) +
t
∫
K1 (t , x, τ , x)u(τ , x)dτ + 0
x
∫
K2 (t , x, t , s)u(t , s)ds + 0
∫ t∫ 0
x
K3 (t , x, τ , s)u(τ , s)dsdτ = f (t , x),
(4)
0
where K1 (·), K2 (·), K3 (·) are (n × n) matrices with continuous elements, f (t , x) is an n-dimensional vector function with continuous elements, has a unique solution. Lemma 1 (See, for Example [2]). Let the sequence ψi+1,j+1 be defined for i = 0, 1, . . . , N − 1 and j = 0, 1, . . . , M − 1, h, q > 0, α, β, γ , l ≥ 0 and
ψi+1,j+1 ≤ l + α h
i −
ψk,j+1 + β q
k=1
j −
ψi+1,p + γ hq
p=1
j i − −
ψk,p .
k=1 p=1
Then the following inequality holds true
ψi+1,j+1 ≤ l exp(α ti + β xj )J0 2I (αβ + γ )ti xj , where ti = ih, xj = jq, h = b/N , q = a/M , J0 (·) is the zero order Bessel function and I stands for the unity of imaginary numbers. Example 5. This example shows that the regularity of the matrix pencil λA(t , x) + K (t , x) does not assure the existence and uniqueness of solution of a IAE, even in the one-dimensional case. Consider the system of one-dimensional equations of the form
1 0
∫ t v(t ) 0 + w(t ) 1 0
t 0
r
2t − τ
v(τ ) 0 dτ = , w(τ ) 0
t ∈ [0, 1],
where r is a scalar. It is easy to verify that for r ̸= 0 the matrix pencil λA(t ) + K (t , t ) is regular, though for r = 2 the given t equation has an infinite set of solutions of the form v(t ) = −t w(t ) − 2 0 w(τ )dτ , where w(t ) is an arbitrary function. Hence, even in the one-dimensional case the regularity (or singularity) of a given matrix pencil λA(t ) + K (t , t ) does not give any information on the existence of a unique solution. Now we will define a condition, independent of the regularity of a matrix pencil, which can be used in the existence and uniqueness theory for TIAE.
M.V. Bulatov, P.M. Lima / Journal of Computational and Applied Mathematics 236 (2011) 132–140
135
Definition 1 (See, for Example, [10]). We say that the matrix pencil λA(t , x) + K (t , x) satisfies the rank-power criterion (or has index one, or has a simple structure) in the domain Ω , if rank A(t , x) = deg det(λA(t , x) + K (t , x)) = k = const, or rank C (t , x) = rank C 2 (t , x) = k = const
∀(t , x) ∈ Ω ,
where C (t , x) = (λA(t , x) + K (t , x)) A(t , x) and deg(·) denotes the degree of a function as a polynomial in λ. If the rank of matrix A(t , x) is constant and equal to k on the domain Ω , it is known [10] that there exists a matrix P (t , x), nonsingular for any (t , x) ∈ Ω , whose elements have the same degree of smoothness as the elements of A(t , x), and which satisfies −1
P (t , x)A(t , x) =
B(t , x) , 0
(5)
where B(t , x) has dimension (k × n), 0 is a null matrix of dimension ((n − k) × n) and rank B(t , x) = k = const∀(t , x) ∈ Ω . Let us denote P (t , x)K (t , x) =
S ( t , x) , Q (t , x)
(6)
where K (t , x) is the same matrix as in Definition 1, P (t , x) is defined by (4), S (t , x) and Q (t , x) are matrices with dimension (k × n) and (n − k × n), respectively. Lemma 2. Let A(t , x) have the form moreover, let K (t , x) =
S (t , x) Q (t , x)
B(t , x) 0
, where B(t , x) is a (k × n) matrix and rank B(t , x) = k = const∀(t , x) ∈ Ω ;
, where S (t , x) and Q (t , x) are (k × n) and (n − k × n) matrices, respectively and the pencil λA(t , x) + K (t , x) satisfies the rank-power criterion. Then B(t , x) det ̸= 0 ∀(t , x) ∈ Ω . Q ( t , x) The proof of this lemma follows from [11]. Let us now formulate the Theorem about existence and uniqueness of solution of a system of the form (1) satisfying condition (2). Theorem 2. Assume that problem (1) satisfies the following conditions: 1. The elements of the matrices A(t , x), K (t , x, τ , s) and f (t , x) are continuous functions, as well as their first order partial derivatives and second order mixed derivatives with respect to t and x. 2. rank A(t , 0) = rank(A(t , 0)|f (t , 0)), rank A(0, x) = rank(A(0, x)|f (0, x)); 3. The matrix pencil λA(t , x) + K (t , x, t , x) satisfies the rank-power criterion in Ω . Then the considered system has a unique solution in Ω . Proof. Substituting into (1) x = 0 or t = 0, we obtain A(t , 0)u(t , 0) = f (t , 0) or A(0, x)u(0, x) = f (0, x), respectively. The solvability of the corresponding systems is guaranteed by condition 2. Multiplying both sides of Eq. (1) by P (t , x) (the matrix considered in (4)) and writing the result in the block form, we obtain
B(t , x) u(t , x) + 0
S (t , x, τ , s) Q (t , x, τ , s)
∫ t∫ x 0
S (t , x, τ , s) u(τ , s)dsdτ = Q (t , x, τ , s)
0
ϕ(t , x) , ψ(t , x)
(7)
= P (t , x)K (t , x, τ , s), (ϕ T (t , x), ψ T (t , x))T = P (t , x)f (t , x). Let us remind that B(·), S (·) and Q (·) are, respectively, (k × n), (k × n) and ((n − k) × n) matrices. Differentiating the second block line of (6) with respect to t and x, where
which is possible according to condition 1, we obtain Q (t , x, t , x)u(t , x) +
x
∫
Q x (t , x, t , s)u(t , s)ds + ′
0
∫ t∫
t
Q ′ t (t , x, τ , x)u(τ , x)dτ
0
x
Q ′′ tx (t , x, τ , s)u(τ , s)dsdτ = ψ ′′ tx (t , x).
+ 0
∫
0
Combining this last equality with the first line of (6) we obtain
B(t , x) u(t , x) + Q (t , x)
∫ t∫ x + 0
0
∫ x
0 ds + Q x (t , x, t , s)u(t , s)
∫ t
0
S (t , x, τ , s) u(τ , s)dsdτ = Q ′′ tx (t , x, τ , s)
0
0
ϕ(t , x) . ′′ ψ tx (t , x)
Q t (t , x, τ , x)u(τ , x) ′
′
ds
136
M.V. Bulatov, P.M. Lima / Journal of Computational and Applied Mathematics 236 (2011) 132–140
From the second condition of the theorem, taking into consideration Lemma 2, we conclude that the matrix nonsingular for (t , x) ∈ Ω . Finally, multiplying both sides of (7) by
B(t , x) Q (t , x)
−1
B(t , x) Q (t , x)
is
, we obtain a system of the form (3); this
system satisfies the conditions of Theorem 1 and hence has a unique solution. This concludes the proof of the theorem.
3. Numerical algorithm We will now describe and analyze a numerical algorithm for the solution of problem (1), assuming that this system satisfies the conditions of Theorem 2. This method is based on the right rectangles cubature formula (implicit Euler Method) and is analogous to the method introduced in [2,3] for the solution of two-dimensional first kind Volterra integral equations (which are a particular case of TIAE, as remarked above). We have decided to start with such an elementary method because very little is known about this class of equations. Moreover, it is very difficult to analyze the stability and convergence of higher order methods when applied to them. Numerical evidence suggests that they are often unstable. Another advantage of the considered method (right rectangles rule) is that it has a regularization effect with respect to rounding-off errors (if the stepsize is correctly chosen). For simplicity, we will illustrate this in the one-dimensional case. Consider for example the problem of computing u from its primitive f : t
∫
u(τ )dτ = f (t ),
t ∈ [0, 1], f (0) = 0, u(t ) = f ′ (t ).
(8)
0
Assume that the right-hand side function f is replaced by g, such that |g (t ) − f (t )| = δ(t ) ≤ δ . Suppose that t
∫
v(τ )dτ = g (t ),
t ∈ [0, 1].
(9)
0
Note that the considered problem is ill-conditioned, since the solution of (9) may be very different from the one of (8), even for small δ (for example, if f (t ) = 0 and g (t ) = δ sin(t /δ 2 )). However, if the perturbed Eq. (9) is solved by the right rectangles formula (analogous to the method, considered in the present article), omitting the algebraic manipulation, we obtain:
|vi − u(ti )| ≤ |(f (ti ) − f (ti−1 ))/h| + 2δ/h ≤ Lh + 2δ/h, where L < ∞. From the last inequality we conclude that the optimal precision (with minimal effect of the input data error) is obtained when h is of the order of δ 1/2 . The fact that the discretization procedure for integral equations of the first kind has a regularization effect, when there is a certain relationship between δ and h, was remarked for the first time in [12]. In the two-dimensional case, for equations that satisfy k(t , x, t , x) ̸= 0, with (t , x) ∈ Ω , a similar result can be proved; in this case, the optimal order for h and q is δ 1/3 [2]. Returning to the construction of a numerical method for TIAE, let us define in the domain Ω the mesh ti = ih, i = 1, 2, . . . , N , h = b/N , xj = jq, j = 1, 2, . . . , M , q = a/M and denote Aij = A(ti , xj ), fij = f (ti , xj ), Kijlm = K (ti , xj , tl , xm ), l ≤ i, m ≤ j; let uij be the approximate value of u(ti , xj ). The numerical method for Eq. (1), based on the right rectangles rule, has the form Aij uij + hq
j i − −
Kijlm ulm = fij .
(10)
l=1 m=1
The main result about the convergence of the described method is given by the following theorem. Theorem 3. Let the following conditions be satisfied: 1. The elements of matrices A(t , x), K (t , x, τ , s) and the vector function f (t , x) are continuous functions, as well as their partial and mixed derivatives, up to the third order; 2. The second and third conditions of Theorem 2 are satisfied. Then the following estimate is true: max ‖uij − u(ti , xj )‖ = O(h + q). i,j
Proof. Replacing in (10) ui,j by the exact value u(ti , xj ), we obtain Aij u(ti , xj ) + hq
j i − − l=1 m=1
Kijlm u(tl , xm ) = fij + Rij ,
(11)
M.V. Bulatov, P.M. Lima / Journal of Computational and Applied Mathematics 236 (2011) 132–140
137
where Rij is the error of the cubature formula of the right rectangles: Rij =
j ∫ i − −
tl
∫
(K (ti , xj , τ , s)u(τ , s) − Kijlm ulm )dsdτ .
(12)
xm−1
tl−1
l=1 m=1
xm
Denoting εij = u(ti , xj ) − uij and subtracting (10) from (11), we conclude that the error εij satisfies the following equation j i − −
Aij εij + hq
Kijlm εlm = Rij .
(13)
l=1 m=1
If we multiply both sides of (13) by Pij = P (ti , xj ) (the same matrix as in (4)), we obtain
Bij 0
εij + hq
j i − − S l=1 m=1
ijlm
Qijlm
εlm = Pij Rij .
(14)
Let us rewrite these equalities in the form of two systems Bij εij + hq
j i − −
Sijlm εlm = ρij ,
(15)
l=1 m=1
hq
j i − −
Qijlm εlm = δij ,
(16)
l=1 m=1
where (ρijT , δijT )T = Pij Rij . Let us introduce the difference operators ∆h and ∆q , which are defined by the equalities
∆h g (ti , xj ) = ∆q g (ti , xj ) =
1 h 1 q
(g (ti , xj ) − g (ti−1 , xj )),
i = 2, 3, . . . , N ,
(g (ti , xj ) − g (ti , xj−1 )),
j = 2, 3, . . . , M ,
for any given vector function g (t , x). These are the simplest difference analogs of the partial derivatives with respect to t and to x, respectively. Applying successively these two operators to both sides of (16), we obtain Qijij εij +
j −1 i−1 − − (Qijlj − Qi−1jlj )εlj + (Qijim − Qij−1im )εim l =1
m=1
i −1 j −1
+
−−
(Qijlm − Qi−1jlm − Qij−1lm + Qi−1j−1lm )εlm =
l=1 m=1
1 h·q
(δij − δij−1 − δi−1j + δi−1j−1 ).
Combining this system with (15) we obtain
Bij + hqSijij Qijij
+
j −1 i −1 − − l=1 m=1
εij +
i −1 − l =1
Qijlj
0 − Qi−1jlj
εlj +
j −1 − m=1
hqSijlm Qijlm − Qi−1jlm − Qij−1lm + Qi−1j−1lm
Qijim
0 − Qij−1im
εlm = 1 hq
εim ρij
, (δij − δij−1 − δi−1j + δi−1j−1 )
(17)
where 0 denotes a (k × n) null matrix. From the third condition of Theorem 2 and Lemma 2 it follows that for sufficiently small h and q the matrix
Bij + hqSijij Qijij
is nonsingular, that is,
−1 Bij + hqSijij ≤ K1 < ∞. Qijij
(18)
Taking into consideration the first condition of the theorem, we conclude that the differences on the left-hand side of (17) satisfy
‖Qijlj − Qi−1jlj ‖ ≤ hK2 ,
K2 < ∞,
(19)
138
M.V. Bulatov, P.M. Lima / Journal of Computational and Applied Mathematics 236 (2011) 132–140
‖Qijim − Qij−1im ‖ ≤ qK3 ,
K3 < ∞,
(20)
‖Qijlm − Qi−1jlm − Qij−1lm − Qi−1j−1lm ‖ ≤ (h + q)K4 ,
K4 < ∞.
(21)
Moreover, from the error estimate for the right rectangles cubature formula it follows that
‖ρij ‖ ≤ (h + q)K5 ,
K5 < ∞.
(22)
According to formula (12) and the first condition of the theorem, after some cumbersome algebraic manipulations (which we omit here) we conclude that 1 hq
‖δij − δij−1 − δi−1j + δi−1j−1 ‖ ≤ (h + q)K6 ,
K6 < ∞.
(23)
Finally, introducing in (17) the estimates (18)–(23), we obtain
‖εij ‖ ≤
K6 K1
(h + q) +
j −1 i−1 i−1 j−1 K2 − K3 − K4 − − ‖εim ‖ + hq ‖εlm ‖. h ‖εlj ‖ + q K1 l = 1 K1 m=1 K1 l=1 m=1
By Lemma 1, from this last inequality we conclude that
‖εij ‖ = O(h + q). This concludes the proof of the theorem.
4. Numerical examples Example 1. Consider the integro-algebraic system A(t , x)u(t , x) +
∫ t∫ 0
x
K (t , x, τ , s)u(τ , s)dsdτ = f (t , x),
(24)
0
where
[
1 A= 0
[ K =
0 1
]
0 , 0
]
1 . 1
The right-hand side function is given by: f 1 ( t , x) = 1 +
x2 t 2 4
,
f2 (t , x) = xt +
x2 t 2 4
.
Let us check that this system satisfies the conditions of Theorem 2. 1. A, K and f are smooth functions of t and x. 2. rank A(t , 0) = rank(A(t , 0)|f (t , 0)) = 1; rank(A(0, x)) = rank(A(0, x)|f (0, x)) = 1. 3. We have
λA + K =
[
λ 1
]
1 1
and therefore det(λA(t , x) + K (t , x, t , x)) = λ − 1. Hence the degree of det(λA + K ) is 1, for t , x > 0, the same as rank A. Since the conditions of Theorem 2 are satisfied, the considered system has a unique continuous solution on any domain
Ω = {(t , x) : 0 ≤ t ≤ a, 0 ≤ x ≤ b}. By direct substitution we can check that this solution is u(t , x) = (v(t , x), w(t , x))T = (1, tx)T . We shall consider in D a uniform mesh, with stepsize h on t and x: ti = ih,
i = 1, . . . , N ,
xj = jh,
j = 1, . . . , M ,
where h = a/N = b/M. With these notations we have applied the numerical method (10). The numerical results are displayed in Table 1. We have considered a = b = 1 and used five different meshes with h = q = 1/4, 1/8, 1/16, 1/32, 1/64. The error norm, given by
‖ϵh,q ‖∞ =
max
i=1,...,N ,j=1,...,M
‖ui,j − u(ti , xj )‖∞ ,
(25)
M.V. Bulatov, P.M. Lima / Journal of Computational and Applied Mathematics 236 (2011) 132–140
139
Table 1 Numerical results for Example 1. h=q
1/4
1/8
1/16
1/32
1/64
ϵhq
0.234375 0.953
0.121094 0.977
0.0615234 0.986
0.0310059 0.994
0.0155064
p
Table 2 Numerical results for Example 2. h=q
1/8
1/16
1/32
1/64
1/128
ϵhq
0.0314114 0.864
0.0172628 0.937
0.00901509 0.970
0.00460285 0.985
0.00232518
p
is given for each different mesh. We also display the estimates of the convergence order: p = log2
‖ϵh,q ‖∞ . ‖ϵh/2,q/2 ‖∞
(26)
The numerical results are in good agreement with the theoretical convergence order, given by Theorem 3. Example 2. Now let A(t , x) =
[
1 ψ(t , x)
K (t , x, τ , s) =
] ψ(t , x) 2 , ψ(t , x)
t +x−τ +s 1
[
(27)
]
0 , 5 exp(t − x + τ + s)
(28)
ψ(t , x) = t + x. In this case the right-hand side function is given by f ( t , x) =
f1 (t , x) 1 + t 2 − x + t 2 x/2 + t 4 x/12 + 3tx2 /2 − t 2 x2 /4 + t 3 x2 /2 − 5tx3 /6 + e−t −x (t + x), = f2 (t , x) tx + 5 exp(t − x)tx + t 3 x/3 − tx2 /2 + (1 + t 2 − x)(t + x) + exp(−t − x)(t + x)2
[
]
[
]
and the exact solution is u(t , x) = (v(t , x), w(t , x))T = (1 + t 2 − x, exp(−t − x))T ,
(29)
as can be verified, by direct substitution. Let us check that this system satisfies the conditions of Theorem 2. 1. A, K and f are smooth functions of t and x. 2. rank A(t , 0) = rank(A(t , 0)|f (t , 0)) = 1; rank(A(0, x)) = rank(A(0, x)|f (0, x)) = 1. 3. We have
] λ(t + x) λ(t + x)2 + 5 exp(2t ) and therefore det(λA(t , x) + K (t , x, t , x)) = λ 5 exp(2t ) + 2x(t + x)2 − (t + x) + 10x exp(2t ). Hence the degree of det(λA(t , x) + K (t , x, t , x)) is 1, for t , x > 0, the same as rank A. λA(t , x) + K (t , x, t , x) =
[
λ + 2x 1 + λ(t + x)
Since the conditions of Theorem 2 are satisfied, we may assure that (29) is the unique solution of the considered system on any domain Ω = {(t , x) : 0 ≤ t ≤ a, 0 ≤ x ≤ b}. The numerical results for this example are displayed in Table 2. As before, we have considered a = b = 1. We have used five different meshes, this time with h = q = 1/8, 1/16, 1/32, 1/64, 1/128. The error norm ϵhq is given by (25) and the empirical estimate of the convergence order p is given by (26). As in the first example, the numerical results confirm the theoretical convergence order. Example 3. Now let A(t , x) =
[
1 ψ(t , x)
K (t , x, τ , s) =
] ψ(t , x) 2 , ψ(t , x)
t +x+τ −s τ + 2s
[
ψ(t , x) = t − x.
(30)
]
0 , 5 exp(t − x − τ + s)
(31)
140
M.V. Bulatov, P.M. Lima / Journal of Computational and Applied Mathematics 236 (2011) 132–140 Table 3 Numerical results for Example 3. h=q
1/8
1/16
1/32
1/64
1/128
ϵhq
0.626502 0.752
0.0371996 0.881
0.020194 0.942
0.0105117 0.971
0.00536162
p
In this case the right-hand side function is given by f (t , x) =
=
f 1 ( t , x) f 2 ( t , x)
[
]
1 + t + e−t −x (t − x) + 3t 2 x/2 + 5t 3 x/6 − x2 + tx2 /2 + t 2 x2 /4 − t 2 x3 /2 − tx4 /12, tx exp(−t − x)(t − x)2 + (t − x)(1 + t − x2 ) + t (3 + 2t ) + 3(2 + t )x − tx2 − 3x3 + 5 exp(−x)x sinh(t ) 6
and the exact solution is u(t , x) = (v(t , x), w(t , x))T = (1 + t − x2 , exp(−t − x))T ,
(32)
as can be verified, by direct substitution. Let us check that this system satisfies the conditions of Theorem 2. 1. A, K and f are smooth functions of t and x. 2. rank A(t , 0) = rank(A(t , 0)|f (t , 0)) = 1; rank(A(0, x)) = rank(A(0, x)|f (0, x)) = 1. 3. We have
λA(t , x) + K (t , x, t , x) =
[
λ + 2t λ(t − x) + t + 2x
λ(t − x) λ(t − x)2 + 5
]
and therefore det(λA(t , x) + K (t , x, t , x)) = λ(5 + 2t 3 + 2x2 + tx(−1 + 2x) − t 2 (1 + 4x)) + 10t. Hence the degree of det(λA(t , x) + K (t , x, t , x)) is 1, for t , x > 0, the same as rank A. Since the conditions of Theorem 2 are satisfied, we may assure that (32) is the unique solution of the considered system on any domain Ω = {(t , x) : 0 ≤ t ≤ a, 0 ≤ x ≤ b}. The numerical results for this example are displayed in Table 3. As in the previous example, we have considered a = b = 1 and used five different meshes, with h = q = 1/8, 1/16, 1/32, 1/64, 1/128. The error norm ϵhq is given again by (25) and the estimated convergence order p is defined by (26). Once again the numerical results are in agreement with Theorem 3. The numerical method (10) was coded as a program in Mathematica. Floating point was used with the default precision (16 digits). All the computations were carried out in a personal computer with a Pentium 2.3 GHz processor. The obtained results confirm that the computational method (10) provides a simple and stable algorithm to approximate the solution of system (1). Its low convergence order is in some sense compensated by the simplicity of computations. We cannot provide here a comparison of methods, since (as far as we know) there were no previous attempts to solve such problems numerically. We believe that certain higher order methods can be applicable. In particular, it seems reasonable to apply a method based on the trapezoidal rule. However, the stability analysis is much more difficult (as remarked above) for higher order methods. This subject deserves further analysis and we leave it as future work. Acknowledgments M. Bulatov acknowledges financial support from RFBR, grant 10-01-00571a; P. Lima acknowledges support from FCT, grant PTDC/MAT/101867/2008. References [1] G.V. Demidenko, S.V. Uspenskii, Partial Differential Equations and Systems not Solvable with Respect to the Highest-Order Derivative, Taylor & Francis, 2005. [2] Ten Men Yan, Approximate solution of linear volterra integral equations of the first kind, Ph.D. Dissertation, Irkutsk, 1985 (in Russian). [3] S. McKee, Tao Tang, T. Diogo, An Euler-type method for two-dimensional Volterra integral equations of the first kind, IMA J. Numer. Anal. 20 (2000) 423–440. [4] H. Brunner, Collocation Methods for Volterra Integral and Related Functional Equations, University Press, Cambridge, 2004. [5] M. Hadizadeh, F. Ghoreishi, S. Pishbin, Jacobi spectral solution for integral algebraic equations of index-2, Appl. Numer. Math. 61 (2011) 131–148. [6] J.-P. Kauthen, The numerical solution of integral-algebraic equations of index-1 by pollinomial spline collocation methods, Math. Comp. 70 (2001) 1503–1514. [7] M.V. Bulatov, V.F. Chistyakov, The properties of differential–algebraic systems and their integral analogs, September, Memorial University of Newfoundland, 1997, Preprint. [8] M.V. Bulatov, Regularization of singular systems of Volterra integral equations, Comput. Math. Math. Phys. 42 (2002) 315–320. [9] V.I. Smirnov, A Course of Higher Mathematics, vol.4, Pergamon Press, Oxford, 1964. [10] V.F. Chistyakov, Differential-Algebraic Operators with Finite-Dimensional Kernel, Nauka, Novosibirsk, 1996 (in Russian). [11] M.V. Bulatov, Transformation of differential–algebraic systems of equations, Comput. Math. Math. Phys. 34 (1994) 301–311. [12] A.S. Apartsyn, A.B. Bakushinsky, Approximate solution of Volterra integral equations of the first kind by the quadratures method, Differential and Integral Equations, V. 1, Irkustk, 1972, pp. 248–258 (in Russian).