AM505a | PDE

Solution of Linear Integral Equations

Robert M. Corless
November 25, 1998

References 2.1 Laplace Transforms

The main references used for this are [2], [6] and [7]. If the integral equation is of convolution type, that
A. Tadesse used [4, 5, 3] for supplementary material is, like Poisson's equation
for this morning's lecture. Zx
See http://www.gams.nist.gov for numerical soft- (x) = K (x ?  )( ) d + f (x) (1)
ware for the solution of integral equations; the pa- 0

per [1] is a pointer to the numerical literature. or Abel's equation

Z x
1 Overview
K (x ?  )( ) d = g(x) ; (2)
0

In this lecture we look at several analytical techniques then the convolution integrals present suggest that
for solving linear one-dimensional integral equations; we take Laplace transforms. Remember
we look at the Maple share library package IntSolve, Z1
written by Honglin Ye [Ph.D. Western]; we look L(f )(s) = e?st f (t) dt (3)
at well-posedness of linear integral equations with 0

smooth kernels; and brie y at approximate and nu- L(f  g) = L(f )L(g) (4)
merical methods. Clearly we're going to have to move
at some speed! where the convolution product f  g is
You have now seen the general Neumann series, Z x
and some examples; you also seen how to solve in- f g = f (x ?  )g( ) d
tegral equations with degenerate kernels. This will Z0 x
be quite useful for us when we look at Marchenko's = f ( )g(x ?  ) d :
equation for soliton problems. 0
We also saw brie y this morning an introduction to Hence, Poisson's integral equation becomes
eigenvalue techniques. We look brie y at some more
examples here.
L() = 1 ?L(Lf()K )
2 Solution Techniques provided that L(K ) =
6 1. Therefore

We now look at some very simple but powerful tech-
 = L?1 L(f )  : (5)
niques for some important classes of problems. 1 ? L(K )
1
For example, if our Poisson's equation is where we want to nd the coecients ak . Substitut-
Z x
ing this into the integral equation gives
(x) = ex?y (y) dy + sin(x) ; (6) X
( ak ?  k ak k ? b k k ) = 0 :
0 k
k
then since L(exp(x)) = 1=(s ? 1) and L(sin(x)) =
1=(s2 + 1), Thus, choosing
! ak = 1 ?bk (12)
(x) = L?1 1 k
(s2 + 1)(1 ? s?1 1 ) solves the problem, unless one of the k = 1.
 For example, consider the Fredholm rst kind in-
?1 
= L?1 (s2 +s1)( tegral equation (the above analysis was for a second-
s ? 2) kind integral equation but it all goes through mutatis
= 5 e ? 5 cos x + 53 sin x :
1 2x 1
(7) mutandis):
Substituting this equation back into the original 1 Z 1? 2
equation shows that this really is a solution. 2 ? 1 ? 2 cos(x ? y) + 2 (y) dy = f (x) :
(13)
We take 0 < < 1 and ?  x  . The eigenfunc-
2.2 Eigenfunction expansions tions are just cos(kx) and sin(kx). The solution can
Suppose we are trying to solve be shown to be
X ?n
Z (y) = 21 a0 + (an cos ny + bn sin ny) ; (14)
(x) = K (x;  )( ) d + f (x) ; (8) n1

where the a and b are the Fourier coecients of f .

and we know a set of eigenfunctions k (x) such that One oftenksees thek eigenfunction expansion method
Z applied to degenerate kernels. For example,
K (x;  ) k ( ) d = k k (x) ; (9) Z 1 
(y) = x y
2 2
1 ? 2 (x) dx + sin y : (15)
for appropriate eigenvalues k . We will concern our- 0

selves here only with the case where the set of eigen- There are only two eigenfunctions. The eigenvalues
values is nite or at most countably in nite, but in are  = 9=20  p889=60, with associated eigenfunc-
the next two weeks we will see examples of integral tions Ai + Bi x2 . The solution is
equations with a continuous spectrum of eigenvalues.
[This is of some importance for physics, and some of 3 61 + 20
2
(y ) = sin y +  3 ? 2 y : (16)
6 2
you will have seen this already.] 5

