- Email: [email protected]
A fast numerical solution method for two dimensional Fredholm integral equations of the second kind
Applied Numerical Mathematics 59 (2009) 1709–1719
Contents lists available at ScienceDirect
Applied Numerical Mathematics www.elsevier.com/locate/apnum
A fast numerical solution method for two dimensional Fredholm integral equations of the second kind ✩ Wen-Jing Xie, Fu-Rong Lin ∗ Department of Mathematics, Shantou University, Shantou Guangdong 515063, China
a r t i c l e
i n f o
Article history: Received 23 August 2006 Received in revised form 3 July 2008 Accepted 30 January 2009 Available online 7 February 2009 MSC: 45L10 65R20 Keywords: Integral equation Polynomial interpolation Approximate matrix Residual correction scheme
a b s t r a c t In this paper we consider numerical solution methods for two dimensional Fredholm integral equation of the second kind
1 1 f ( x, y ) −
a(x, y , u , v ) f (u , v ) du dv = g (x, y ),
(x, y ) ∈ [−1, 1] × [−1, 1],
−1 −1
where a(x, y , u , v ) is smooth and g (x, y ) is in L 2 [−1, 1]2 . We discuss polynomial interpolation methods for four-variable functions and then use the interpolating polynomial to approximate the kernel function a(x, y , u , v ). Based on the approximation we deduce fast matrix-vector multiplication algorithms and efficient preconditioners for the above two dimensional integral equations. The residual correction scheme is used to solve the discretization linear system. © 2009 IMACS. Published by Elsevier B.V. All rights reserved.
1. Introduction Many problems in engineering and mechanics can be transformed into two dimensional Fredholm integral equations of the second kind. For example, it is usually required to solve Fredholm integral equations in the calculation of plasma physics [11]. There are many works on developing and analyzing numerical methods for solving Fredholm integral equations of the second kind (see [1–4,6,10]). In recent years, a number of algorithms for the fast numerical solution of integral equations of the second kind have been developed; see, e.g., [1,8,13,14,18–20]. The fast multipole method proposed in [13] combines the use of low-degree polynomial interpolation of the kernel functions with a divide-and-conquer strategy. For kernel functions that are Coulombic or gravitational in nature, it results in an order O( N ) algorithm for the matrix-vector multiplications, where N is the number of quadrature points used in the discretization. In [8], an O( N log N ) algorithm was developed by exploiting the connections between the use of wavelets and their applications on Calderon–Zygmund operators. In [18], we proposed an approximation scheme to obtain low rank approximation of the discretization matrix. By using preconditioned iterative methods such as the residual correction (RC) scheme or the preconditioned conjugate gradient (PCG) method, the total cost for solving the integral equation with smooth kernel is O( N ) operations (ops). The papers [14,19,20] are concerned with fast algorithms for two dimensional integral equations with weakly singular kernel functions. In [14] and [19], different versions of fast multipole methods are proposed to solve integral equations with weakly singular kernels. In [20], the authors developed a fast wavelet collocation method for the integral equations defined ✩
*
The research was supported by NSF of China No. 10271070. Corresponding author. E-mail address: [email protected] (F.-R. Lin).
0168-9274/$30.00 © 2009 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.apnum.2009.01.009
1710
W.-J. Xie, F.-R. Lin / Applied Numerical Mathematics 59 (2009) 1709–1719
on polygons. Besides fast algorithms, there are many investigations on improving the accuracy of numerical solutions, see for instance [15] and [21]. In [15], a Nyström method for two dimensional Fredholm integral equation of second kind was used. The asymptotic expansion of numerical solution is obtained and then the Richardson extrapolation method is used to enhance the precision of the numerical solution. In [21], based on the finite-element solution, an iterative correction method was proposed. It shows that if the kernel function is smooth, the precision could be improved considerably. In this paper, we consider the numerical solution of two dimensional Fredholm integral equation of the second kind with smooth kernel:
1 1 f (x, y ) −
a(x, y , u , v ) f (u , v ) du dv = g (x, y ),
(x, y ) ∈ [−1, 1] × [−1, 1],
(1.1)
−1 −1
where the kernel function a(x, y , u , v ) is smooth in [−1, 1]4 , the unknown function f (x, y ) and the right-hand side function g (x, y ) are in L 2 [−1, 1]2 . For the case where the integration domain is not [−1, 1]2 , say [α , β]2 , we can use the linear transformation
x = (2x − α − β)/(β − α ), y = (2 y − α − β)/(β − α )
to transform the domain into [−1, 1]2 . (N ) (N ) (N ) We discretize the equation by using certain numerical quadrature. Let −1 t 1 < t 2 < · · · < t N 1 be the quadrature (N )
points and w i
, i = 1, 2, . . . , N, be the corresponding weights, a numerical quadrature is defined as
1 f (t ) dt ≈
N
(N )
wi
(N )
f ti
(1.2)
.
i =1
−1
Typical examples include Gaussian rules, repeated Gaussian rules, and repeated Newton–Cotes rules, etc., see for instance [10]. Extending the above formula to two dimensional case, we get the following approximate system for (1.1): f (x, y ) −
N N
(N )
wi
(N )
(N )
w j a x, y , t i
, t (jN ) f t i(N ) , t (jN ) = g (x, y ),
x, y ∈ [−1, 1].
i =1 j =1
In particular, we get the discretization linear system of (1.1):
˜ )f = g, (I − A W
(1.3) (N )
where f and g are vectors obtained by reordering the matrices [ f (t i
(N )
(N )
, t j )]iN, j =1 and [ g (t i
(N )
, t j )]iN, j=1 row-by-row respec-
˜ = W ⊗ W with W the weight matrix diag( w , w , . . . , w ), and tively, I is the identity matrix, W N 1 2 (N )
⎛
A (1,1)
⎜ (2,1) ⎜ A ⎜ A=⎜ . ⎜ .. ⎝ A ( N , 1) (N )
with A (i , j ) = [a(t i
···
A (1, N )
A (2,2)
A (2, N ) ⎟
.. .
··· .. .
A ( N , 2)
···
A (N ,N )
(N )
, tl
(N )
(N )
⎞
A (1,2)
.. .
(N )
⎟ ⎟ ⎟ ⎟ ⎠
(1.4)
(N )
, t j , tm )]lN,m=1 . Here ⊗ denotes the Kronecker tensor product.
To solve the linear system (1.3) by classical direct methods such as Gaussian elimination method requires O( N 6 ) operations (ops). For iterative methods such as the conjugate gradient (CG) method (see [12]), each iteration requires O( N 4 ) ops. Therefore even for well-conditioned problems, the CG method requires O( N 4 ) ops, which for large-scale problems is often prohibitive. In this paper, we propose a fast numerical solution method for the linear system (1.3) by using the polynomial interpolation technique. That is, rather than solving the discretization linear system (1.3), we solve a Nyström approximation of the following equation
f
[n]
1 1 (x, y ) −
an (x, y , u , v ) f [n] (u , v ) du dv = g (x, y ),
−1 x, y 1,
(1.5)
−1 −1
where an (x, y , u , v ) is an interpolation based approximation to a(x, y , u , v ) in which n Chebyshev nodes are used in each space variable. More precisely, we solve the following approximation equation of (1.3):
˜ )f[n] = g, ( I − Aa W
(1.6)
W.-J. Xie, F.-R. Lin / Applied Numerical Mathematics 59 (2009) 1709–1719
(N )
(N )
(N )
1711
(N )
where A a = ([ A a ](i , j ) )iN, j =1 with [ A a ](i , j ) = [an (t i , tl , t j , tm )]lN,m=1 . The outline of this paper is as follows. In Section 2, we discuss how to interpolate multi-variable functions efficiently and analyze the error in the interpolating polynomials. In Section 3, we propose fast algorithms for solving (1.6) by using the results obtained in Section 2 and discuss the complexity of our algorithms. Numerical results are given in Section 4 to illustrate the efficiency and the accuracy of our algorithms. Finally, concluding remarks are given in Section 5. 2. Interpolation of multi-variable functions In this section we discuss the approximation of two-variable and four-variable smooth functions h(x, y ) (−1 x, y 1) and a(x, y , u , v ) (−1 x, y , u , v 1) by polynomial interpolation at Chebyshev grids. We discuss how to obtain the interpolating polynomials efficiently and analyze the error in the approximation. 2.1. Interpolation (n)
Let c j (x), j = 0, 1, . . . , n − 1, be the scaled Chebyshev polynomials defined by
(n)
c j (x) = δi cos j arccos(x) ,
√
Here δ0 = 1/ n and δ j =
j = 1, 2, . . . , n − 1, x ∈ [−1, 1].
√
2δ0 for j = 1, . . . , n − 1. Let (2i − 1)π (n) p i = cos , i = 1, 2, . . . , n,
2n
the roots of the degree-n Chebyshev polynomial cos(n arccos(x)), we call
(n)
(n)
pi , p j
i , j = 1, 2, . . . , n
the n × n Chebyshev grid of the domain [−1, 1]2 . Let the interpolating polynomial of h(x, y ) at the Chebyshev grid of [−1, 1]2 be denoted by (C n h)(x, y ), we have
(Cn h)(x, y ) =
n n p =1 q=1
(n)
(n)
λ pq c p −1 (x)cq−1 ( y ),
−1 x, y 1,
(2.1)
where the coefficients λ pq are determined by the interpolating conditions:
(n)
(n)
h pi , p j
=
n n p =1 q=1
λ pq c (pn−) 1 p (i n) cq(n−)1 p (jn) ,
i , j = 1, 2, . . . , n.
(2.2)
(n) (n) Let Λ = [λi j ]ni, j =1 and c(n) (·) = (c 0 (·), . . . , cn−1 (·)) T , the formula (2.1) can be written as
(Cn h)(x, y ) = c(n) (x)T Λc(n) ( y ). (n)
(n)
Let H g = [h( p i , p j )]ni, j =1 , and
(n) n i , j =1
C = c i −1 p j
n (i − 1)(2 j − 1)π = δi −1 cos 2n
,
i , j =1
the conditions (2.2) can be written in compact form as H g = C T ΛC . Note that C is the cosine transform matrix, which is orthogonal, we have Λ = C H g C T and it follows that Λ can be obtained in O(n2 log n) ops by using fast cosine transform. In summary, we get an approximation of h(x, y ): h(x, y ) ≈ (Cn h)(x, y ) =
n n
C HgCT
p =1 q=1
(n) (n) c (x)cq−1 ( y ). p ,q p −1 (n)
(n)
(n)
(n)
Similarly, the interpolating polynomial of h(x, y , u , v ) at the Chebyshev grid of ({( p i , p j , pl , pm ) | i , j , l, m = 1, . . . , n}) is given by
(Cn h)(x, y , u , v ) =
n n n n p =1 q=1 r =1 s=1
λ pqrs c (pn−) 1 (x)cq(n−)1 ( y )cr(n−)1 (u )c (sn−)1 ( v )
(2.3)
for −1 x, y , u , v 1, where the coefficients λ pqrs are determined by the interpolating conditions:
(n)
(n)
(n)
(n)
h p i , p j , pl , pm
=
n n n n p =1 q=1 r =1 s=1
(n) (n) (n) (n) (n) (n) (n) (n) λ pqrs c p −1 p i cq−1 p j cr −1 pl c s−1 pm ,
(2.4)
1712
W.-J. Xie, F.-R. Lin / Applied Numerical Mathematics 59 (2009) 1709–1719
for i , j , l, m = 1, 2, . . . , n. Similar to the two-variable case, we have
T (Cn h)(x, y , u , v ) = c(n) (x) ⊗ c(n) ( y ) Λ c(n) (u ) ⊗ c(n) ( v ) ,
where
⎛
Λ(1,1) Λ(1,2) · · · Λ(1,n)
(2.5)
⎞
⎟ ⎜ (2,1) ⎜Λ Λ(2,2) · · · Λ(2,n) ⎟ ⎟ ⎜ Λ = ⎜ .. .. .. .. ⎟ ⎟ ⎜ . . . . ⎠ ⎝ ( n , 1) ( n , 2) ( n,n) Λ Λ ··· Λ
with Λ(i , j ) = [λiljm ]ln,m=1 .
(i , j ) n (i , j ) ]i , j=1 with n × n blocks H g
Let H g be the n × n block matrix H g = [ H g get
(n)
(n)
(n)
(n)
= [h( p i , pl , p j , pm )]ln,m=1 . From (2.4) we
H g = (C ⊗ C ) T Λ(C ⊗ C ). Therefore,
Λ = (C ⊗ C ) H g (C ⊗ C ) T ,
(2.6) 4
which can be obtained by using fast cosine transform in O(n log n) ops. 2.2. Error analysis In this subsection, we discuss the error in the interpolating polynomial of a(x, y , u , v ) at the Chebyshev grid ((Cn a)(x, y , u , v )). We first discuss the error in the interpolating polynomials for one-variable functions. Lemma 1 (Jackson’s Fifth Theorem, see for instance [16, p. 95]). Suppose f ∈ C [−1, 1] and f (k) ∈ C [−1, 1] for some positive integer k. Let (Bn f )(x) be the best approximation polynomial of degree-(n − 1) for f (x) ∈ C [−1, 1], i.e.,
Bn f − f =
min
qn :degree-(n−1) polynomial
qn − f .
Here and throughout the paper, f = maxx∈[−1,1] | f (x)| is the supremum norm. Then for n k, we have
E n ( f ) ≡ Bn f − f
π
k
2
(k) f / (n + 1)n · · · (n − k + 2)
π
k
2
(k) f /(n − k + 2)k .
The Lebesgue function and Lebesgue constant are useful tools for analyzing the error of interpolating polynomials. For a given interpolating grid {x1 , . . . , xn }, the Lebesgue function is defined as L n (x) =
n hk (x), k=1
where hk (x) =
n
x−xi i =1,i =k xk −xi
is the kth Lagrange basis polynomial, and the Lebesgue constant is defined as
Ωn = Ln , see for instance [9]. For Chebyshev grids, we have (see [17])
Ωn < 2 +
2
π
From (Cn f )(x) =
log n.
n
k=1
(2.7) (n)
f ( pk )hk (x), it follows that
n (n) (Cn f )(x) f p hk (x) Ωn f .
(2.8)
k
k=1
Lemma 2. Suppose f (x) ∈ C k [−1, 1]. If n > k, then
f (x) − (Cn f )(x) δn,k f (k) = O n−k log n ,
where
δn,k =
π 2
k 3+
2
π
x ∈ [−1, 1],
log n (n − k + 2)−k .
(2.9)
W.-J. Xie, F.-R. Lin / Applied Numerical Mathematics 59 (2009) 1709–1719
1713
Proof. It is easily seen that Cn f is a linear operator and Cn Bn f = Bn f . Thus
(Cn f )(x) − f (x) f (x) − (Bn f )(x) + (Cn Bn f )(x) − (Cn f )(x) E n ( f ) + Cn (Bn f − f ) (x) E n ( f ) + Ωn Bn f − f by (2.8) = (1 + Ωn ) E n ( f ). Using (2.7) and the result of Lemma 1, the result of the lemma follows from the above inequality.
2
From the above lemma, we see that Chebyshev grids are very powerful in high degree polynomial approximation. In [7], continuous functions defined in [−1, 1] are approximated very accurately by using the polynomial interpolant with sufficiently many Chebyshev points. We now extend the above result to interpolating polynomials of two-variable functions. Lemma 3. Suppose the function a(u , v ) is such that ∂ k a/∂ uk and ∂ k a/∂ v k are continuous in the domain [−1, 1]2 , where k is a positive integer. Let Cn a be the degree-(n − 1) interpolating polynomial of a(u , v ) at the Chebyshev grid (cf. (2.1)). Then
a − Cn a
π
k
2
3+
2
π
2 log n
(n − k + 2)−k c = O n−k log2 n
for n > k, where c = max{∂ k a/∂ uk , ∂ k a/∂ v k }. Proof. First of all, we define the one-variable interpolating operators for two-variable function: L v a (u , v ) =
n
(n) (n)
a u, pi
hi ( v )
(2.10)
i =1
and L u a (u , v ) =
n
(n)
(n)
a p i , v h i (u ),
i =1
(n)
where h i
(n)
(n)
(n)
is the ith Lagrange basis polynomial at the Chebyshev grid { p 1 , p 2 , . . . , pn }. Obviously, the interpolating (n)
(n)
polynomial (Cn a)(u , v ) of a(u , v ) at the Chebyshev grid {( p i , p j ) | i , j = 1, 2, . . . , n} is given by
(Cn a)(u , v ) =
n n
(n)
(n) (n)
a pi , p j
(n)
h i (u )h j (u ) = ( L u L v a)(u , v ) = ( L v L u a)(u , v ).
i =1 j =1
By the triangular inequality a − Cn a a − L v a + L v a − L u L v a, we only require to estimate a − L v a and L v a − L u L v a respectively. By using Lemma 2 and (2.10), we have
k k a(u , v ) − ( L v a)(u , v ) δn,k max ∂ a (u , t ) δn,k ∂ a , ∂ vk −1t 1 ∂ v k
Thus
−1 u , v 1.
k ∂ a a − L v a δn,k ∂ v k .
(2.11)
As for L v a − L u L v a, by using (2.8), we have
k ∂ a L v a − L u L v a = L v (a − L u a) Ωn a − L u a Ωn δn,k ∂ uk .
Let c = max{∂ k a/∂ uk , ∂ k a/∂ v k }, it follows from (2.11) and (2.12) that
a − Cn a δn,k (1 + Ωn )c . Therefore, from (2.7) and (2.9), the result follows.
2
Now we can prove our main theorem on the error of interpolating polynomials for four-variable functions.
(2.12)
1714
W.-J. Xie, F.-R. Lin / Applied Numerical Mathematics 59 (2009) 1709–1719
Theorem 1. Let a(x, y , u , v ) be a function defined in the domain [−1, 1]4 such that ∂ k a/∂ xk , ∂ k a/∂ yk , ∂ k a/∂ uk and ∂ k a/∂ v k are continuous in the domain [−1, 1]4 , where k is a positive integer. Let (Cn a)(x, y , u , v ) be the interpolating polynomial of a(x, y , u , v ) at the Chebyshev grid (cf. (2.3)). Then for n > k, we have
a − Cn a = O n−k log4 n .
(2.13)
Proof. Let L x , L y , L u , and L v be the one-variable interpolating operators which are defined as in the proof of Lemma 3, we have Cn a = L x L y L u L v a. Therefore
a − Cn a = a − L x L y L u L v a a − L v a + L v a − L u L v a + L u L v a − L y L u L v a + L y L u L v a − L x L y L u L v a. By using Lemma 2, we have a − L v a δn,k ∂ k a/∂ v k . Moreover, by using (2.8), we get
L v a − L u L v a L v a − L u a Ωn δn,k ∂ k a/∂ uk , L y L u L v a − L u L v a L u L v ( L y a − a) Ωn2 δn,k ∂ k a/∂ yk , L y L u L v a − L x L y L u L v a = L y L u L v (a − L x a) Ωn3 δn,k ∂ k a/∂ xk , where Ωn is the Lebesgue constant. It follows that
a − Cn a δn,k 1 + Ωn + Ωn2 + Ωn3 max ∂ k a/∂ xk , ∂ k a/∂ yk , ∂ k a/∂ uk , ∂ k a/∂ v k . Using (2.7) and (2.9), we obtain (2.13).
2
3. Fast algorithms for integral equations In this section, we propose new algorithms for the fast numerical solution of two dimensional Fredholm equations of the second kind by using the interpolation method discussed in Subsection 2.1. We approximate the matrix A defined by (1.4) with A a such that
( I − A W ˜ <
˜ ) − ( I − Aa W ˜ ) = ( A − A a ) W for given > 0, so that the solution of the approximation linear system (1.6) approximates the solution of (1.3) accurately. ˜ ) as We then approximate A is by a less accurate matrix B a and split ( I − A a W
˜ = (I − Ba W ˜ ) − ( Aa − B a )W ˜. I − Aa W Thus, Eq. (1.6) is equivalent to
˜ f[n] + g. ˜ )f[n] = ( A a − B a ) W (I − Ba W In this way, we come to the following iterative method:
˜ )f(q+1) = ( A a − B a ) W ˜ f(q) + g, (I − Ba W
q = 0, 1, 2, . . . ,
(3.1)
or equivalently
˜ )f(q) , r(q) = g − ( I − A a W
(3.2)
˜ )d(q) = r(q) , (I − Ba W
(3.3)
(q+1)
f
(q)
=f
(q)
+d .
Eq. (3.3) is called the preconditioning system. The above iterative method is often called a residual correction (RC) scheme or a method of iterative improvement. We note that the main cost per iteration of the RC scheme is the matrix-vector ˜ )f(q) and the solution of the preconditioning system (3.3), and the method converges if and only if multiplication ( I − A a W
ρ = ( I − B a W˜ )−1 ( A a − B a ) W˜ < 1. Moreover, the smaller the
ρ , the faster the convergence.
W.-J. Xie, F.-R. Lin / Applied Numerical Mathematics 59 (2009) 1709–1719
1715
3.1. The approximate matrices By (2.5), we have
(N )
a ti
(N ) (N ) , t (jN ) , tl(N ) , tm ≈ (Cn a) t i(N ) , t (jN ) , tl(N ) , tm (N ) (N ) T (n) (N ) (N ) = c(n) t i Λ c tl ⊗ c(n) t j ⊗ c(n) tm
for i , j , l, m = 1, . . . , N. It is follows that
˜ = ( P ⊗ P )T Λ( P ⊗ P ) A≈A where
(n)
(3.4)
( N ) j =1: N . i =1:n
P = c i −1 t j
(3.5)
Based on the discussion in Subsection 2.1, we have the following algorithm for computing Λ and P . Algorithm 1 (The construction of approximation matrix). (n)
Step 1. Choose n, form the Chebyshev points: p i , i = 1, 2, . . . , n and compute the function values of a(x, y , u , v ) at the Chebyshev grid of the domain [−1, 1]4 :
⎛
(1,1)
Hg
⎜ (2,1) ⎜ Hg Hg = ⎜ ⎜ .. ⎝ .
(1,2)
· · · H (g1,n)
⎞
(2,2)
··· Hg .. .. . .
⎟ ⎟ ⎟ ⎠
Hg Hg
.. .
(n,1)
Hg
(n)
(n)
Hg (i , j )
(n,2)
(2,n) ⎟
(n,n)
··· Hg
(n)
(n)
with H g = [a( p i , pl , p j , pm )]ln,m=1 . Step 2. Compute Λ by using two dimensional fast cosine transform on H g (cf. (2.6)). Step 3. Compute the matrix P which is defined in (3.5). We note that both A a and B a can be obtained by using Algorithm 1. We choose a large n (denoted by na ) for computing A a and a moderate n (denoted by nb ) for computing B a , see numerical results in Section 4. In the following, we will denote A a and B a as A a = ( P A ⊗ P A ) T Λ A ( P A ⊗ P A ) and
B a = ( P B ⊗ P B )T Λ B ( P B ⊗ P B )
(3.6)
respectively. It is easily seen that P A , P B , Λ A , and Λ B are of the size na × N, nb × N, na2 × na2 , and
˜ we have the following theorem. About the error in the approximate matrix A,
nb2
× nb2
respectively.
Theorem 2. Let a(x, y , u , v ) be a function defined in the domain [−1, 1]4 such that ∂ k a/∂ xk , ∂ k a/∂ yk , ∂ k a/∂ uk and ∂ k a/∂ v k are ˜ be the approximate matrix of A given by (3.4). continuous, where k is a positive integer. Let A (N )
(i) If all weights w i
are positive and formula (1.2) exactly holds for constant functions, then we have
( A˜ − A ) W ˜ = O n−k log4 n
for n > k, where · denotes the infinity matrix norm. (ii) If there exist two constants α and β satisfying 0 α β such that exactly holds for constant functions, then we have
α / N w (i N ) β/ N for i = 1, 2, . . . , N and formula (1.2)
( A˜ − A ) W ˜ = O n−k log4 n F
for n > k, where · F denotes the Frobenius norm. Proof. (i) By setting f (x) = 1 in (1.2) we get
N
l=1
˜ respectively. Note that A max
i , j ,l,m=1,..., N
(i , j ) (i , j ) A˜ −A = l,m
l,m
max
i , j ,l,m=1,..., N
max
(N )
wl
(i , j ) (i , j ) = 2. Let Al,m and A˜ l,m be (l, m)-entry of the (i , j )-block of A and
(N ) ( N ) (Cn a) t (N ) , t (N ) , t (N ) , tm − a t ( N ) , t ( N ) , t ( N ) , tm
−1x, y ,u , v 1
i
j
l
i
(C n a)(x, y , u , v ) − a(x, y , u , v ),
j
l
1716
W.-J. Xie, F.-R. Lin / Applied Numerical Mathematics 59 (2009) 1709–1719
by using (2.13) in Theorem 1, we get max
i , j ,l,m=1,..., N
(i , j ) (i , j ) A˜ −A = O n−k log4 n . l,m
l,m
Therefore N (i , j ) (i , j ) ( N ) ( N ) A˜ −A w wm
( A˜ − A ) W ˜ = max
1i ,l N
l,m
j
l,m
j ,m=1
max
1i , j ,l,mn
N (i , j ) (i , j ) (N ) (N ) A˜ − A w j wm l,m l,m j ,m=1
=
max
1i , j ,l,mn
N (i , j ) (i , j ) (N ) A˜ − A wj l,m l,m j =1
= O n−k log4 n . (ii) From
2
α / N W (j N ) β/ N, it follows [ w (jN ) ]2 β 2 / N 2 , j = 1, . . . , N. Therefore
( A˜ − A ) W ˜ 2 = F
N (i , j ) (i , j ) ( N ) ( N ) 2 A˜ −A w wm l,m
l,m
j
i , j ,l,m=1
max
1i , j ,l,mn
max
1i , j ,l,mn
=O n
−2k
(i , j ) (i , j ) 2 A˜ −A l,m
N
l,m
(N )
wj
( N ) 2
wm
i , j ,l,m=1 N (i , j ) 2 2 2 (i , j ) 2 2 A˜ − A N β /N l,m l,m j ,l=1
log n . 8
˜ − A)W ˜ F = O(n−k log4 n). That is, ( A
2
To end this subsection, we make the following remarks. (N )
1. For many quadrature rules, the conditions about W j ( j = 1, 2, . . . , N) are satisfied (for example, repeated low-order Newton–Cotes rules). ˜ is independent of the number of quadrature points N, that is, na can be chosen 2. We observe that the accurate in A a W ˜ )f = g will be very without considering the size of the linear system. If na is large enough, the solution of ( I − A a W ˜ )f = g. close to that of ( I − A W ˜ = O(n−k log4 nb ) will be small, and therefore the convergence of the RC 3. By choosing a moderate nb , ( A a − B a ) W b scheme (3.1) will be very fast. 3.2. Implementation and complexity analysis
˜ )−1 , which In this subsection, we first discuss how to compute the matrices P A and Λ A , P B and Λ B , and ( I − B a W are required in the iterative method (3.1) and then we discuss the implementation of the RC scheme. Details are listed as follows. The computation of Λ A and Λ B . The matrix Λ A can be obtained in two steps: 1) evaluate the function values of a(x, y , u , v ) at the Chebyshev grid to form the na2 × na2 matrix H g ; 2) apply two dimensional fast cosine transform on H g . The cost is O(na4 + na4 log na ) = O(na4 log na ). Similarly, the cost of obtaining Λ B is O(nb4 log nb ). The computation of P A and P B . Since that P A and P B are of size na × N and nb × N, the cost of obtaining them is O((na + nb ) N ). The computation of the residual vector. From (3.2) and (3.6), we get
r(q+1) = r(q) + I − ( P A ⊗ P A ) T Λ A ( P A ⊗ P A )( W ⊗ W ) d(q) . The computation of e = ( W ⊗ W )d(q) requires N 2 ops, the computation of ( P A ⊗ P A )e can be done in 2na N 2 ops by using ( P A ⊗ P A )e = vec( P A mat(e) P TA ), where mat(e) is the N × N matrix obtained by reordering the vector e (the first N entries form the first column of mat(e) and so on) and vec( A ) is the vector obtained by column-by-column ordering of A. Thus the cost of computing the residual vector is about N 2 + 4N 2 na + na4 ops (multiplications).
W.-J. Xie, F.-R. Lin / Applied Numerical Mathematics 59 (2009) 1709–1719
1717
A fast solution method for the preconditioning system. In the RC scheme (3.1), we require to solve the preconditioning system (3.3). We show in the following that by using the decomposition of B a , we can solve the equation efficiently. By using the Sherman–Morrison formula (see for instance [5, p. 74]), we see that
( I M − U V )−1 = I M + U ( I m − V U )−1 V , where I M and I m are identity matrices, and U and V are M × m and m × M matrices respectively. Note that
−1
I N 2 − B a (W ⊗ W )
−1 = I N 2 − ( P B ⊗ P B )T Λ B ( P B ⊗ P B )( W ⊗ W ) .
Let U = ( P B ⊗ P B ) T and V = Λ B ( P B ⊗ P B )( W ⊗ W ) = Λ B (( P B W ) ⊗ ( P B W )), we have
−1
I N 2 − B a (W ⊗ W )
−1 = I N 2 + U I n2 − Λ B P B W P BT ⊗ P B W P BT V. b
(3.7)
The matrix ( I n2 − Λ B (( P B W P BT ) ⊗ ( P B W P BT )))−1 can be obtained efficiently in O( N (nb + nb2 ) + nb6 ) ops by first computing b
P B W P BT (about N (nb + nb2 ) ops). The preconditioning system can be solved by using (3.7) in N 2 + 2N 2 nb + 2Nnb2 + 2nb4 ops.
The computational costs. In summary, the cost for constructing the required matrices is O(na4 log na + N (na + nb2 ) + nb6 ), the cost of the computation of the residual vector is O( N 2 na + na4 ), and the cost for solving the preconditioning system O( N 2 nb + Nnb2 + nb4 ). Storage requirement. We require to store the six matrices: P A , P B , Λ A , Λ B , W , and ( I − Λ B (( P B W P BT ) ⊗ ( P B W P BT )))−1 , which contain na N, nb N, na4 , nb4 , N and nb4 entries respectively. We also require to store three N 2 × 1 vectors r(q) , d(q) , f(q) . Therefore, the total Storage requirement is O( N 2 + Nna + na4 ). 4. Numerical examples In this section, we use our algorithms to solve one dimensional and two dimensional Fredholm equations of the second kind. We use the composite 4-point Gauss–Legendre quadrature to discretize the integral equations. Let the roots of the degree-4 Legendre polynomial be denoted by ξi (i = 1, 2, 3, 4), then the 4-point Gauss–Legendre quadrature is given by
1 f (x) dx =
4
ηi f (ξi ) + E 4 ( f ),
i =1
−1
where ηi (i = 1, 2, 3, 4) are the weights and E 4 ( f ) is the error in the quadrature. For the one dimensional equation
1 f (x) −
a(x, y ) f ( y ) dy = g (x),
−1 x 1,
−1
we partition the interval [−1, 1] into n subintervals of length d = 2/n. In this case, the weight matrix W is given by W =
d 2
I n ⊗ diag(η1 , η2 , η3 , η4 )
and the quadrature points are
(N )
ti
= −1 +
i−1 4
+
1 2
d+
(4.1)
d
ξ i −1 , 2 i −4[ 4 ]
i = 1, 2, . . . , N = 4n,
(4.2)
where [x] denotes the maximal integer that is not larger than x. For the two dimensional equation
1 1 f (x, y ) −
a(x, y , u , v ) f (u , v ) du dv = g (x, y ),
−1 x, y 1,
−1 −1
˜ = W ⊗ W with W given by (4.1) and the quadrature points are (t (N ) , t (N ) ), i , j = 1, 2, . . . , N, where the weight matrix is W i j (N )
for i = 1, 2, . . . , N are given by (4.2). The RC scheme is used to solve the discretization systems of one dimensional and two dimensional equation respectively. The initial guesses are set to zero vectors and the stopping criterion is
ti
(q) r = g − ( I − A a W )f(q) < 10−6 (for one dimensional equation)
or
(q) r = g − ( I − A a W ˜ )f(q) < 10−6 (for two dimensional equation),
1718
W.-J. Xie, F.-R. Lin / Applied Numerical Mathematics 59 (2009) 1709–1719
where f(q) is the approximate solution at the qth iteration. The error of the numerical solution f approx is computed by e ( f approx ) = f approx − f true ∞ . Two examples are tested. The first example is a one dimensional Fredholm equation:
1 f (x) −
a(τ ; x, t ) f (t ) dt = 1 − arctan
τ (1 − x) 2
−1
where
τ > 0 is a parameter and a(τ ; x, t ) = ∂ j a(τ ; x, t ) = ∂xj
j +1
τ
2
∂ j a(2; x, t ) ∂xj
− arctan
τ /2 1+(τ (x−t )/2)2
2
x ∈ [−1, 1],
,
(4.3)
. The true solution of (4.3) is f true (x) = 1 for x ∈ [−1, 1]. Since
∂ j a(τ ; x, t ) = ∂t j
and
τ (1 + x)
j +1
τ
2
∂ j a(2; x, t ) , ∂t j
a large τ will lead to large partial differentials of a(τ ; x, t ). Therefore, for a given increase in order that ( A a − A ) W < . The second example is a two dimensional Fredholm equation:
> 0, the value of na increases as τ
1 1 f (x, y ) −
a(τ ; x, y , u , v ) f (u , v ) du dv = g (τ ; x, y ),
(x, y ) ∈ [−1, 1] × [−1, 1],
(4.4)
−1 −1
where
τ > 0 is a parameter, a(τ ; x, y , u , v ) =
τ2 4
exp −
τ (1 + x) τ (1 + u ) 2
2
−
τ (1 + y ) τ (1 + v ) 2
2
and g (τ ; x, y ) = 1 −
2 2 2 e −τ (1+x+ y )/2 − e −τ (1+ y )/2 − e −τ (1+x)/2 + 1 .
4
τ 2 (1 + x)(1 + y )
The true solution of (4.4) is f true (x, y ) = 1. Similar to the one dimensional case, we have
∂ j a(τ ; t 1 , t 2 , t 3 , t 4 ) j ∂ ti
=
j +2
τ
∂ j a(2; t 1 , t 2 , t 3 , t 4 )
2
∂ ti
j
,
i = 1, 2, 3, 4,
and it follows that a large τ leads to a large na for a fixed accuracy . In Tables 1 and 2, we show the numbers of iterations which are denoted by #I and the errors in the numerical solutions which are denoted by Error for Eqs. (4.3) and (4.4) with different τ and N, respectively. We observe from the numerical results that: 1. The RC scheme converges very fast. 2. Given τ , for fixed na and nb , the numbers of iterations are the same for all N. This result is consistent with the ˜ depend only on na and nb . ˜ )−1 ( A a − B a ) W theoretical result that ( I − B a W )−1 ( A a − B a ) W and ( I − B a W 3. Given τ , for fixed na and nb , the error in the numerical solution will decrease as N increases when N is small, which indicates that the error in the quadrature rule is larger than that in the interpolating polynomial. When N is large enough, increasing the value of N will not improve the accuracy of the numerical solution, which indicates that the error in the interpolating polynomial is dominating in this case. Finally, we make a briefly remark on the relation of na , nb and N. In order that the error in the numerical solution decreases at N increase for large N, we should increase the value of na while fix the value of nb . Suppose the error in the numerical solution of (1.3), the discretization linear system of the original integral equation, is of order O( N −l ). Table 1 Numbers of iterations for the RC scheme and errors in the numerical solutions for one dimensional Fredholm equation (4.3) for the cases na = 32, nb = 8 (left), τ = 16, na = 128, nb = 32 (middle), and τ = 64, na = 320, nb = 120 (right). N
256 512 1024 2048 4096
(τ , na , nb ) = (4, 32, 8)
τ = 4,
(τ , na , nb ) = (16, 128, 32)
(τ , na , nb ) = (64, 320, 120)
#I
Error
#I
Error
#I
Error
6 6 6 6 6
1.80e−09 1.80e−09 3.23e−10 3.24e−10 3.23e−10
6 6 6 6 6
1.59e−09 1.08e−09 2.33e−10 6.29e−11 2.85e−10
6 6 6 6 6
1.56e−5 3.16e−7 9.12e−8 2.30e−8 1.68e−8
W.-J. Xie, F.-R. Lin / Applied Numerical Mathematics 59 (2009) 1709–1719
1719
Table 2 Numbers of iterations for the RC scheme and errors in the numerical solutions for two dimensional Fredholm equation (4.4) for the cases τ = 4, na = 16, nb = 6 (left) and τ = 8, na = 32, nb = 12 (right).
(τ , na , nb ) = (4, 16, 6)
N
32 64 128 256 512
(τ , na , nb ) = (8, 32, 12)
#I
Error
#I
Error
7 7 7 7 7
1.30e−06 1.01e−06 1.31e−06 1.30e−06 1.31e−06
8 8 8 8 8
1.63e−03 2.07e−05 6.74e−07 4.87e−07 6.74e−07
From Theorem 2 we see that the error in the approximate matrix is O(na−k log4 na ). Thus, in order that the error in the solution of (1.3) and that in the approximate matrix are of the same order, we require O( N −l ) = O(na−k log4 na ). Since that log4 na = o(na ) for large na , it follows that if we choose na such that
na = O N l/(k−1) ,
(4.5)
then the error in solution of (1.6) is also O( N −l ). Recall that the value of l depends on the quadrature rule used, and the smoothness of the solution and the kernel function. If the solution of (1.1) has some kind of singularity, then l and l/(k − 1) are small and it follows that na will increase slowly as N increases. On the other hand, if the solution of (1.1) is very smooth (as our examples) and a high order quadrature is used, then na will increase quickly as N increases. 5. Concluding remarks Based on interpolating polynomial approximation, we proposed new algorithms for the fast numerical solution of two dimensional Fredholm equation of the second kind with smooth kernel. According to (4.5) and the discussions in Section 3.2, we see that the complexity of our algorithm is: 1. Storage requirement: O( N 2 + Nna + na4 ) = O( N 2 + N (k+l−1)/(k−1) + N 4l/(k−1) ). 2. The cost per iteration of the RC scheme: 1) the cost of the computational of residual vector is N 2 + 4N 2 na + na4 = O( N (2k+l−2)/(k−1) + N 4l/(k−1) ); 2) the cost of the solution of the preconditioning system is N 2 + 2N 2 nb + 2Nnb2 + 2nb4 = O( N 2 ). We observe that the main work of the iteration is the computation of the residual vector. Numerical results were given to illustrate the efficiency of our methods. Acknowledgement We thank the anonymous referees for their detailed and valuable comments and suggestions. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21]
B. Alpert, G. Beylkin, R. Coifman, V. Rokhlin, Wavelets for the fast solution of second-kind integral equations, SIAM J. Sci. Comp. 14 (1993) 159–184. P.M. Anselone, Collectively Compact Operator Approximation Theory, Prentice-Hall, Englewood Cliffs, NJ, 1971. K.E. Atkinson, A Survey of Numerical Methods for the Solution of Fredholm Integral Equations of the Second Kind, SIAM, Philadelphia, PA, 1976. K.E. Atkinson, The Numerical Solution of Integral Equations of the Second Kind, Cambridge University Press, Cambridge, 1997. O. Axelsson, Iterative Solution Methods, Cambridge University Press, Cambridge, 1996. C.T. Baker, The Numerical Treatment of Integral Equations, Clarendon Press, Oxford, 1977. Z. Battles, L.N. Trefethen, An extension of MATLAB to continuous functions and operators, SIAM J. Sci. Comp. 25 (2004) 1743–1770. G. Beylkin, R. Coifman, V. Rokhlin, Fast wavelet transforms and numerical algorithms I, Comm. Pure Appl. Math. 46 (1991) 141–183. L. Brutman, Lebesgue functions for polynomial interpolation – a survey, Ann. Numer. Math. 4 (1977) 111–127. L.M. Delves, J.L. Mohamed, Computational Methods for Integral Equations, Cambridge University Press, New York, 1985. R. Farengo, Y.C. Lee, P.N. Guzdar, An electromagnetic integral equation: application to microtearing modes, Phys. Fluids 26 (1983) 3515–3523. G. Golub, C. Van Loan, Matrix Computations, second ed., Johns Hopkins University Press, Baltimore, 1989. L. Greengard, V. Rokhlin, A fast algorithm for particle simulations, J. Comput. Phys. 73 (1987) 325–348. L. Greengard, V. Rokhlin, A new version of the fast multipole method for the Laplace equation in three dimensions, Acta Numerica 6 (1997) 229–269. G.Q. Han, L.Q. Zhang, Asymptotic error expansion and extrapolation of the Nystrom eethord for two-dimensional Fredholm integral equations, Acta Math. Appl. Sinica 18 (2) (1995) 218–224 (in Chinese). L.C. Hsu, Y.S. Zhou, Y.B. Sun, Approximation Theory, National Defencial Industry Publishing Company, China, 1985 (in Chinese). V.V. Ivanov, V.K. Zadiraka, Some new results in theory of approximation, in: Computational Mathematics, Kiev University Press, Kiev, 1966, pp. 30–36 (in Russian). F.R. Lin, Preconditioned iterative methods for the numerical solution of Fredholm equations of the second kind, Calcolo 40 (2003) 231–248. J. Tausch, The variable order fast multipole method for boundary integral equations of the second kind, Computing 72 (2004) 267–291. Y. Wang, Y. Xu, A fast wavelet collocation method for integral equations on polygons, J. Integral Equations Appl. 17 (2005) 277–330. A.S. Xia, L. Wang, Iterative corrections for finite element solutions of two dimensional second kind Fredholm integral equations, J. Inst. Command Technology 11 (2) (2000) 105–108 (in Chinese).
Contents lists available at ScienceDirect
Applied Numerical Mathematics www.elsevier.com/locate/apnum
A fast numerical solution method for two dimensional Fredholm integral equations of the second kind ✩ Wen-Jing Xie, Fu-Rong Lin ∗ Department of Mathematics, Shantou University, Shantou Guangdong 515063, China
a r t i c l e
i n f o
Article history: Received 23 August 2006 Received in revised form 3 July 2008 Accepted 30 January 2009 Available online 7 February 2009 MSC: 45L10 65R20 Keywords: Integral equation Polynomial interpolation Approximate matrix Residual correction scheme
a b s t r a c t In this paper we consider numerical solution methods for two dimensional Fredholm integral equation of the second kind
1 1 f ( x, y ) −
a(x, y , u , v ) f (u , v ) du dv = g (x, y ),
(x, y ) ∈ [−1, 1] × [−1, 1],
−1 −1
where a(x, y , u , v ) is smooth and g (x, y ) is in L 2 [−1, 1]2 . We discuss polynomial interpolation methods for four-variable functions and then use the interpolating polynomial to approximate the kernel function a(x, y , u , v ). Based on the approximation we deduce fast matrix-vector multiplication algorithms and efficient preconditioners for the above two dimensional integral equations. The residual correction scheme is used to solve the discretization linear system. © 2009 IMACS. Published by Elsevier B.V. All rights reserved.
1. Introduction Many problems in engineering and mechanics can be transformed into two dimensional Fredholm integral equations of the second kind. For example, it is usually required to solve Fredholm integral equations in the calculation of plasma physics [11]. There are many works on developing and analyzing numerical methods for solving Fredholm integral equations of the second kind (see [1–4,6,10]). In recent years, a number of algorithms for the fast numerical solution of integral equations of the second kind have been developed; see, e.g., [1,8,13,14,18–20]. The fast multipole method proposed in [13] combines the use of low-degree polynomial interpolation of the kernel functions with a divide-and-conquer strategy. For kernel functions that are Coulombic or gravitational in nature, it results in an order O( N ) algorithm for the matrix-vector multiplications, where N is the number of quadrature points used in the discretization. In [8], an O( N log N ) algorithm was developed by exploiting the connections between the use of wavelets and their applications on Calderon–Zygmund operators. In [18], we proposed an approximation scheme to obtain low rank approximation of the discretization matrix. By using preconditioned iterative methods such as the residual correction (RC) scheme or the preconditioned conjugate gradient (PCG) method, the total cost for solving the integral equation with smooth kernel is O( N ) operations (ops). The papers [14,19,20] are concerned with fast algorithms for two dimensional integral equations with weakly singular kernel functions. In [14] and [19], different versions of fast multipole methods are proposed to solve integral equations with weakly singular kernels. In [20], the authors developed a fast wavelet collocation method for the integral equations defined ✩
*
The research was supported by NSF of China No. 10271070. Corresponding author. E-mail address: [email protected] (F.-R. Lin).
0168-9274/$30.00 © 2009 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.apnum.2009.01.009
1710
W.-J. Xie, F.-R. Lin / Applied Numerical Mathematics 59 (2009) 1709–1719
on polygons. Besides fast algorithms, there are many investigations on improving the accuracy of numerical solutions, see for instance [15] and [21]. In [15], a Nyström method for two dimensional Fredholm integral equation of second kind was used. The asymptotic expansion of numerical solution is obtained and then the Richardson extrapolation method is used to enhance the precision of the numerical solution. In [21], based on the finite-element solution, an iterative correction method was proposed. It shows that if the kernel function is smooth, the precision could be improved considerably. In this paper, we consider the numerical solution of two dimensional Fredholm integral equation of the second kind with smooth kernel:
1 1 f (x, y ) −
a(x, y , u , v ) f (u , v ) du dv = g (x, y ),
(x, y ) ∈ [−1, 1] × [−1, 1],
(1.1)
−1 −1
where the kernel function a(x, y , u , v ) is smooth in [−1, 1]4 , the unknown function f (x, y ) and the right-hand side function g (x, y ) are in L 2 [−1, 1]2 . For the case where the integration domain is not [−1, 1]2 , say [α , β]2 , we can use the linear transformation
x = (2x − α − β)/(β − α ), y = (2 y − α − β)/(β − α )
to transform the domain into [−1, 1]2 . (N ) (N ) (N ) We discretize the equation by using certain numerical quadrature. Let −1 t 1 < t 2 < · · · < t N 1 be the quadrature (N )
points and w i
, i = 1, 2, . . . , N, be the corresponding weights, a numerical quadrature is defined as
1 f (t ) dt ≈
N
(N )
wi
(N )
f ti
(1.2)
.
i =1
−1
Typical examples include Gaussian rules, repeated Gaussian rules, and repeated Newton–Cotes rules, etc., see for instance [10]. Extending the above formula to two dimensional case, we get the following approximate system for (1.1): f (x, y ) −
N N
(N )
wi
(N )
(N )
w j a x, y , t i
, t (jN ) f t i(N ) , t (jN ) = g (x, y ),
x, y ∈ [−1, 1].
i =1 j =1
In particular, we get the discretization linear system of (1.1):
˜ )f = g, (I − A W
(1.3) (N )
where f and g are vectors obtained by reordering the matrices [ f (t i
(N )
(N )
, t j )]iN, j =1 and [ g (t i
(N )
, t j )]iN, j=1 row-by-row respec-
˜ = W ⊗ W with W the weight matrix diag( w , w , . . . , w ), and tively, I is the identity matrix, W N 1 2 (N )
⎛
A (1,1)
⎜ (2,1) ⎜ A ⎜ A=⎜ . ⎜ .. ⎝ A ( N , 1) (N )
with A (i , j ) = [a(t i
···
A (1, N )
A (2,2)
A (2, N ) ⎟
.. .
··· .. .
A ( N , 2)
···
A (N ,N )
(N )
, tl
(N )
(N )
⎞
A (1,2)
.. .
(N )
⎟ ⎟ ⎟ ⎟ ⎠
(1.4)
(N )
, t j , tm )]lN,m=1 . Here ⊗ denotes the Kronecker tensor product.
To solve the linear system (1.3) by classical direct methods such as Gaussian elimination method requires O( N 6 ) operations (ops). For iterative methods such as the conjugate gradient (CG) method (see [12]), each iteration requires O( N 4 ) ops. Therefore even for well-conditioned problems, the CG method requires O( N 4 ) ops, which for large-scale problems is often prohibitive. In this paper, we propose a fast numerical solution method for the linear system (1.3) by using the polynomial interpolation technique. That is, rather than solving the discretization linear system (1.3), we solve a Nyström approximation of the following equation
f
[n]
1 1 (x, y ) −
an (x, y , u , v ) f [n] (u , v ) du dv = g (x, y ),
−1 x, y 1,
(1.5)
−1 −1
where an (x, y , u , v ) is an interpolation based approximation to a(x, y , u , v ) in which n Chebyshev nodes are used in each space variable. More precisely, we solve the following approximation equation of (1.3):
˜ )f[n] = g, ( I − Aa W
(1.6)
W.-J. Xie, F.-R. Lin / Applied Numerical Mathematics 59 (2009) 1709–1719
(N )
(N )
(N )
1711
(N )
where A a = ([ A a ](i , j ) )iN, j =1 with [ A a ](i , j ) = [an (t i , tl , t j , tm )]lN,m=1 . The outline of this paper is as follows. In Section 2, we discuss how to interpolate multi-variable functions efficiently and analyze the error in the interpolating polynomials. In Section 3, we propose fast algorithms for solving (1.6) by using the results obtained in Section 2 and discuss the complexity of our algorithms. Numerical results are given in Section 4 to illustrate the efficiency and the accuracy of our algorithms. Finally, concluding remarks are given in Section 5. 2. Interpolation of multi-variable functions In this section we discuss the approximation of two-variable and four-variable smooth functions h(x, y ) (−1 x, y 1) and a(x, y , u , v ) (−1 x, y , u , v 1) by polynomial interpolation at Chebyshev grids. We discuss how to obtain the interpolating polynomials efficiently and analyze the error in the approximation. 2.1. Interpolation (n)
Let c j (x), j = 0, 1, . . . , n − 1, be the scaled Chebyshev polynomials defined by
(n)
c j (x) = δi cos j arccos(x) ,
√
Here δ0 = 1/ n and δ j =
j = 1, 2, . . . , n − 1, x ∈ [−1, 1].
√
2δ0 for j = 1, . . . , n − 1. Let (2i − 1)π (n) p i = cos , i = 1, 2, . . . , n,
2n
the roots of the degree-n Chebyshev polynomial cos(n arccos(x)), we call
(n)
(n)
pi , p j
i , j = 1, 2, . . . , n
the n × n Chebyshev grid of the domain [−1, 1]2 . Let the interpolating polynomial of h(x, y ) at the Chebyshev grid of [−1, 1]2 be denoted by (C n h)(x, y ), we have
(Cn h)(x, y ) =
n n p =1 q=1
(n)
(n)
λ pq c p −1 (x)cq−1 ( y ),
−1 x, y 1,
(2.1)
where the coefficients λ pq are determined by the interpolating conditions:
(n)
(n)
h pi , p j
=
n n p =1 q=1
λ pq c (pn−) 1 p (i n) cq(n−)1 p (jn) ,
i , j = 1, 2, . . . , n.
(2.2)
(n) (n) Let Λ = [λi j ]ni, j =1 and c(n) (·) = (c 0 (·), . . . , cn−1 (·)) T , the formula (2.1) can be written as
(Cn h)(x, y ) = c(n) (x)T Λc(n) ( y ). (n)
(n)
Let H g = [h( p i , p j )]ni, j =1 , and
(n) n i , j =1
C = c i −1 p j
n (i − 1)(2 j − 1)π = δi −1 cos 2n
,
i , j =1
the conditions (2.2) can be written in compact form as H g = C T ΛC . Note that C is the cosine transform matrix, which is orthogonal, we have Λ = C H g C T and it follows that Λ can be obtained in O(n2 log n) ops by using fast cosine transform. In summary, we get an approximation of h(x, y ): h(x, y ) ≈ (Cn h)(x, y ) =
n n
C HgCT
p =1 q=1
(n) (n) c (x)cq−1 ( y ). p ,q p −1 (n)
(n)
(n)
(n)
Similarly, the interpolating polynomial of h(x, y , u , v ) at the Chebyshev grid of ({( p i , p j , pl , pm ) | i , j , l, m = 1, . . . , n}) is given by
(Cn h)(x, y , u , v ) =
n n n n p =1 q=1 r =1 s=1
λ pqrs c (pn−) 1 (x)cq(n−)1 ( y )cr(n−)1 (u )c (sn−)1 ( v )
(2.3)
for −1 x, y , u , v 1, where the coefficients λ pqrs are determined by the interpolating conditions:
(n)
(n)
(n)
(n)
h p i , p j , pl , pm
=
n n n n p =1 q=1 r =1 s=1
(n) (n) (n) (n) (n) (n) (n) (n) λ pqrs c p −1 p i cq−1 p j cr −1 pl c s−1 pm ,
(2.4)
1712
W.-J. Xie, F.-R. Lin / Applied Numerical Mathematics 59 (2009) 1709–1719
for i , j , l, m = 1, 2, . . . , n. Similar to the two-variable case, we have
T (Cn h)(x, y , u , v ) = c(n) (x) ⊗ c(n) ( y ) Λ c(n) (u ) ⊗ c(n) ( v ) ,
where
⎛
Λ(1,1) Λ(1,2) · · · Λ(1,n)
(2.5)
⎞
⎟ ⎜ (2,1) ⎜Λ Λ(2,2) · · · Λ(2,n) ⎟ ⎟ ⎜ Λ = ⎜ .. .. .. .. ⎟ ⎟ ⎜ . . . . ⎠ ⎝ ( n , 1) ( n , 2) ( n,n) Λ Λ ··· Λ
with Λ(i , j ) = [λiljm ]ln,m=1 .
(i , j ) n (i , j ) ]i , j=1 with n × n blocks H g
Let H g be the n × n block matrix H g = [ H g get
(n)
(n)
(n)
(n)
= [h( p i , pl , p j , pm )]ln,m=1 . From (2.4) we
H g = (C ⊗ C ) T Λ(C ⊗ C ). Therefore,
Λ = (C ⊗ C ) H g (C ⊗ C ) T ,
(2.6) 4
which can be obtained by using fast cosine transform in O(n log n) ops. 2.2. Error analysis In this subsection, we discuss the error in the interpolating polynomial of a(x, y , u , v ) at the Chebyshev grid ((Cn a)(x, y , u , v )). We first discuss the error in the interpolating polynomials for one-variable functions. Lemma 1 (Jackson’s Fifth Theorem, see for instance [16, p. 95]). Suppose f ∈ C [−1, 1] and f (k) ∈ C [−1, 1] for some positive integer k. Let (Bn f )(x) be the best approximation polynomial of degree-(n − 1) for f (x) ∈ C [−1, 1], i.e.,
Bn f − f =
min
qn :degree-(n−1) polynomial
qn − f .
Here and throughout the paper, f = maxx∈[−1,1] | f (x)| is the supremum norm. Then for n k, we have
E n ( f ) ≡ Bn f − f
π
k
2
(k) f / (n + 1)n · · · (n − k + 2)
π
k
2
(k) f /(n − k + 2)k .
The Lebesgue function and Lebesgue constant are useful tools for analyzing the error of interpolating polynomials. For a given interpolating grid {x1 , . . . , xn }, the Lebesgue function is defined as L n (x) =
n hk (x), k=1
where hk (x) =
n
x−xi i =1,i =k xk −xi
is the kth Lagrange basis polynomial, and the Lebesgue constant is defined as
Ωn = Ln , see for instance [9]. For Chebyshev grids, we have (see [17])
Ωn < 2 +
2
π
From (Cn f )(x) =
log n.
n
k=1
(2.7) (n)
f ( pk )hk (x), it follows that
n (n) (Cn f )(x) f p hk (x) Ωn f .
(2.8)
k
k=1
Lemma 2. Suppose f (x) ∈ C k [−1, 1]. If n > k, then
f (x) − (Cn f )(x) δn,k f (k) = O n−k log n ,
where
δn,k =
π 2
k 3+
2
π
x ∈ [−1, 1],
log n (n − k + 2)−k .
(2.9)
W.-J. Xie, F.-R. Lin / Applied Numerical Mathematics 59 (2009) 1709–1719
1713
Proof. It is easily seen that Cn f is a linear operator and Cn Bn f = Bn f . Thus
(Cn f )(x) − f (x) f (x) − (Bn f )(x) + (Cn Bn f )(x) − (Cn f )(x) E n ( f ) + Cn (Bn f − f ) (x) E n ( f ) + Ωn Bn f − f by (2.8) = (1 + Ωn ) E n ( f ). Using (2.7) and the result of Lemma 1, the result of the lemma follows from the above inequality.
2
From the above lemma, we see that Chebyshev grids are very powerful in high degree polynomial approximation. In [7], continuous functions defined in [−1, 1] are approximated very accurately by using the polynomial interpolant with sufficiently many Chebyshev points. We now extend the above result to interpolating polynomials of two-variable functions. Lemma 3. Suppose the function a(u , v ) is such that ∂ k a/∂ uk and ∂ k a/∂ v k are continuous in the domain [−1, 1]2 , where k is a positive integer. Let Cn a be the degree-(n − 1) interpolating polynomial of a(u , v ) at the Chebyshev grid (cf. (2.1)). Then
a − Cn a
π
k
2
3+
2
π
2 log n
(n − k + 2)−k c = O n−k log2 n
for n > k, where c = max{∂ k a/∂ uk , ∂ k a/∂ v k }. Proof. First of all, we define the one-variable interpolating operators for two-variable function: L v a (u , v ) =
n
(n) (n)
a u, pi
hi ( v )
(2.10)
i =1
and L u a (u , v ) =
n
(n)
(n)
a p i , v h i (u ),
i =1
(n)
where h i
(n)
(n)
(n)
is the ith Lagrange basis polynomial at the Chebyshev grid { p 1 , p 2 , . . . , pn }. Obviously, the interpolating (n)
(n)
polynomial (Cn a)(u , v ) of a(u , v ) at the Chebyshev grid {( p i , p j ) | i , j = 1, 2, . . . , n} is given by
(Cn a)(u , v ) =
n n
(n)
(n) (n)
a pi , p j
(n)
h i (u )h j (u ) = ( L u L v a)(u , v ) = ( L v L u a)(u , v ).
i =1 j =1
By the triangular inequality a − Cn a a − L v a + L v a − L u L v a, we only require to estimate a − L v a and L v a − L u L v a respectively. By using Lemma 2 and (2.10), we have
k k a(u , v ) − ( L v a)(u , v ) δn,k max ∂ a (u , t ) δn,k ∂ a , ∂ vk −1t 1 ∂ v k
Thus
−1 u , v 1.
k ∂ a a − L v a δn,k ∂ v k .
(2.11)
As for L v a − L u L v a, by using (2.8), we have
k ∂ a L v a − L u L v a = L v (a − L u a) Ωn a − L u a Ωn δn,k ∂ uk .
Let c = max{∂ k a/∂ uk , ∂ k a/∂ v k }, it follows from (2.11) and (2.12) that
a − Cn a δn,k (1 + Ωn )c . Therefore, from (2.7) and (2.9), the result follows.
2
Now we can prove our main theorem on the error of interpolating polynomials for four-variable functions.
(2.12)
1714
W.-J. Xie, F.-R. Lin / Applied Numerical Mathematics 59 (2009) 1709–1719
Theorem 1. Let a(x, y , u , v ) be a function defined in the domain [−1, 1]4 such that ∂ k a/∂ xk , ∂ k a/∂ yk , ∂ k a/∂ uk and ∂ k a/∂ v k are continuous in the domain [−1, 1]4 , where k is a positive integer. Let (Cn a)(x, y , u , v ) be the interpolating polynomial of a(x, y , u , v ) at the Chebyshev grid (cf. (2.3)). Then for n > k, we have
a − Cn a = O n−k log4 n .
(2.13)
Proof. Let L x , L y , L u , and L v be the one-variable interpolating operators which are defined as in the proof of Lemma 3, we have Cn a = L x L y L u L v a. Therefore
a − Cn a = a − L x L y L u L v a a − L v a + L v a − L u L v a + L u L v a − L y L u L v a + L y L u L v a − L x L y L u L v a. By using Lemma 2, we have a − L v a δn,k ∂ k a/∂ v k . Moreover, by using (2.8), we get
L v a − L u L v a L v a − L u a Ωn δn,k ∂ k a/∂ uk , L y L u L v a − L u L v a L u L v ( L y a − a) Ωn2 δn,k ∂ k a/∂ yk , L y L u L v a − L x L y L u L v a = L y L u L v (a − L x a) Ωn3 δn,k ∂ k a/∂ xk , where Ωn is the Lebesgue constant. It follows that
a − Cn a δn,k 1 + Ωn + Ωn2 + Ωn3 max ∂ k a/∂ xk , ∂ k a/∂ yk , ∂ k a/∂ uk , ∂ k a/∂ v k . Using (2.7) and (2.9), we obtain (2.13).
2
3. Fast algorithms for integral equations In this section, we propose new algorithms for the fast numerical solution of two dimensional Fredholm equations of the second kind by using the interpolation method discussed in Subsection 2.1. We approximate the matrix A defined by (1.4) with A a such that
( I − A W ˜ <
˜ ) − ( I − Aa W ˜ ) = ( A − A a ) W for given > 0, so that the solution of the approximation linear system (1.6) approximates the solution of (1.3) accurately. ˜ ) as We then approximate A is by a less accurate matrix B a and split ( I − A a W
˜ = (I − Ba W ˜ ) − ( Aa − B a )W ˜. I − Aa W Thus, Eq. (1.6) is equivalent to
˜ f[n] + g. ˜ )f[n] = ( A a − B a ) W (I − Ba W In this way, we come to the following iterative method:
˜ )f(q+1) = ( A a − B a ) W ˜ f(q) + g, (I − Ba W
q = 0, 1, 2, . . . ,
(3.1)
or equivalently
˜ )f(q) , r(q) = g − ( I − A a W
(3.2)
˜ )d(q) = r(q) , (I − Ba W
(3.3)
(q+1)
f
(q)
=f
(q)
+d .
Eq. (3.3) is called the preconditioning system. The above iterative method is often called a residual correction (RC) scheme or a method of iterative improvement. We note that the main cost per iteration of the RC scheme is the matrix-vector ˜ )f(q) and the solution of the preconditioning system (3.3), and the method converges if and only if multiplication ( I − A a W
ρ = ( I − B a W˜ )−1 ( A a − B a ) W˜ < 1. Moreover, the smaller the
ρ , the faster the convergence.
W.-J. Xie, F.-R. Lin / Applied Numerical Mathematics 59 (2009) 1709–1719
1715
3.1. The approximate matrices By (2.5), we have
(N )
a ti
(N ) (N ) , t (jN ) , tl(N ) , tm ≈ (Cn a) t i(N ) , t (jN ) , tl(N ) , tm (N ) (N ) T (n) (N ) (N ) = c(n) t i Λ c tl ⊗ c(n) t j ⊗ c(n) tm
for i , j , l, m = 1, . . . , N. It is follows that
˜ = ( P ⊗ P )T Λ( P ⊗ P ) A≈A where
(n)
(3.4)
( N ) j =1: N . i =1:n
P = c i −1 t j
(3.5)
Based on the discussion in Subsection 2.1, we have the following algorithm for computing Λ and P . Algorithm 1 (The construction of approximation matrix). (n)
Step 1. Choose n, form the Chebyshev points: p i , i = 1, 2, . . . , n and compute the function values of a(x, y , u , v ) at the Chebyshev grid of the domain [−1, 1]4 :
⎛
(1,1)
Hg
⎜ (2,1) ⎜ Hg Hg = ⎜ ⎜ .. ⎝ .
(1,2)
· · · H (g1,n)
⎞
(2,2)
··· Hg .. .. . .
⎟ ⎟ ⎟ ⎠
Hg Hg
.. .
(n,1)
Hg
(n)
(n)
Hg (i , j )
(n,2)
(2,n) ⎟
(n,n)
··· Hg
(n)
(n)
with H g = [a( p i , pl , p j , pm )]ln,m=1 . Step 2. Compute Λ by using two dimensional fast cosine transform on H g (cf. (2.6)). Step 3. Compute the matrix P which is defined in (3.5). We note that both A a and B a can be obtained by using Algorithm 1. We choose a large n (denoted by na ) for computing A a and a moderate n (denoted by nb ) for computing B a , see numerical results in Section 4. In the following, we will denote A a and B a as A a = ( P A ⊗ P A ) T Λ A ( P A ⊗ P A ) and
B a = ( P B ⊗ P B )T Λ B ( P B ⊗ P B )
(3.6)
respectively. It is easily seen that P A , P B , Λ A , and Λ B are of the size na × N, nb × N, na2 × na2 , and
˜ we have the following theorem. About the error in the approximate matrix A,
nb2
× nb2
respectively.
Theorem 2. Let a(x, y , u , v ) be a function defined in the domain [−1, 1]4 such that ∂ k a/∂ xk , ∂ k a/∂ yk , ∂ k a/∂ uk and ∂ k a/∂ v k are ˜ be the approximate matrix of A given by (3.4). continuous, where k is a positive integer. Let A (N )
(i) If all weights w i
are positive and formula (1.2) exactly holds for constant functions, then we have
( A˜ − A ) W ˜ = O n−k log4 n
for n > k, where · denotes the infinity matrix norm. (ii) If there exist two constants α and β satisfying 0 α β such that exactly holds for constant functions, then we have
α / N w (i N ) β/ N for i = 1, 2, . . . , N and formula (1.2)
( A˜ − A ) W ˜ = O n−k log4 n F
for n > k, where · F denotes the Frobenius norm. Proof. (i) By setting f (x) = 1 in (1.2) we get
N
l=1
˜ respectively. Note that A max
i , j ,l,m=1,..., N
(i , j ) (i , j ) A˜ −A = l,m
l,m
max
i , j ,l,m=1,..., N
max
(N )
wl
(i , j ) (i , j ) = 2. Let Al,m and A˜ l,m be (l, m)-entry of the (i , j )-block of A and
(N ) ( N ) (Cn a) t (N ) , t (N ) , t (N ) , tm − a t ( N ) , t ( N ) , t ( N ) , tm
−1x, y ,u , v 1
i
j
l
i
(C n a)(x, y , u , v ) − a(x, y , u , v ),
j
l
1716
W.-J. Xie, F.-R. Lin / Applied Numerical Mathematics 59 (2009) 1709–1719
by using (2.13) in Theorem 1, we get max
i , j ,l,m=1,..., N
(i , j ) (i , j ) A˜ −A = O n−k log4 n . l,m
l,m
Therefore N (i , j ) (i , j ) ( N ) ( N ) A˜ −A w wm
( A˜ − A ) W ˜ = max
1i ,l N
l,m
j
l,m
j ,m=1
max
1i , j ,l,mn
N (i , j ) (i , j ) (N ) (N ) A˜ − A w j wm l,m l,m j ,m=1
=
max
1i , j ,l,mn
N (i , j ) (i , j ) (N ) A˜ − A wj l,m l,m j =1
= O n−k log4 n . (ii) From
2
α / N W (j N ) β/ N, it follows [ w (jN ) ]2 β 2 / N 2 , j = 1, . . . , N. Therefore
( A˜ − A ) W ˜ 2 = F
N (i , j ) (i , j ) ( N ) ( N ) 2 A˜ −A w wm l,m
l,m
j
i , j ,l,m=1
max
1i , j ,l,mn
max
1i , j ,l,mn
=O n
−2k
(i , j ) (i , j ) 2 A˜ −A l,m
N
l,m
(N )
wj
( N ) 2
wm
i , j ,l,m=1 N (i , j ) 2 2 2 (i , j ) 2 2 A˜ − A N β /N l,m l,m j ,l=1
log n . 8
˜ − A)W ˜ F = O(n−k log4 n). That is, ( A
2
To end this subsection, we make the following remarks. (N )
1. For many quadrature rules, the conditions about W j ( j = 1, 2, . . . , N) are satisfied (for example, repeated low-order Newton–Cotes rules). ˜ is independent of the number of quadrature points N, that is, na can be chosen 2. We observe that the accurate in A a W ˜ )f = g will be very without considering the size of the linear system. If na is large enough, the solution of ( I − A a W ˜ )f = g. close to that of ( I − A W ˜ = O(n−k log4 nb ) will be small, and therefore the convergence of the RC 3. By choosing a moderate nb , ( A a − B a ) W b scheme (3.1) will be very fast. 3.2. Implementation and complexity analysis
˜ )−1 , which In this subsection, we first discuss how to compute the matrices P A and Λ A , P B and Λ B , and ( I − B a W are required in the iterative method (3.1) and then we discuss the implementation of the RC scheme. Details are listed as follows. The computation of Λ A and Λ B . The matrix Λ A can be obtained in two steps: 1) evaluate the function values of a(x, y , u , v ) at the Chebyshev grid to form the na2 × na2 matrix H g ; 2) apply two dimensional fast cosine transform on H g . The cost is O(na4 + na4 log na ) = O(na4 log na ). Similarly, the cost of obtaining Λ B is O(nb4 log nb ). The computation of P A and P B . Since that P A and P B are of size na × N and nb × N, the cost of obtaining them is O((na + nb ) N ). The computation of the residual vector. From (3.2) and (3.6), we get
r(q+1) = r(q) + I − ( P A ⊗ P A ) T Λ A ( P A ⊗ P A )( W ⊗ W ) d(q) . The computation of e = ( W ⊗ W )d(q) requires N 2 ops, the computation of ( P A ⊗ P A )e can be done in 2na N 2 ops by using ( P A ⊗ P A )e = vec( P A mat(e) P TA ), where mat(e) is the N × N matrix obtained by reordering the vector e (the first N entries form the first column of mat(e) and so on) and vec( A ) is the vector obtained by column-by-column ordering of A. Thus the cost of computing the residual vector is about N 2 + 4N 2 na + na4 ops (multiplications).
W.-J. Xie, F.-R. Lin / Applied Numerical Mathematics 59 (2009) 1709–1719
1717
A fast solution method for the preconditioning system. In the RC scheme (3.1), we require to solve the preconditioning system (3.3). We show in the following that by using the decomposition of B a , we can solve the equation efficiently. By using the Sherman–Morrison formula (see for instance [5, p. 74]), we see that
( I M − U V )−1 = I M + U ( I m − V U )−1 V , where I M and I m are identity matrices, and U and V are M × m and m × M matrices respectively. Note that
−1
I N 2 − B a (W ⊗ W )
−1 = I N 2 − ( P B ⊗ P B )T Λ B ( P B ⊗ P B )( W ⊗ W ) .
Let U = ( P B ⊗ P B ) T and V = Λ B ( P B ⊗ P B )( W ⊗ W ) = Λ B (( P B W ) ⊗ ( P B W )), we have
−1
I N 2 − B a (W ⊗ W )
−1 = I N 2 + U I n2 − Λ B P B W P BT ⊗ P B W P BT V. b
(3.7)
The matrix ( I n2 − Λ B (( P B W P BT ) ⊗ ( P B W P BT )))−1 can be obtained efficiently in O( N (nb + nb2 ) + nb6 ) ops by first computing b
P B W P BT (about N (nb + nb2 ) ops). The preconditioning system can be solved by using (3.7) in N 2 + 2N 2 nb + 2Nnb2 + 2nb4 ops.
The computational costs. In summary, the cost for constructing the required matrices is O(na4 log na + N (na + nb2 ) + nb6 ), the cost of the computation of the residual vector is O( N 2 na + na4 ), and the cost for solving the preconditioning system O( N 2 nb + Nnb2 + nb4 ). Storage requirement. We require to store the six matrices: P A , P B , Λ A , Λ B , W , and ( I − Λ B (( P B W P BT ) ⊗ ( P B W P BT )))−1 , which contain na N, nb N, na4 , nb4 , N and nb4 entries respectively. We also require to store three N 2 × 1 vectors r(q) , d(q) , f(q) . Therefore, the total Storage requirement is O( N 2 + Nna + na4 ). 4. Numerical examples In this section, we use our algorithms to solve one dimensional and two dimensional Fredholm equations of the second kind. We use the composite 4-point Gauss–Legendre quadrature to discretize the integral equations. Let the roots of the degree-4 Legendre polynomial be denoted by ξi (i = 1, 2, 3, 4), then the 4-point Gauss–Legendre quadrature is given by
1 f (x) dx =
4
ηi f (ξi ) + E 4 ( f ),
i =1
−1
where ηi (i = 1, 2, 3, 4) are the weights and E 4 ( f ) is the error in the quadrature. For the one dimensional equation
1 f (x) −
a(x, y ) f ( y ) dy = g (x),
−1 x 1,
−1
we partition the interval [−1, 1] into n subintervals of length d = 2/n. In this case, the weight matrix W is given by W =
d 2
I n ⊗ diag(η1 , η2 , η3 , η4 )
and the quadrature points are
(N )
ti
= −1 +
i−1 4
+
1 2
d+
(4.1)
d
ξ i −1 , 2 i −4[ 4 ]
i = 1, 2, . . . , N = 4n,
(4.2)
where [x] denotes the maximal integer that is not larger than x. For the two dimensional equation
1 1 f (x, y ) −
a(x, y , u , v ) f (u , v ) du dv = g (x, y ),
−1 x, y 1,
−1 −1
˜ = W ⊗ W with W given by (4.1) and the quadrature points are (t (N ) , t (N ) ), i , j = 1, 2, . . . , N, where the weight matrix is W i j (N )
for i = 1, 2, . . . , N are given by (4.2). The RC scheme is used to solve the discretization systems of one dimensional and two dimensional equation respectively. The initial guesses are set to zero vectors and the stopping criterion is
ti
(q) r = g − ( I − A a W )f(q) < 10−6 (for one dimensional equation)
or
(q) r = g − ( I − A a W ˜ )f(q) < 10−6 (for two dimensional equation),
1718
W.-J. Xie, F.-R. Lin / Applied Numerical Mathematics 59 (2009) 1709–1719
where f(q) is the approximate solution at the qth iteration. The error of the numerical solution f approx is computed by e ( f approx ) = f approx − f true ∞ . Two examples are tested. The first example is a one dimensional Fredholm equation:
1 f (x) −
a(τ ; x, t ) f (t ) dt = 1 − arctan
τ (1 − x) 2
−1
where
τ > 0 is a parameter and a(τ ; x, t ) = ∂ j a(τ ; x, t ) = ∂xj
j +1
τ
2
∂ j a(2; x, t ) ∂xj
− arctan
τ /2 1+(τ (x−t )/2)2
2
x ∈ [−1, 1],
,
(4.3)
. The true solution of (4.3) is f true (x) = 1 for x ∈ [−1, 1]. Since
∂ j a(τ ; x, t ) = ∂t j
and
τ (1 + x)
j +1
τ
2
∂ j a(2; x, t ) , ∂t j
a large τ will lead to large partial differentials of a(τ ; x, t ). Therefore, for a given increase in order that ( A a − A ) W < . The second example is a two dimensional Fredholm equation:
> 0, the value of na increases as τ
1 1 f (x, y ) −
a(τ ; x, y , u , v ) f (u , v ) du dv = g (τ ; x, y ),
(x, y ) ∈ [−1, 1] × [−1, 1],
(4.4)
−1 −1
where
τ > 0 is a parameter, a(τ ; x, y , u , v ) =
τ2 4
exp −
τ (1 + x) τ (1 + u ) 2
2
−
τ (1 + y ) τ (1 + v ) 2
2
and g (τ ; x, y ) = 1 −
2 2 2 e −τ (1+x+ y )/2 − e −τ (1+ y )/2 − e −τ (1+x)/2 + 1 .
4
τ 2 (1 + x)(1 + y )
The true solution of (4.4) is f true (x, y ) = 1. Similar to the one dimensional case, we have
∂ j a(τ ; t 1 , t 2 , t 3 , t 4 ) j ∂ ti
=
j +2
τ
∂ j a(2; t 1 , t 2 , t 3 , t 4 )
2
∂ ti
j
,
i = 1, 2, 3, 4,
and it follows that a large τ leads to a large na for a fixed accuracy . In Tables 1 and 2, we show the numbers of iterations which are denoted by #I and the errors in the numerical solutions which are denoted by Error for Eqs. (4.3) and (4.4) with different τ and N, respectively. We observe from the numerical results that: 1. The RC scheme converges very fast. 2. Given τ , for fixed na and nb , the numbers of iterations are the same for all N. This result is consistent with the ˜ depend only on na and nb . ˜ )−1 ( A a − B a ) W theoretical result that ( I − B a W )−1 ( A a − B a ) W and ( I − B a W 3. Given τ , for fixed na and nb , the error in the numerical solution will decrease as N increases when N is small, which indicates that the error in the quadrature rule is larger than that in the interpolating polynomial. When N is large enough, increasing the value of N will not improve the accuracy of the numerical solution, which indicates that the error in the interpolating polynomial is dominating in this case. Finally, we make a briefly remark on the relation of na , nb and N. In order that the error in the numerical solution decreases at N increase for large N, we should increase the value of na while fix the value of nb . Suppose the error in the numerical solution of (1.3), the discretization linear system of the original integral equation, is of order O( N −l ). Table 1 Numbers of iterations for the RC scheme and errors in the numerical solutions for one dimensional Fredholm equation (4.3) for the cases na = 32, nb = 8 (left), τ = 16, na = 128, nb = 32 (middle), and τ = 64, na = 320, nb = 120 (right). N
256 512 1024 2048 4096
(τ , na , nb ) = (4, 32, 8)
τ = 4,
(τ , na , nb ) = (16, 128, 32)
(τ , na , nb ) = (64, 320, 120)
#I
Error
#I
Error
#I
Error
6 6 6 6 6
1.80e−09 1.80e−09 3.23e−10 3.24e−10 3.23e−10
6 6 6 6 6
1.59e−09 1.08e−09 2.33e−10 6.29e−11 2.85e−10
6 6 6 6 6
1.56e−5 3.16e−7 9.12e−8 2.30e−8 1.68e−8
W.-J. Xie, F.-R. Lin / Applied Numerical Mathematics 59 (2009) 1709–1719
1719
Table 2 Numbers of iterations for the RC scheme and errors in the numerical solutions for two dimensional Fredholm equation (4.4) for the cases τ = 4, na = 16, nb = 6 (left) and τ = 8, na = 32, nb = 12 (right).
(τ , na , nb ) = (4, 16, 6)
N
32 64 128 256 512
(τ , na , nb ) = (8, 32, 12)
#I
Error
#I
Error
7 7 7 7 7
1.30e−06 1.01e−06 1.31e−06 1.30e−06 1.31e−06
8 8 8 8 8
1.63e−03 2.07e−05 6.74e−07 4.87e−07 6.74e−07
From Theorem 2 we see that the error in the approximate matrix is O(na−k log4 na ). Thus, in order that the error in the solution of (1.3) and that in the approximate matrix are of the same order, we require O( N −l ) = O(na−k log4 na ). Since that log4 na = o(na ) for large na , it follows that if we choose na such that
na = O N l/(k−1) ,
(4.5)
then the error in solution of (1.6) is also O( N −l ). Recall that the value of l depends on the quadrature rule used, and the smoothness of the solution and the kernel function. If the solution of (1.1) has some kind of singularity, then l and l/(k − 1) are small and it follows that na will increase slowly as N increases. On the other hand, if the solution of (1.1) is very smooth (as our examples) and a high order quadrature is used, then na will increase quickly as N increases. 5. Concluding remarks Based on interpolating polynomial approximation, we proposed new algorithms for the fast numerical solution of two dimensional Fredholm equation of the second kind with smooth kernel. According to (4.5) and the discussions in Section 3.2, we see that the complexity of our algorithm is: 1. Storage requirement: O( N 2 + Nna + na4 ) = O( N 2 + N (k+l−1)/(k−1) + N 4l/(k−1) ). 2. The cost per iteration of the RC scheme: 1) the cost of the computational of residual vector is N 2 + 4N 2 na + na4 = O( N (2k+l−2)/(k−1) + N 4l/(k−1) ); 2) the cost of the solution of the preconditioning system is N 2 + 2N 2 nb + 2Nnb2 + 2nb4 = O( N 2 ). We observe that the main work of the iteration is the computation of the residual vector. Numerical results were given to illustrate the efficiency of our methods. Acknowledgement We thank the anonymous referees for their detailed and valuable comments and suggestions. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21]
B. Alpert, G. Beylkin, R. Coifman, V. Rokhlin, Wavelets for the fast solution of second-kind integral equations, SIAM J. Sci. Comp. 14 (1993) 159–184. P.M. Anselone, Collectively Compact Operator Approximation Theory, Prentice-Hall, Englewood Cliffs, NJ, 1971. K.E. Atkinson, A Survey of Numerical Methods for the Solution of Fredholm Integral Equations of the Second Kind, SIAM, Philadelphia, PA, 1976. K.E. Atkinson, The Numerical Solution of Integral Equations of the Second Kind, Cambridge University Press, Cambridge, 1997. O. Axelsson, Iterative Solution Methods, Cambridge University Press, Cambridge, 1996. C.T. Baker, The Numerical Treatment of Integral Equations, Clarendon Press, Oxford, 1977. Z. Battles, L.N. Trefethen, An extension of MATLAB to continuous functions and operators, SIAM J. Sci. Comp. 25 (2004) 1743–1770. G. Beylkin, R. Coifman, V. Rokhlin, Fast wavelet transforms and numerical algorithms I, Comm. Pure Appl. Math. 46 (1991) 141–183. L. Brutman, Lebesgue functions for polynomial interpolation – a survey, Ann. Numer. Math. 4 (1977) 111–127. L.M. Delves, J.L. Mohamed, Computational Methods for Integral Equations, Cambridge University Press, New York, 1985. R. Farengo, Y.C. Lee, P.N. Guzdar, An electromagnetic integral equation: application to microtearing modes, Phys. Fluids 26 (1983) 3515–3523. G. Golub, C. Van Loan, Matrix Computations, second ed., Johns Hopkins University Press, Baltimore, 1989. L. Greengard, V. Rokhlin, A fast algorithm for particle simulations, J. Comput. Phys. 73 (1987) 325–348. L. Greengard, V. Rokhlin, A new version of the fast multipole method for the Laplace equation in three dimensions, Acta Numerica 6 (1997) 229–269. G.Q. Han, L.Q. Zhang, Asymptotic error expansion and extrapolation of the Nystrom eethord for two-dimensional Fredholm integral equations, Acta Math. Appl. Sinica 18 (2) (1995) 218–224 (in Chinese). L.C. Hsu, Y.S. Zhou, Y.B. Sun, Approximation Theory, National Defencial Industry Publishing Company, China, 1985 (in Chinese). V.V. Ivanov, V.K. Zadiraka, Some new results in theory of approximation, in: Computational Mathematics, Kiev University Press, Kiev, 1966, pp. 30–36 (in Russian). F.R. Lin, Preconditioned iterative methods for the numerical solution of Fredholm equations of the second kind, Calcolo 40 (2003) 231–248. J. Tausch, The variable order fast multipole method for boundary integral equations of the second kind, Computing 72 (2004) 267–291. Y. Wang, Y. Xu, A fast wavelet collocation method for integral equations on polygons, J. Integral Equations Appl. 17 (2005) 277–330. A.S. Xia, L. Wang, Iterative corrections for finite element solutions of two dimensional second kind Fredholm integral equations, J. Inst. Command Technology 11 (2) (2000) 105–108 (in Chinese).