- Email: [email protected]
Lp-approximation method for the numerical solution of singular integral equations
il NORTH-~
Lp-Approximation Method for the Numerical Solution of Singular Integral Equations J. H. De Klerk, D. Eyre, and L. M. Venter
Department of Mathematics and Applied Mathematics Potchefstroom University for Christian Higher Education Potchefstroom Republic of South Africa
Transmitted by Melvin Scott
ABSTRACT This paper describes a numerical approach for obtaining the solution of singular integral equations using an Lp-approximation method. The basic elements of this approach are a minimization with respect to a set of parameters in a function chosen to represent the solution of the integral equation and a numerical integration to evaluate the Lp-norm. Attention is given to different values of p > 1. The performance of the numerical procedures is investigated by solving a number of test problems.
1.
INTRODUCTION
T h e term "singular" in "singular integral equations" can be used rather widely. According to the use of Delves and M o h a m e d [1], it refers to any lack of analyticity in an integral equation problem. A m o n g others, these authors distinguish between the following types of singular integral equations: (i) integral equations with a semi-infinite or infinite range, (ii) integral equations with a discontinuous derivative in either the kernel or driving term, and (iii) integral equations with either an infinite or nonexisting derivative of some finite order. Of course, any combination of these phenom-
APPLIED MA THEMATICS AND COMPUTATION 72:285-300 (1995) © Elsevier Science Inc., 1995. 655 Avenue of the Americas, New York, NY 10010
0096-3003/95/$9.50 SSDI 0096-3003(95)00072-P
286
J.H. DE KLERK, D. EYRE, AND L. M. VENTER
ena may also occur. In this paper, attention will be given to problems exhibiting one or more of these "difficulties." Different numerical techniques for the solution of integral equations, including singular integral equations, are discussed, for example, in Atkinson [2] and in Delves and Mohamed [1]. Some of the more recent techniques for the solution of singular integral equations include those of Cuminato [3, 4] Diogo et al. [5], Kaneko and Xu [6], and Riley [7]. Most numerical techniques are based on collocation methods or are developed for some definite kind of singular integral equation. The Lp-approximation technique described in the present paper is related to well known methods in the sense that the minimization of the Le-norm of the residue is equivalent to a weighted Galerkin technique. Minimization of the residue in other norms will reduce to the Galerkin method only in special cases, however, and in general these norms will lead to a different result for the numerical solution on a given basis set. The aim of this study is to show how the Lp-approximation method can be used to solve singular integral equations. Emphasis is placed on the numerical treatment of singular integral equations. The theory of Lp-spaces is used to support the numerical procedures and to interpret the results of the Lp-approximation method. We begin with a brief discussion of the relevant facts concerning the theory of Lp-spaces. (For a thorough discussion of the theory of Lp-spaces consult the monograph by KSthe [8].) Consider Lebesgue measure on the subset X = [ a, b] of the real line. The space Lp(X), 1 < p < ~, is defined as the set of all measurable real (or complex) valued functions f defined on X (except possibly on a subset of measure zero) such that r f(t)l p is integrable over X. To make Lp(X) into a seperable complete normed vector space, functions that are equal up to a set of measure 0 are regarded as the same function. The norm is defined by 1
To define the space Loo(X), again consider the measurable functions f defined on X (except possibly on a subset of measure zero). One says that f is essentially bounded if Jf(t)l is bounded from above for all t outside a set of measure 0. The essential supremum of f is the smallest of these upper bounds, denoted by ess sup If(t)f. L~(X) consists of all essentially bounded measurable functions defined on X, where again one identifies functions that are equal almost everywhere. If a norm on L~(X) is defined by IIflL = ess supl f ( t) l,
(1.2)
Lp-Approximation Method
287
then it is a complete normed vector space, but it is not seperable as in the case of Lp(X), 1 _< p < ~. From HSlder's inequality, 1
Ilfglll ~ Ilfllpllgllq,
1
- + - = 1, p q
p > 1,
one can deduce that:
THEOREM.
If 1 < Pl
<
P2 < oo, then L~ c Lp~ c Lpl C L1, and if f
Lp2, then 1
Ilfllpl ~< Ilftlp~(b- a) pl /ff~
1
p~
(1.3)
L~(X), then Ilfll~ = p---* limc~ IIfL[p.
If b - a = 1, one has a particularly interesting case; the sequence (11fll p)~= 1 is then monotonically increasing, and in fact IIfll p 1' IIfll~. The present paper is organized as follows: In Section 2 the Lp-approximation method for solving integral equations is discussed; in Section 3 attention is given to a number of test problems that exhibit singular behavior in the kernel and in the solution of the integral equation, as well as to a real physical problem (in this case from the field of quantum scattering theory); in Section 4 the results are discussed, and in Section 5 some concluding remarks are made. 2.
T H E L,-APPROXIMATION METHOD
This section describes the Lp-approximation method for solving integral equations. The discussion here will be restricted to linear integral equations of the second kind; however, the same procedure may be applied to more general integral equations, including integral equations of the first kind. Consider the integral equation
x(t) = y(t) - fabK(t,s) x(s)ds,
t ~ [a,b].
(2.1)
288
J . H . DE KLERK, D. EYRE, AND L. M. VENTER
Here the kernel K and driving function y are known, and x is to be calculated. In general the kernel is not continuous; this is the case for integral equations with singular kernels. Let (~h0, ~bl,.-- , ~hN) be a set of linearly independent trial functions. The solution x(t) is approximated in some way which can be either linear or nonlinear. For example, in the case of a linear combination, by N
F( A, t) --- E a i 0 ~ ( t ) ,
(2.2)
i=0
where F( A, t) depends on the N + 1 coefficients A = (a0, al,..., r( A, t) be the residual function r(A,t)
= ]~
0,(t) +
-
O~N).
Let
,
i=0
(2.3) Clearly a good choice of the coefficients, A, is one that will make the residual small in some sense. The idea behind the Lp-approximation method is to minimize the Lp-norm of the residual function with respect to the coefficients A. The approximation problem then becomes 1
fblT"( A, t)] p dt) -~.
min(A\aa
(2.4)
The approximate solution of the integral equation, ~(t) = F( A, t), is obtained from the coefficients A resulting from the minimization problem. The two major steps in the numerical technique are (i) a numerical search procedure to calculate the coefficients .A, and (ii) a numerical integration procedure to evaluate the Lp-norm. Some additional information on these procedures will be given in the next section. The following general remarks can be made concerning these procedures: • The initial values in the minimization procedure can be chosen from previously calculated values using smaller N. • When moving from a smaller to larger value of p, the solution, ,4, calculated for the smaller p provides a good initial value in the minimization procedure for the larger p. • For integration over absolute value, an adaptive integration procedure should be used so that the zeros can be determined numerically.
L p-A pprozirnation Method
289
In De Klerk [9, 10] some aspects of this numerical technique for non-singular integral equations are discussed in more detail for the special case p = 1; among others, the ease of a nonlinear integral equation and of bifurcations in the solution of the integral equation. We remark that, for linear integral equations of the second kind and for linear combinations of basis functions, the minimization of the L2-norm of the residue function is equivalent to the Petrov-Galerkin method
E °=
t) ~,(t) dt +
% ( t ) K( t, ~),~(~) d~ dt ~,
= f ) % ( t ) v(t) et,
(2.~)
where the test functions, ~k(t), have the special form
% ( t ) = ~o~(t) + fa%~( S) K( ~, t) e~. 3.
(2.6)
NUMERICAL EXAMPLES
As a demonstration of the above-mentioned technique, attention will first be given to some test problems. As a final example, a real physical problem will be considered. In all cases different values of p will be employed in the approximation method. In constructing the approximating function, it is important to make a good choice for the trial functions in the sense that these functions should have the same characteristics, as far as possible, as the unknown solution of the original integral equation. In practice this can be a problem. A reasonable choice are functions that can represent, for example, the driving term or some known a s y m p t o t i c behavior of the analytical solution. With this in mind an appropriate choice of trial functions will be made for each of the test examples. The numerical integration procedure used for evaluating integrals over absolute values of functions is performed via a technique developed by Laurie and Venter [11]. For this technique the standard adaptive integration routine QAGE (from the Quadpack subroutine package [12]) was modified to allow for the efficient evaluation of integrals of the general form f~[ f( x)l dz. Using function values available during the integration process, the zeros of the fimction are found and used as subdivision points.
290
J.H. DE KLERK, D. EYRE, AND L. M. VENTER
The minimization was performed using standard optimization algorithms (UMPOL in the IMSL Library and the downhill simplex method of Nelder and Mead (see, for example, Press et al. [13])). The results were insensitive to the choice of these optimization algorithms.
EXAMPLE 1. The first example is a singular integral equation with a weak singularity in the kernel. The equation (Kaneko and Xu [6]) is 1 x(s) ds
y(t),
(3.1)
with
y(t) = 3[t(1 - 0 ] 3/4
3
-Ir[1 + 4t(1
-
8
-
t)],
(3.2)
and the solution x(t) = 2v~[t(1
-
0] 3/4 ,
(3.3)
for t ~ [0, 1]. (Note that reference [6] has a misprint that has been corrected in the above expression.) Setting as the approximating function, N
x( t) ~ F( A, t) = ~., a,t i,
(3.4)
i=O
the residual function becomes
r(A,t)
= Ea i
- 4(t)
- y(t),
(3.5)
i=0
with
Ji( t) --- f l J0
sid8
Iv17- sl
(3.6)
L p-Approximation Method
291
Two different techniques can be used for calculating the integral in (3.6)--the "inner" integral--namely, (a) a numerical, or (b) an analytical technique. It should be noted that it would generally not be possible to calculate the inner integral using analytical methods. However, whenever it is possible to do the integration analytically, this can often lead to a considerable saving in computation time. In this particular example the latter technique is employed. For Ji(t), one finds the recurrence relation J o ( t ) = 2(vr2 + ~/1 - t ) , 2
J'(tt
2 i + 1 (itJ~-l(t) + l~-C~-t), i = 1 . . . . . g.
(3.7 t
The approximation problem then becomes 1
rain
t),P dt)
(3.8)
Calculating the minimum value of the above approximation problem, the best results in the Lp-norm for different values of p and for the above choice of approximation function are shown in Tables l a and lb.
EXAMPLE 2. The second example is a singular integral equation with a strong singularity in the kernel. In this case the integral equation (Cuminato [3]) is 1 ~
y(t)=~
1+
1 x ( s ) ds 1 s t
y(t),
(3.9t
Vu~_t21n v ~ _ t 2 + t - 1
'
TABLE la CALCULATED COEFFICIENTS OF THE APPROXIMATION FUNCTION
N
p= 1
p= 2
p= 4
p=8
p = 16
1
0.729875 0.000001
0.743128 0.000000
0.781953 0.000000
0.823502 0.000000
0.853554 0.000000
2
0.073169 3.868774 -3.868776
0.037902 4.067628 -4.067625
0.005074 4.249036 -4.249036
- 0.013321 4.347200 -4.347200
- 0.021027 4.386523 -4.386523
292
J. H. DE KLERK, D. EYRE, AND L. M. VENTER TABLE lb CALCULATED MINIMUM Lp-NORIvlS OF THE RESIDUAL FUNCTION
N
p= 1
1 2
0.036689 0.005923
p=2
p=4
0.079033 0.008215
0.134422 0.010976
p= 8
p = 16
0.189447 0.013203
0.228436 0.014582
where t ~ [ - 1, 1] and the integral in (3.9) is defined in the principal-value sense. This integral equation has the solution
tltl x(t) = ~
,
(3.11)
if the additional condition f 1__1 x ( t ) dt = 0 is imposed. For this example the underlying approximation is chosen to be a rational (i.e., a nonlinear) function. The solution can then be represented by a function with simple poles, thus providing an improved numerical solution close to the singularities of the analytical solution. The approximation function is the rational function, n
Ei=oOli t
i
x( t) = F( A , t) ~- F~"i=ofli t i ' with fl0 = 1. This function has
N=
m+
(3.12)
1 coefficients with
n+
A =
(ao, a 1. . . . . a , , fll . . . . . tim). T h e residual function then becomes
~( A,
1 1 1 Ei~=o%s ' t) = -~[l~- s - t E:~ot3~s ~ ds - y ( t ) .
(3.13)
Again calculating the inner integral analytically, and considering the case m = n = 2, one arrives at the residual function,
r(A,t)
= --
1[
~-/32
A*ln
+ B*ln ~
blXl
1 + X1
1-x ]I -
+ c* in ~
y(t)'
(3.14)
Lv-Approximation Method
293
with
X, =
2/32
X2 =
2&
,
(3.15(a))
,
(3.15(b))
a o + a~ t + a 2 t 2
A* = ( t -
X,)(t-
X2)'
(a.15(e))
Ot0 + Ot1 X 1 + a 2 X 2 B* =
C* =
(3.15(d))
( x~ - t ) ( x~ - x ~ ) ' a o + a 1 X2 + a 2 X 2 (X2 _ t)(X 2 - X~)"
(3.15(e))
T o d e t e r m i n e t h e coefficients, A, for a n i n t e g r a l w i t h l e n g t h 1 ( c o m p a r e t h e r e m a r k for t h e special case b - a = 1 m e n t i o n e d a f t e r t h e t h e o r e m in S e c t i o n 1), t h e v a r i a b l e t ~ [ - 1, 1] is m a p p e d o n t o u ~ [ - 1, !]2 u s i n g t h e t t r a n s f o r m a t i o n u = 3. T h e L K a p p r o x i m a t i o n p r o b l e m t h e n b e c o m e s 1
r( A, u)] p du
min m
(3.16)
.
__ 2
C a l c u l a t i n g t h e m i n i m u m v a l u e of t h e a b o v e a p p r o x i m a t i o n p r o b l e m , t h e b e s t r e s u l t s in t h e L K n o r m for d i f f e r e n t v a l u e s of p a n d for t h e a b o v e c h o i c e of a p p r o x i m a t i o n f u n c t i o n a r e g i v e n in T a b l e s 2a a n d 2b.
T A B L E 2a CALCULATED COEFFICIENTSOF THE APPROXIMATIONFUNCTION
p=l
p=2
p=4
p=8
p=16
- 0.015901 0.499960 0.009663 0.000199 - 0.952559
0.000132 0.457612 -0.000090 -0.000002 - 0.961707
- 0.000002 0.405397 0.000001 0.000000 - 0.968475
-0.057666 0.374331 0.043482 0.000682 -0.971525
- 0.012619 0.366593 0.009586 0.000144 -0.971796
294
J. H. DE KLERK, D. EYRE, AND L. M. VENTER
TABLE 2b CALCULATED MINIMUM /p-NORMS OF THE RESIDUAL FUNCTION
p=l
p=2
p=4
p=8
p=16
0.057091
0.090140
0.138462
0.174507
0.189342
EXAMPLE 3. T h e last example is a singular integral equation from the theory of q u a n t u m scattering (Reid [14])
m( q, K ) = v( q, K ) - 2 ~ v (
q, ~/ ) m( ¢, K ) - -
at2 _ /£2"
(3.17)
The p a r a m e t e r K is defined by the scattering energy
2 = 27rM > 0,
(3.18)
where M is the particle mass and h is Planck's constant. T h e function v( q, q') is the Fourier transformed potential, while m(q, K) is a real valued function, normalized such t h a t m( K, K) = - [K cot 8( K)] - 1 ,
(3.19)
where 8(K) is the scattering phase shift. In nuclear scales the ratio h2/(2 ~-M) = 41.47 MeV.fm 2. The calculations are carried out at a scattering energy e = 24 MeV. T h e potential function v(q, q') is chosen to have a form t h a t is typical of phenomenological two-nucleon potentials. For this purpose the Fourier transform of the Reid potential [14] is taken,
v( q, ~ )
~M 3 2 h2~l q~/ ~., V~ in i=1
( q + 4) 2 + (q_
] (3.20)
where ~[~1 ---~ 0 . 7 f m - 1 , ~$.L2 = 4#1 , IZ3 = 7/zl, V1 = - 1 0 . 4 6 3 MeV.fm -3, V2 = - 1650.6 MeV.fm -3, and V3 = 6484.2 MeV.fm -3. This last example represents a realistic problem t h a t requires a numerical procedure to evaluate the integrals. A detailed discussion of the approximation functions and the numerical procedures is now given.
Lp-Approximation Method
295
The trial functions are chosen as Chebyshev polynomials of the second kind, U~(t). In the q-variable the trial functions are obtained by mapping t onto the semi-infinite domain 1
q = ~" - ] - - ~ ]
,
(3.21)
where ~ > 0 is some appropriately chosen mapping parameter. The value of remains fixed throughout the calculation. This parameter was chosen to be ~"= 5. In the mapped variable the trial functions become 3
~i(q) = ( 2 ~ ) ~ (q2 + ~2) Ui q2 + ~.2 •
(3.22)
The trial functions satisfy the orthogonality relation
O,(q)~bj(q) q2 + ~2 dq = 6,j,
(3.23)
and behave asymptotically like q-:. This asymptotic behavior of the functions qJ, makes the trial functions particularly well suited for approximating the solution of the quantum scattering integral equation. For this example the approximation function is the linear combination N
m(q,K) = F( A(K),q) = ~_.a,(K)~O,(q).
(3.24)
i=0
The principal value singularity is evaluated using the method of subtraction. Making use of the identity
f
d4
0,
(3.25)
296
J.H. DE KLERK, D. EYRE, AND L. M. VENTER
the residual function becomes
N[
r ( A , q ) = ~ ~O~(q) 4=0
+ --Tr
q~2
K2
× a i ( K ) - v(q, K).
(3.26)
The integral is evaluated numerically using a standard Gauss-Legendre quadrature formula. The points and weights are mapped onto a semi-infinite domain q' ~ [0, ~) yielding the set {q,, w~}L=] where L is the number of quadrature points. An additional point qL+l K is included in this set of points. The residual function is now approximated by the double summation =
2 L+I ~b~(q) + -
~ ( A , q ) = ~] i=O
17"
]
F, v(q, q)~,(Owj ~,(,,1- v(q,K), 3=1 (3.27)
where qj2coj 2 K2, qj--
j= 1,...,L, (3.28)
w~ =
L
K2(.Ok
-Fl= qkff---2'K j=L+I To determine the coefficients, A, the variable q is mapped onto the t The final step is interval [ - 1, ½] using the transformation u = 7.
1 min A
]~( A, u)l p du
,
(3.29)
L p-Approximation Method
297
to obtain the coefficients A and the approximate solution of (3.17), N
q, ,,) = E
(3.30)
i=0
Calculating this m i n i m u m value, the best results in the Lp-norm for different values of p and for the mentioned choice of approximation function are given in Tables 3a, 3b, and 3c. 4.
DISCUSSION OF RESULTS
T h e results given in the previous section represent the " b e s t " results for a given approximating function in the sense t h a t they produce the smallest residual function in the Lp-norm. For a small n u m b e r of trial functions it can be seen t h a t each Lp-approximation produces a different solution. T h e usual p-norms such as p = 1, 2 and ~ will normally suffice for most solutions; however, other values of p, namely 2 < p < ~, can be used to obtain approximations for the L~-solution where this is difficult to obtain. Some general features of the approximation for a given p can be seen in the numerical results. For example, it is noted in Example 1, for the case p = 1, N = 2, t h a t the absolute value of the greatest difference between the analytical solution and the Ll-approximate solution is 0.073169 at the end
TABLE 3a RESULT USING N = 2 TRIALFUNCTIONSSHOWINGTHE CALCULATED PHASESHIFT (~AND Lp-NORMOF THE RESIDUALFUNCTION
1 2 3 4 5 10 15 20 40 80 160 :¢
26.45 31.72 32.78 33.53 33.60 34.54 34.82 34.87 35.10 35.22 35.24 35.35
0.092 O.114 0.125 0.133 0.139 0.153 0.159 0.163 0.170 0.175 0.177 0.180
298
J . H . DE KLERK, D. EYRE, AND L. M. VENTER TABLE 3b CALCULATED PHASE SHIFT ~ USING NTRIAL FUNCTIONS
N
p=l
p=2
p=3
p=4
2 3 4 5 6 7
26.45 39.12 38.58 38.98 39.05 38.84
31.72 38.68 38.73 38.70 38.70 38.84
32.78 38.39 38.75 38.70 38.74 38.85
32.53 38.22 38.77 38.71 38.74 39.02
points. This relatively large error at the boundaries of the interval is typical of approximations in the Ll-norm. Disregarding these extremities, and considering the interval [0.01, 0.99], the greatest difference is only about 0.04. In general the approximating function used in Example 2 can have poles on the real line, or in the complex plane close to the real line, that are not present in the analytical solution. The reason for choosing m - - - 2 is to produce a pole structure in the approximate solution that simulates the end-point singularities in the analytical solution. In all the examples, the values calculated for the Lp-norms of the residual function increase monotonically as p increases. This is consistent with the theory of Lp-spaces; under the conditions imposed in the various examples (the length of the interval of integration is 1) the Lv-norm of a fixed function would increase as p increases. Because the Lp-norm is less than the L~-norm, one can also reasonably expect the norm of the L~-residual (if it exists) to be larger than the norm of the Lp-residual. This is illustrated by Example 3. In fact, from Table 3a it is clear that the norms of the residuals form an increasing sequence that converges to the norm of the L~-residual.
TABLE 3c CALCULATED MINIMUM Lp-NORMS OF THE RESIDUAL FUNCTION
N 2 3 4 5 6 7
p=l 0.09195 0.02141 0.00855 0.00305 0.00285 0.00184
p=2 0.11399 0.02526 0.00969 0.00521 0.00509 0.00268
p=3 0.12542 0.02731 0.01040 0.00638 0.00637 0.00373
p=4 0.13313 0.02867 0.01087 0.00727 0.00723 0.00447
L p-Approximation Method
299
In Example 2 the analytical solution is not in L~, because it is unbounded near the end points of the interval. Even in this case, an approximate solution might exist in the L~-sense. If such a solution exists, the norm of the residual will be larger than the values calculated for large p.
5.
CONCLUSIONS
An important feature of the Lp-approximation method is a good choice of approximating function. The choice of approximating function is guided by general information that can be obtained about the solution of the integral equation such as its singular or asymptotic behavior. Where a good choice of the approximating function is possible, it is found that the Lp-approximation method is an efficient and accurate method for solving the integral equation. This paper has shown using examples of integral equations that exhibit different singular behavior, that solutions of singular integral equations can be obtained from the Lp-approximation method via a technique developed previously. The relationship between Lp-norms allows bounds to be placed on the residual function for different values of p. The Lp-approximation method can be used to construct solutions in the L~-norm even when these solutions are difficult to obtain directly by L~-methods. REFERENCES 1 L.M. Delves and J. L. Mohamed, Computational Methods for Integral Equations, Cambridge University Press, Cambridge, 1985. 2 K.E. Atkinson, A Survey of Numerical Methods for the Solution of Fredholm Integral Equations of the Second Kind, SIAM, Philadelphia, 1978. 3 J.A. Cuminato, On the uniform convergence of a collocation method for a class of singular integral equations, BIT 27:190-202 (1987). 4 J.A. Cuminato, Uniform convergence of a collocation method for the numerical solution of Cauchy-type singular integral equations: a generalization, IMA J. Numer. Anal. 12:31-45 (1992). 5 T. Diogo, S. McKee, and T. Tang, A Hermite-type collocation method for the solution of an integral equation with a certain weakly singular kernel, IMA J. Numer. Anal. 11:595-605 (1991). 6 H. Kaneko and Y. Xu, Numerical solutions for weakly singular Fredholm integral equations of the second kind, Appl. Numer. Math. 7:167-177 (1991). 7 B.V. Riley, The numerical solution of Volterra integral equations with nonsmooth solutions based on sinc approximation, Appl. Numer. Anal. 9:249-257 (1992). 8 G. KSthe, Topological Vector Spaces I, Springer-Verlag, Berlin, 1969.
300
J . H . DE KLERK, D. EYRE, AND L. M. VENTER
9 J.H. De Klerk, Solutions of integral equations via Ll-approximations, Int. J. Control 48:2121-2128 (1988). 10 J.H. De Klerk, The sensitivity of a computational Ll-approximation technique for integral equations, Appl. Math. Comput. 62:171-182 (1994). 11 D.P. Laurie and A. Venter, Automatic Quadrature of Functions of the Form If(x)l, Technical note T/15, Potchefstroom University, 1993. 12 R. Piessens, E. De Doncker-Kapenga, C. W. Uberhuber, and D. K. Kahaner, QUADPACK--A Subroutine Package for Automatic Integration, SpringerVerlag, Berlin, 1983. 13 W.H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling, Numerical Recipes in Pascal: The Art of Scientific Programming, Cambridge University Press, Cambridge, 1989. 14 R. V. Reid, Local phenomenological nucleon-nucleon potentials, Ann. Phys. 50:411-448 (1968).
Lp-Approximation Method for the Numerical Solution of Singular Integral Equations J. H. De Klerk, D. Eyre, and L. M. Venter
Department of Mathematics and Applied Mathematics Potchefstroom University for Christian Higher Education Potchefstroom Republic of South Africa
Transmitted by Melvin Scott
ABSTRACT This paper describes a numerical approach for obtaining the solution of singular integral equations using an Lp-approximation method. The basic elements of this approach are a minimization with respect to a set of parameters in a function chosen to represent the solution of the integral equation and a numerical integration to evaluate the Lp-norm. Attention is given to different values of p > 1. The performance of the numerical procedures is investigated by solving a number of test problems.
1.
INTRODUCTION
T h e term "singular" in "singular integral equations" can be used rather widely. According to the use of Delves and M o h a m e d [1], it refers to any lack of analyticity in an integral equation problem. A m o n g others, these authors distinguish between the following types of singular integral equations: (i) integral equations with a semi-infinite or infinite range, (ii) integral equations with a discontinuous derivative in either the kernel or driving term, and (iii) integral equations with either an infinite or nonexisting derivative of some finite order. Of course, any combination of these phenom-
APPLIED MA THEMATICS AND COMPUTATION 72:285-300 (1995) © Elsevier Science Inc., 1995. 655 Avenue of the Americas, New York, NY 10010
0096-3003/95/$9.50 SSDI 0096-3003(95)00072-P
286
J.H. DE KLERK, D. EYRE, AND L. M. VENTER
ena may also occur. In this paper, attention will be given to problems exhibiting one or more of these "difficulties." Different numerical techniques for the solution of integral equations, including singular integral equations, are discussed, for example, in Atkinson [2] and in Delves and Mohamed [1]. Some of the more recent techniques for the solution of singular integral equations include those of Cuminato [3, 4] Diogo et al. [5], Kaneko and Xu [6], and Riley [7]. Most numerical techniques are based on collocation methods or are developed for some definite kind of singular integral equation. The Lp-approximation technique described in the present paper is related to well known methods in the sense that the minimization of the Le-norm of the residue is equivalent to a weighted Galerkin technique. Minimization of the residue in other norms will reduce to the Galerkin method only in special cases, however, and in general these norms will lead to a different result for the numerical solution on a given basis set. The aim of this study is to show how the Lp-approximation method can be used to solve singular integral equations. Emphasis is placed on the numerical treatment of singular integral equations. The theory of Lp-spaces is used to support the numerical procedures and to interpret the results of the Lp-approximation method. We begin with a brief discussion of the relevant facts concerning the theory of Lp-spaces. (For a thorough discussion of the theory of Lp-spaces consult the monograph by KSthe [8].) Consider Lebesgue measure on the subset X = [ a, b] of the real line. The space Lp(X), 1 < p < ~, is defined as the set of all measurable real (or complex) valued functions f defined on X (except possibly on a subset of measure zero) such that r f(t)l p is integrable over X. To make Lp(X) into a seperable complete normed vector space, functions that are equal up to a set of measure 0 are regarded as the same function. The norm is defined by 1
To define the space Loo(X), again consider the measurable functions f defined on X (except possibly on a subset of measure zero). One says that f is essentially bounded if Jf(t)l is bounded from above for all t outside a set of measure 0. The essential supremum of f is the smallest of these upper bounds, denoted by ess sup If(t)f. L~(X) consists of all essentially bounded measurable functions defined on X, where again one identifies functions that are equal almost everywhere. If a norm on L~(X) is defined by IIflL = ess supl f ( t) l,
(1.2)
Lp-Approximation Method
287
then it is a complete normed vector space, but it is not seperable as in the case of Lp(X), 1 _< p < ~. From HSlder's inequality, 1
Ilfglll ~ Ilfllpllgllq,
1
- + - = 1, p q
p > 1,
one can deduce that:
THEOREM.
If 1 < Pl
<
P2 < oo, then L~ c Lp~ c Lpl C L1, and if f
Lp2, then 1
Ilfllpl ~< Ilftlp~(b- a) pl /ff~
1
p~
(1.3)
L~(X), then Ilfll~ = p---* limc~ IIfL[p.
If b - a = 1, one has a particularly interesting case; the sequence (11fll p)~= 1 is then monotonically increasing, and in fact IIfll p 1' IIfll~. The present paper is organized as follows: In Section 2 the Lp-approximation method for solving integral equations is discussed; in Section 3 attention is given to a number of test problems that exhibit singular behavior in the kernel and in the solution of the integral equation, as well as to a real physical problem (in this case from the field of quantum scattering theory); in Section 4 the results are discussed, and in Section 5 some concluding remarks are made. 2.
T H E L,-APPROXIMATION METHOD
This section describes the Lp-approximation method for solving integral equations. The discussion here will be restricted to linear integral equations of the second kind; however, the same procedure may be applied to more general integral equations, including integral equations of the first kind. Consider the integral equation
x(t) = y(t) - fabK(t,s) x(s)ds,
t ~ [a,b].
(2.1)
288
J . H . DE KLERK, D. EYRE, AND L. M. VENTER
Here the kernel K and driving function y are known, and x is to be calculated. In general the kernel is not continuous; this is the case for integral equations with singular kernels. Let (~h0, ~bl,.-- , ~hN) be a set of linearly independent trial functions. The solution x(t) is approximated in some way which can be either linear or nonlinear. For example, in the case of a linear combination, by N
F( A, t) --- E a i 0 ~ ( t ) ,
(2.2)
i=0
where F( A, t) depends on the N + 1 coefficients A = (a0, al,..., r( A, t) be the residual function r(A,t)
= ]~
0,(t) +
-
O~N).
Let
,
i=0
(2.3) Clearly a good choice of the coefficients, A, is one that will make the residual small in some sense. The idea behind the Lp-approximation method is to minimize the Lp-norm of the residual function with respect to the coefficients A. The approximation problem then becomes 1
fblT"( A, t)] p dt) -~.
min(A\aa
(2.4)
The approximate solution of the integral equation, ~(t) = F( A, t), is obtained from the coefficients A resulting from the minimization problem. The two major steps in the numerical technique are (i) a numerical search procedure to calculate the coefficients .A, and (ii) a numerical integration procedure to evaluate the Lp-norm. Some additional information on these procedures will be given in the next section. The following general remarks can be made concerning these procedures: • The initial values in the minimization procedure can be chosen from previously calculated values using smaller N. • When moving from a smaller to larger value of p, the solution, ,4, calculated for the smaller p provides a good initial value in the minimization procedure for the larger p. • For integration over absolute value, an adaptive integration procedure should be used so that the zeros can be determined numerically.
L p-A pprozirnation Method
289
In De Klerk [9, 10] some aspects of this numerical technique for non-singular integral equations are discussed in more detail for the special case p = 1; among others, the ease of a nonlinear integral equation and of bifurcations in the solution of the integral equation. We remark that, for linear integral equations of the second kind and for linear combinations of basis functions, the minimization of the L2-norm of the residue function is equivalent to the Petrov-Galerkin method
E °=
t) ~,(t) dt +
% ( t ) K( t, ~),~(~) d~ dt ~,
= f ) % ( t ) v(t) et,
(2.~)
where the test functions, ~k(t), have the special form
% ( t ) = ~o~(t) + fa%~( S) K( ~, t) e~. 3.
(2.6)
NUMERICAL EXAMPLES
As a demonstration of the above-mentioned technique, attention will first be given to some test problems. As a final example, a real physical problem will be considered. In all cases different values of p will be employed in the approximation method. In constructing the approximating function, it is important to make a good choice for the trial functions in the sense that these functions should have the same characteristics, as far as possible, as the unknown solution of the original integral equation. In practice this can be a problem. A reasonable choice are functions that can represent, for example, the driving term or some known a s y m p t o t i c behavior of the analytical solution. With this in mind an appropriate choice of trial functions will be made for each of the test examples. The numerical integration procedure used for evaluating integrals over absolute values of functions is performed via a technique developed by Laurie and Venter [11]. For this technique the standard adaptive integration routine QAGE (from the Quadpack subroutine package [12]) was modified to allow for the efficient evaluation of integrals of the general form f~[ f( x)l dz. Using function values available during the integration process, the zeros of the fimction are found and used as subdivision points.
290
J.H. DE KLERK, D. EYRE, AND L. M. VENTER
The minimization was performed using standard optimization algorithms (UMPOL in the IMSL Library and the downhill simplex method of Nelder and Mead (see, for example, Press et al. [13])). The results were insensitive to the choice of these optimization algorithms.
EXAMPLE 1. The first example is a singular integral equation with a weak singularity in the kernel. The equation (Kaneko and Xu [6]) is 1 x(s) ds
y(t),
(3.1)
with
y(t) = 3[t(1 - 0 ] 3/4
3
-Ir[1 + 4t(1
-
8
-
t)],
(3.2)
and the solution x(t) = 2v~[t(1
-
0] 3/4 ,
(3.3)
for t ~ [0, 1]. (Note that reference [6] has a misprint that has been corrected in the above expression.) Setting as the approximating function, N
x( t) ~ F( A, t) = ~., a,t i,
(3.4)
i=O
the residual function becomes
r(A,t)
= Ea i
- 4(t)
- y(t),
(3.5)
i=0
with
Ji( t) --- f l J0
sid8
Iv17- sl
(3.6)
L p-Approximation Method
291
Two different techniques can be used for calculating the integral in (3.6)--the "inner" integral--namely, (a) a numerical, or (b) an analytical technique. It should be noted that it would generally not be possible to calculate the inner integral using analytical methods. However, whenever it is possible to do the integration analytically, this can often lead to a considerable saving in computation time. In this particular example the latter technique is employed. For Ji(t), one finds the recurrence relation J o ( t ) = 2(vr2 + ~/1 - t ) , 2
J'(tt
2 i + 1 (itJ~-l(t) + l~-C~-t), i = 1 . . . . . g.
(3.7 t
The approximation problem then becomes 1
rain
t),P dt)
(3.8)
Calculating the minimum value of the above approximation problem, the best results in the Lp-norm for different values of p and for the above choice of approximation function are shown in Tables l a and lb.
EXAMPLE 2. The second example is a singular integral equation with a strong singularity in the kernel. In this case the integral equation (Cuminato [3]) is 1 ~
y(t)=~
1+
1 x ( s ) ds 1 s t
y(t),
(3.9t
Vu~_t21n v ~ _ t 2 + t - 1
'
TABLE la CALCULATED COEFFICIENTS OF THE APPROXIMATION FUNCTION
N
p= 1
p= 2
p= 4
p=8
p = 16
1
0.729875 0.000001
0.743128 0.000000
0.781953 0.000000
0.823502 0.000000
0.853554 0.000000
2
0.073169 3.868774 -3.868776
0.037902 4.067628 -4.067625
0.005074 4.249036 -4.249036
- 0.013321 4.347200 -4.347200
- 0.021027 4.386523 -4.386523
292
J. H. DE KLERK, D. EYRE, AND L. M. VENTER TABLE lb CALCULATED MINIMUM Lp-NORIvlS OF THE RESIDUAL FUNCTION
N
p= 1
1 2
0.036689 0.005923
p=2
p=4
0.079033 0.008215
0.134422 0.010976
p= 8
p = 16
0.189447 0.013203
0.228436 0.014582
where t ~ [ - 1, 1] and the integral in (3.9) is defined in the principal-value sense. This integral equation has the solution
tltl x(t) = ~
,
(3.11)
if the additional condition f 1__1 x ( t ) dt = 0 is imposed. For this example the underlying approximation is chosen to be a rational (i.e., a nonlinear) function. The solution can then be represented by a function with simple poles, thus providing an improved numerical solution close to the singularities of the analytical solution. The approximation function is the rational function, n
Ei=oOli t
i
x( t) = F( A , t) ~- F~"i=ofli t i ' with fl0 = 1. This function has
N=
m+
(3.12)
1 coefficients with
n+
A =
(ao, a 1. . . . . a , , fll . . . . . tim). T h e residual function then becomes
~( A,
1 1 1 Ei~=o%s ' t) = -~[l~- s - t E:~ot3~s ~ ds - y ( t ) .
(3.13)
Again calculating the inner integral analytically, and considering the case m = n = 2, one arrives at the residual function,
r(A,t)
= --
1[
~-/32
A*ln
+ B*ln ~
blXl
1 + X1
1-x ]I -
+ c* in ~
y(t)'
(3.14)
Lv-Approximation Method
293
with
X, =
2/32
X2 =
2&
,
(3.15(a))
,
(3.15(b))
a o + a~ t + a 2 t 2
A* = ( t -
X,)(t-
X2)'
(a.15(e))
Ot0 + Ot1 X 1 + a 2 X 2 B* =
C* =
(3.15(d))
( x~ - t ) ( x~ - x ~ ) ' a o + a 1 X2 + a 2 X 2 (X2 _ t)(X 2 - X~)"
(3.15(e))
T o d e t e r m i n e t h e coefficients, A, for a n i n t e g r a l w i t h l e n g t h 1 ( c o m p a r e t h e r e m a r k for t h e special case b - a = 1 m e n t i o n e d a f t e r t h e t h e o r e m in S e c t i o n 1), t h e v a r i a b l e t ~ [ - 1, 1] is m a p p e d o n t o u ~ [ - 1, !]2 u s i n g t h e t t r a n s f o r m a t i o n u = 3. T h e L K a p p r o x i m a t i o n p r o b l e m t h e n b e c o m e s 1
r( A, u)] p du
min m
(3.16)
.
__ 2
C a l c u l a t i n g t h e m i n i m u m v a l u e of t h e a b o v e a p p r o x i m a t i o n p r o b l e m , t h e b e s t r e s u l t s in t h e L K n o r m for d i f f e r e n t v a l u e s of p a n d for t h e a b o v e c h o i c e of a p p r o x i m a t i o n f u n c t i o n a r e g i v e n in T a b l e s 2a a n d 2b.
T A B L E 2a CALCULATED COEFFICIENTSOF THE APPROXIMATIONFUNCTION
p=l
p=2
p=4
p=8
p=16
- 0.015901 0.499960 0.009663 0.000199 - 0.952559
0.000132 0.457612 -0.000090 -0.000002 - 0.961707
- 0.000002 0.405397 0.000001 0.000000 - 0.968475
-0.057666 0.374331 0.043482 0.000682 -0.971525
- 0.012619 0.366593 0.009586 0.000144 -0.971796
294
J. H. DE KLERK, D. EYRE, AND L. M. VENTER
TABLE 2b CALCULATED MINIMUM /p-NORMS OF THE RESIDUAL FUNCTION
p=l
p=2
p=4
p=8
p=16
0.057091
0.090140
0.138462
0.174507
0.189342
EXAMPLE 3. T h e last example is a singular integral equation from the theory of q u a n t u m scattering (Reid [14])
m( q, K ) = v( q, K ) - 2 ~ v (
q, ~/ ) m( ¢, K ) - -
at2 _ /£2"
(3.17)
The p a r a m e t e r K is defined by the scattering energy
2 = 27rM > 0,
(3.18)
where M is the particle mass and h is Planck's constant. T h e function v( q, q') is the Fourier transformed potential, while m(q, K) is a real valued function, normalized such t h a t m( K, K) = - [K cot 8( K)] - 1 ,
(3.19)
where 8(K) is the scattering phase shift. In nuclear scales the ratio h2/(2 ~-M) = 41.47 MeV.fm 2. The calculations are carried out at a scattering energy e = 24 MeV. T h e potential function v(q, q') is chosen to have a form t h a t is typical of phenomenological two-nucleon potentials. For this purpose the Fourier transform of the Reid potential [14] is taken,
v( q, ~ )
~M 3 2 h2~l q~/ ~., V~ in i=1
( q + 4) 2 + (q_
] (3.20)
where ~[~1 ---~ 0 . 7 f m - 1 , ~$.L2 = 4#1 , IZ3 = 7/zl, V1 = - 1 0 . 4 6 3 MeV.fm -3, V2 = - 1650.6 MeV.fm -3, and V3 = 6484.2 MeV.fm -3. This last example represents a realistic problem t h a t requires a numerical procedure to evaluate the integrals. A detailed discussion of the approximation functions and the numerical procedures is now given.
Lp-Approximation Method
295
The trial functions are chosen as Chebyshev polynomials of the second kind, U~(t). In the q-variable the trial functions are obtained by mapping t onto the semi-infinite domain 1
q = ~" - ] - - ~ ]
,
(3.21)
where ~ > 0 is some appropriately chosen mapping parameter. The value of remains fixed throughout the calculation. This parameter was chosen to be ~"= 5. In the mapped variable the trial functions become 3
~i(q) = ( 2 ~ ) ~ (q2 + ~2) Ui q2 + ~.2 •
(3.22)
The trial functions satisfy the orthogonality relation
O,(q)~bj(q) q2 + ~2 dq = 6,j,
(3.23)
and behave asymptotically like q-:. This asymptotic behavior of the functions qJ, makes the trial functions particularly well suited for approximating the solution of the quantum scattering integral equation. For this example the approximation function is the linear combination N
m(q,K) = F( A(K),q) = ~_.a,(K)~O,(q).
(3.24)
i=0
The principal value singularity is evaluated using the method of subtraction. Making use of the identity
f
d4
0,
(3.25)
296
J.H. DE KLERK, D. EYRE, AND L. M. VENTER
the residual function becomes
N[
r ( A , q ) = ~ ~O~(q) 4=0
+ --Tr
q~2
K2
× a i ( K ) - v(q, K).
(3.26)
The integral is evaluated numerically using a standard Gauss-Legendre quadrature formula. The points and weights are mapped onto a semi-infinite domain q' ~ [0, ~) yielding the set {q,, w~}L=] where L is the number of quadrature points. An additional point qL+l K is included in this set of points. The residual function is now approximated by the double summation =
2 L+I ~b~(q) + -
~ ( A , q ) = ~] i=O
17"
]
F, v(q, q)~,(Owj ~,(,,1- v(q,K), 3=1 (3.27)
where qj2coj 2 K2, qj--
j= 1,...,L, (3.28)
w~ =
L
K2(.Ok
-Fl= qkff---2'K j=L+I To determine the coefficients, A, the variable q is mapped onto the t The final step is interval [ - 1, ½] using the transformation u = 7.
1 min A
]~( A, u)l p du
,
(3.29)
L p-Approximation Method
297
to obtain the coefficients A and the approximate solution of (3.17), N
q, ,,) = E
(3.30)
i=0
Calculating this m i n i m u m value, the best results in the Lp-norm for different values of p and for the mentioned choice of approximation function are given in Tables 3a, 3b, and 3c. 4.
DISCUSSION OF RESULTS
T h e results given in the previous section represent the " b e s t " results for a given approximating function in the sense t h a t they produce the smallest residual function in the Lp-norm. For a small n u m b e r of trial functions it can be seen t h a t each Lp-approximation produces a different solution. T h e usual p-norms such as p = 1, 2 and ~ will normally suffice for most solutions; however, other values of p, namely 2 < p < ~, can be used to obtain approximations for the L~-solution where this is difficult to obtain. Some general features of the approximation for a given p can be seen in the numerical results. For example, it is noted in Example 1, for the case p = 1, N = 2, t h a t the absolute value of the greatest difference between the analytical solution and the Ll-approximate solution is 0.073169 at the end
TABLE 3a RESULT USING N = 2 TRIALFUNCTIONSSHOWINGTHE CALCULATED PHASESHIFT (~AND Lp-NORMOF THE RESIDUALFUNCTION
1 2 3 4 5 10 15 20 40 80 160 :¢
26.45 31.72 32.78 33.53 33.60 34.54 34.82 34.87 35.10 35.22 35.24 35.35
0.092 O.114 0.125 0.133 0.139 0.153 0.159 0.163 0.170 0.175 0.177 0.180
298
J . H . DE KLERK, D. EYRE, AND L. M. VENTER TABLE 3b CALCULATED PHASE SHIFT ~ USING NTRIAL FUNCTIONS
N
p=l
p=2
p=3
p=4
2 3 4 5 6 7
26.45 39.12 38.58 38.98 39.05 38.84
31.72 38.68 38.73 38.70 38.70 38.84
32.78 38.39 38.75 38.70 38.74 38.85
32.53 38.22 38.77 38.71 38.74 39.02
points. This relatively large error at the boundaries of the interval is typical of approximations in the Ll-norm. Disregarding these extremities, and considering the interval [0.01, 0.99], the greatest difference is only about 0.04. In general the approximating function used in Example 2 can have poles on the real line, or in the complex plane close to the real line, that are not present in the analytical solution. The reason for choosing m - - - 2 is to produce a pole structure in the approximate solution that simulates the end-point singularities in the analytical solution. In all the examples, the values calculated for the Lp-norms of the residual function increase monotonically as p increases. This is consistent with the theory of Lp-spaces; under the conditions imposed in the various examples (the length of the interval of integration is 1) the Lv-norm of a fixed function would increase as p increases. Because the Lp-norm is less than the L~-norm, one can also reasonably expect the norm of the L~-residual (if it exists) to be larger than the norm of the Lp-residual. This is illustrated by Example 3. In fact, from Table 3a it is clear that the norms of the residuals form an increasing sequence that converges to the norm of the L~-residual.
TABLE 3c CALCULATED MINIMUM Lp-NORMS OF THE RESIDUAL FUNCTION
N 2 3 4 5 6 7
p=l 0.09195 0.02141 0.00855 0.00305 0.00285 0.00184
p=2 0.11399 0.02526 0.00969 0.00521 0.00509 0.00268
p=3 0.12542 0.02731 0.01040 0.00638 0.00637 0.00373
p=4 0.13313 0.02867 0.01087 0.00727 0.00723 0.00447
L p-Approximation Method
299
In Example 2 the analytical solution is not in L~, because it is unbounded near the end points of the interval. Even in this case, an approximate solution might exist in the L~-sense. If such a solution exists, the norm of the residual will be larger than the values calculated for large p.
5.
CONCLUSIONS
An important feature of the Lp-approximation method is a good choice of approximating function. The choice of approximating function is guided by general information that can be obtained about the solution of the integral equation such as its singular or asymptotic behavior. Where a good choice of the approximating function is possible, it is found that the Lp-approximation method is an efficient and accurate method for solving the integral equation. This paper has shown using examples of integral equations that exhibit different singular behavior, that solutions of singular integral equations can be obtained from the Lp-approximation method via a technique developed previously. The relationship between Lp-norms allows bounds to be placed on the residual function for different values of p. The Lp-approximation method can be used to construct solutions in the L~-norm even when these solutions are difficult to obtain directly by L~-methods. REFERENCES 1 L.M. Delves and J. L. Mohamed, Computational Methods for Integral Equations, Cambridge University Press, Cambridge, 1985. 2 K.E. Atkinson, A Survey of Numerical Methods for the Solution of Fredholm Integral Equations of the Second Kind, SIAM, Philadelphia, 1978. 3 J.A. Cuminato, On the uniform convergence of a collocation method for a class of singular integral equations, BIT 27:190-202 (1987). 4 J.A. Cuminato, Uniform convergence of a collocation method for the numerical solution of Cauchy-type singular integral equations: a generalization, IMA J. Numer. Anal. 12:31-45 (1992). 5 T. Diogo, S. McKee, and T. Tang, A Hermite-type collocation method for the solution of an integral equation with a certain weakly singular kernel, IMA J. Numer. Anal. 11:595-605 (1991). 6 H. Kaneko and Y. Xu, Numerical solutions for weakly singular Fredholm integral equations of the second kind, Appl. Numer. Math. 7:167-177 (1991). 7 B.V. Riley, The numerical solution of Volterra integral equations with nonsmooth solutions based on sinc approximation, Appl. Numer. Anal. 9:249-257 (1992). 8 G. KSthe, Topological Vector Spaces I, Springer-Verlag, Berlin, 1969.
300
J . H . DE KLERK, D. EYRE, AND L. M. VENTER
9 J.H. De Klerk, Solutions of integral equations via Ll-approximations, Int. J. Control 48:2121-2128 (1988). 10 J.H. De Klerk, The sensitivity of a computational Ll-approximation technique for integral equations, Appl. Math. Comput. 62:171-182 (1994). 11 D.P. Laurie and A. Venter, Automatic Quadrature of Functions of the Form If(x)l, Technical note T/15, Potchefstroom University, 1993. 12 R. Piessens, E. De Doncker-Kapenga, C. W. Uberhuber, and D. K. Kahaner, QUADPACK--A Subroutine Package for Automatic Integration, SpringerVerlag, Berlin, 1983. 13 W.H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling, Numerical Recipes in Pascal: The Art of Scientific Programming, Cambridge University Press, Cambridge, 1989. 14 R. V. Reid, Local phenomenological nucleon-nucleon potentials, Ann. Phys. 50:411-448 (1968).