Here, when we expand everything in sight, Of course, I did that with Maple, not by hand!
X
f (x) = bk k (x) (10)
k 3 Maple code for solution of
where the bk are presumed known because f (x) and linear integral equations
the eigenfunctions k (x) are known, and
X
In [7] we nd a description of IntSolve, which resides
(x) = ak k (x) (11) in the share library of Maple. Watch out for bugs.
k The code has \rotted", to use a technical term|that

2
is, Maple has evolved since the package was written, > p2 := IntSolve( example2, phi(y),
eigenfunc );
and it is not clear that the package still works. In
particular I have found some bugs in the eigenfunc p2 := sin( y) +
3 61 2 + 20 ? 6 y2
method. 5 3 
> P2 := unapply(p2,y);

P2 := y ! sin( y) + 3 61  3+ 20 ? 6 y
2 2
4 Solution of Linear Integral
Equations in Maple
5  
> with(share); > value( eval( subs( phi=P2, example2 )
) );
See ?share and ?share,contents for

sin( y) + 35 61 3+ 20 ? 6 y =
information about the share library 2 2

[]
? 53 ?61  + 10 y2 2 ? 20 + sin( y)
2

3
> with(IntSolve);
> normal( lhs(%) - rhs(%) );
Share Library: IntSolve
0
Authors: Ye, Honglin and Corless, Robert. > example3 := phi(x) = lambda*Int(
exp(k*(x-y))*phi(y), y=0..1) + sin(x);
Description: Version 1 of
Z 1
an integral equation solver. example3 := (x) =  e(k (x?y)) (y) dy + sin(x)
0

[\IntSolve"] Two terms of the Neumann series:

> IntSolve( example3, phi(x), 2 );
> ?IntSolve
> example1 := phi(x) = Int(
(sin(x) k2 + sin(x) ?  e(k (x?1)) k sin(1)
exp(x-y)*phi(y), y=0..x) + sin(x);
Z
? e(k (x?1)) 2 k sin(1) ? e(k (x?1)) 2 cos(1)
example1 := (x) =
x
e(x?y) (y) dy + sin(x) ?  e(k (x?1)) cos(1) + 2 e(k x) +  e(k x) )=(
0 k2 + 1)
> p1 := IntSolve( example1, phi(x), Can try to convert to a di erential equation:
Laplace ); > p3 := IntSolve( example3, phi(x),

p1 := 1 e(2 x) ? 1 cos(x) + 3 sin(x)

differentiate);
5 5 5 (
p3 := D()(x) ? k (x) ? cos(x) + k sin(x) = 0;
> example2 := phi(y) = Int(
(1-x^2*y^2/2)*phi(x), x=0..1) + Z 1
)
sin(Pi*y);
(0) ?  e(?k y) (y) dy = 0
0
example2 :=
Z
(1 ? 12 x2 y2 ) (x) dx + sin( y)
1
>
(y) = P3sol := dsolve(p3[1], phi(x));
0 P3sol := (x) = sin(x) + e(k x) C1

3
 > P3 := unapply( rhs(P3sol), x ); Then we have found the exact solution to the related
P3 := x ! sin(x) + e (k x)
C1 integral equation
Z
> eval(subs(phi=P3,p3[2])); (x) = K (x;  ) ( ) d + f^(x) ; (18)
C1 ? (?cos(1) e(?k) ? k e(?k) sin(1) + C1 k2 where f^(x) = f (x) + r(x). This is trivial. But it's
+ C1 + 1)=(k2 + 1) = 0 also profound. If r(x) is small compared to the ap-
proximations already made in deriving the integral
> map(normal, readlib(isolate)(%,_C1) equation, or to physically reasonable perturbations
); in f , then we are done.
(?k) (?k)
C1 =  (cos(1) e 2 + k e 2 sin(1) ? 1)
Perhaps instead we have exactly computed a so-
?k ? 1 +  k +  lution satisfying an equation with a modi ed kernel
That's very impressive, but is it right? (one technique is to replace K (x;  ) with a degener-
ate kernel K^ that approximates K , for example). Is
> assign(%); K^ ? K smaller than terms neglected in the model?
> value( eval( subs(phi=P3, example3 ) Yes? Done.
) ); This raises the question of how sensitive the solu-
tion of the integral equation is to changes in f or K .
sin(x) + e  (cos(1)
(k x)
e(?k) + k e(?k) sin(1) ? 1) = In general this is a deep subject (with singular ker-
?k2 ? 1 +  k2 +  nels, for example). But for our case, we nd in [2,
(k e(k x?k) sin(1) ?  k sin(1) e(k x?k) p. 156] that there are computable constants N  and
C such that
+  k e(?k) sin(1) e(k x) 
+  cos(1) e(?k) e(k x) + cos(1) e(k x?k) j^ ? j  ("N1 ?+"(1)(1
+
+ C)
C ) (19)
?  cos(1) e(k x?k) ? e(k x) )=((k2 + 1) (?1 + )
) + sin(x) where " is a bound on the distance jjK^ ? K jj and 
is a bound on jjf^ ? f jj. This is as good a behaviour
> simplify( normal( lhs(%) - rhs(%) ) as can be expected, and shows that not only is the
); problem well posed, we have a linear dependence on
0 the perturbations.

