Linear Methods of Applied Mathematics

An example of an integral equation

(c) Copyright 1994,1997,2000 by Evans M. Harrell II and James V. Herod. 
     All rights reserved.

Notes for the instructor.

This contains calculations and examples which correlate with chapter 14 of the WWW text by Harrell and Herod.

Students can be encouraged to cut and paste from this notebook to do homework.

Instructions

This notebook uses Mathematica  to perform calculations for Harrell and Herod's hypertext book, Linear Methods of Applied Mathematics.   The student needs only a basic knowledge of Mathematica  to use the notebook, which is designed both to show how to work problems in the text and to provide a template for using Mathematica  to work other problems of the student's own design.  

Calculations will be performed when the reader presses enter in a given calculation cell (bold print).  It is best to activate the cells in order, so that Mathematica  will be able to call on operators and functions defined in earlier cells.  Red color coding is used to warn the reader when a given calculation relies on an earlier one.

Example of an Equation with a Small Integral Operator

Let's consider the integral equation

         y(x) =  K y  + f(x)
where K g   denotes the integral from 0 to x of Sin[Pi(x-t)] g(t) dt, we take  f(x) = x, and the variable x runs from 0 to 1.

This kernel is not separable, even though a trigonometric identity can be used to write Sin[Pi(x-t)] as a sum in the form ap(x) bp(t), because the limits of the integral depend on x.  (The equation is said to be of Volterra  type.)

We first ask whether this integral kernel is small in the two senses we have learned about.  This means checking whether either of two calculations is less than 1.  Notice that since the integral runs only from 0 to x, the integral kernel is
          K(x,t) = Sin[Pi(x-t)] , 0 < t < x < 1
          K(x,t) = 0 , 0 < x < t < 1
The integrals we shall calculate will thus only run from 0 to x in the variable t.  

To see if the kernel is small in the sense appropriate for a continuous  solution y(x), we calculate:

Integrate[Sin[Pi(x-t)], {t,0,x}]

This is maximized by 2/π, which is less than 1:

N[2/Pi]

We are thus guaranteed that the integral equation has a continuous solution, and that we can calculate it by iteration.

To see if the kernel is small in the sense appropriate for a square-integrable  solution y(x), we calculate:

Integrate[Sin[Pi(x-t)]^2, {t,0,x}]

Integrate[%, {x,0,1}]

N[%]

This is also less than 1, so we are guaranteed that iteration will lead us to the solution in the r.m.s. sense.  Normally this is a weaker statement than the previous one, so we are not surprised.

Next we calculate the solution

Solving the Equation

Clear[K]

K[f_] := Integrate[Sin[Pi(x-t)] (f /. x -> t), {t,0,x}]

Explanation of the syntax

The possibly puzzling command embedded in this, (f /. x -> t), has the effect of changing the variable in the function f from x to t, which we need to do to integrate properly.

y1[x_] = K[x]

y2[x_] = K[y1[x]]

y3[x_] = K[y2[x]]

{Plot[x, {x,0,1}], Plot[y1[x]+x, {x,0,1}],     Plot[y2[x]+y1[x]+x, {x,0,1}],     Plot[y3[x]+y2[x]+y1[x]+x, {x,0,1}]}

Error Analysis

These last two plots look somewhat similar, which is an indication that we are getting close to the exact solution.  But how far away are we really?  Well, if you stop the series at third order,

        f + Kf + K2 f + K3 f,
then the error is   K4 f + K5 f +...

We know that the calculations we must check for convergence are actually bounds 
of the form  ||K g|| <= c || g||,
where for the infinity norm (for continuous functions), c = 2/π, and for the r.m.s. norm (for square-integrable functions) c = 1/2.  (The calculation we did gave 1/4, but that is actually c2).

The right scale for this problem is   ||x||infinity = 1   or, respectively, ||x||2 = 1/Sqrt[3]

N[1/Sqrt[3]]

The norm of Kn f is at most  Kn  times the norm of f (either norm), so the size of the error is
c4 + c5 + ... = c4 (1 + c + c2 + ...) = c4/(1-c).

Conclusion:  The maximum possible error
													 |x+y1+y2+y3   -  yexact(x)|
is at most

N[(2/Pi)^4 /(1-2/Pi)]

The maximum possible error in the rms sense is

(1/2)^4 / (1-1/2)

Which is fairly small.  If we want more accuracy, we need more terms.  With 8 terms,

N[(2/Pi)^8 /(1-2/Pi)]

N[(1/2)^7]

So we would be within 5% uniformly and about 1% in the rms sense.  Actually, we are probably doingmuch better, for already at third order,

Plot[y3[x], {x,0,1/2}]
Plot[y3[x], {x,1/2,1}]

y4[x_] = K[y3[x]]
y5[x_] = K[y4[x]]
y6[x_] = K[y5[x]]
y7[x_] = K[y6[x]]
y8[x_] = K[y7[x]]

What beautiful functions!

You can hardly see the difference between the really accurate solution and the third order solution:

Plot[{x+y1[x]+y2[x]+y3[x]+y4[x]+y5[x]+y6[x]     + y7[x]+y8[x],      x+y1[x]+y2[x]+y3[x]},    {x,0,1}]

Plot[{x+y1[x]+y2[x]+y3[x]+y4[x]+y5[x]+y6[x]     + y7[x]+y8[x],      x+y1[x]+y2[x]+y3[x]},    {x,.99,1}]