6 Collocation
5 Well-posedness
This is a reasonable numerical method to solve many
There are many methods to approximate the solution kinds of integral equations, not just the simple ones
of an integral equation. There are many methods to here. See however the many numerical papers, start-
solve them numerically. But before we show how to ing perhaps from those of Hermann Brunner (Memo-
solve problems approximately, we observe that lin- rial U.), e.g. [1], or the software available through
ear integral equations with smooth kernels are well- http://www.gams.nist.gov, for real methods. The
posed. Suppose for example we have somehow com- following is just a sketch.
puted an approximate solution (x) to the Fredholm Choose points (\collocation points") xk on a  x 
2nd kind equation (8). De ne the residual b. Then
Z Z b
r(x) := (x) ? K (x;  )( ) d ? f (x) : (17) (x) = K (x;  )( ) d + f (x)
a

4
implies > phi(x) = Int( exp(x*y)*phi(y),
Z b y=0..1) + 1/(1+x^2);
(xk ) = K (xk ;  )( ) d + f (xk ) : (20) Z
e(x y) (y) dy + 1 +1 x2
1
a (x) =
Now if we approximate (y) by a linear combina- 0

tion of some chosen \gauge functions" j (y) (say B- > um := %:

splines, or piecewise constants, or whatever), namely >
X IntSolve( um, phi(x),2 );
(y) = j j (y ) ; (21) Warning, computation interrupted
j
then substitution into (20) gives a set of linear equa- > restart;
tions for the unknown
Rb
coecients j , which depend
on the integrals a K (x;  ) j ( ) d . We choose the > p := proc(N::posint) local A, b;
guage functions to make these simple to evaluate, of > A := matrix(N,N,
course. > (i,j)->evalf((exp((i-1/2)*j/N^2)
> -exp((i-1/2)*(j-1)/N^2))*N/(i-1/2)));

6.1 Example >

>
b := vector(N, i->evalf(1/(1+(i-1/2)^2/N^2)) );
linalg[linsolve](evalm(-A+&*()),b)
Suppose we wish to solve > end:
Z
exy (y) dy + 1 +1 x2
1 > p(2);
(x) = (22)
0 [?1:595537339; ?2:795225413]
Choose xk = (k ? 1=2)=N for 1  k  N , equally
spaced on 0  x  1. Choose our gauge functions to > p(4);
be piecewise constant, j (x) = j on (k ? 1)=N  [?1:248020298; ?1:721998520; ?2:319643386;
x  k=N and zero otherwise. Then our collocation
equations become ?3:002682256]
XN Z k=N
ey(j?1=2)=N dy +  1 2 ;
> p(10);
j= k
k=1 (k?1)=N
1 + j?N1=2 [?1:103197136; ?1:253032875; ?1:430046122;
(23) ?1:631232552; ?1:852947221;
and rearranging this gives a simple linear system
8 9
?2:091770186; ?2:345087651;
>
< >
=
?2:611336610; ?2:889992145;
~
(A ? I ) = >  ?1 ?3:181416068]
2 >
: 1 + j?1=2 ;
N > p(30);
where the matrix entries are Aij =
[?1:053813402; ?1:097363590; ?1:144113251;
N exp  (i ? 1=2)j  ? exp  (i ? 1=2)(j ? 1)  : ?1:194031617; ?1:247062256;
i ? 1=2 N2 N2 ?1:303125609; ?1:362122403;
(24) ?1:423937411; ?1:488443464;
We hope that as N ! 1, the j ! (xj ). As the ?1:555505672; ?1:624985325;
following Maple session shows, it works. ?1:696743735; ?1:770645567;
> restart; ?1:846561767; ?1:924371993;
?2:003966513; ?2:085247612;
5
 ?2:168130508; ?2:252543882; ?1:924371993; x < 158 ; ?2:003966513;
?2:338430027; ?2:425744663;
?2:514456593; ?2:604547123; x < 17 ; ?2:085247612; x < 3 ; ?2:168130508;
?2:696009313; ?2:788847278; 30 5
?2:883075290; ?2:978716985; x < 19 ; ?2:252543882; x < 2 ; ?2:338430027;
?3:075804529; ?3:174377822; 30 3
?3:274483757] 7 ; ?2:425744663; x < 11 ;
x < 10 15
23 ; ?2:604547123; x < 4 ;
?2:514456593; x < 30
> p30 := %: 5
> plot([seq( [(k-1/2)/30,
p30[k]],k=1..30 )] );
?2:696009313; x < 56 ; ?2:788847278; x < 15
13 ;

?2:883075290; x < 109 ; ?2:978716985;

x < 14
15 ; ?3:075804529; x < 29 ;
30
–1.5 ?3:174377822; x < 1; ?3:274483757)
–2
> term := unapply( int(
exp(x*y),y=(k-1)/30..k/30), x,k);
?e ?1))
term := (x; k) ! e
(1=30 x k) (1=30 x (k
–2.5
x
–3 > res := x -> evalf( ph(x) - add(
p30[k]*term(x,k), k=1..30) - 1/(1+x^2)):
>
0 0.2 0.4 0.6 0.8 1
res(0.1);

> ph := unapply( piecewise( seq(

?:0253523379
op([x<k/30,p30[k]]), k=1..30 )), x);
> plot(res, 0..1, -0.05..0.05);

ph := x ! piecewise(x < 1 ; ?1:053813402; x < 1 ;

30 15
1
?1:097363590; x < 10 ; ?1:144113251; 0.04

2 ; ?1:194031617; x < 1 ; ?1:247062256;

x < 15 6 0.02

1 7
x < 5 ; ?1:303125609; x < 30 ; ?1:362122403;
4 ; ?1:423937411; x < 3 ;
0 0.2 0.4 0.6 0.8 1

x < 15 10
1 11 ;
–0.02

?1:488443464; x < 3 ; ?1:555505672; x < 30

?1:624985325; x < 52 ; ?1:696743735; x < 30
13 ; –0.04

?1:770645567; x < 157 ; ?1:846561767; x < 21 ;

6
References

[1] J. Blom and H.Brunner, ACM TOMS 17 (1991)

pp. 167-177.
[2] Ll. G. Chambers, Integral Equations: A Short
Course, International Textbook Company Ltd.,
1976.
[3] M. A. McKiernan and J. Wainwright, Advanced
analytical methods in applied mathematics, lec-
ture notes for AM964a, 1981.
[4] A. C. Pipkin, A course on integral equations,
Texts in Applied Mathematics, Springer-Verlag,
1991.
[5] W. Pogorzelski, Integral equations and their ap-
plications, vol. 1, Pergamon, Oxford, 1966.
[6] S. L. Sobolev, Partial Di erential Equations of
Mathematical Physics, Dover, 1964.
[7] Honglin Ye and Robert M. Corless, \Solving lin-
ear integral equations with Maple", Proc. ISSAC
Berkeley, CA, July 27{29 (1992) 95{103