Download Numerical Solution of Fredholm Integral Equations using Trapezoidal and Simpson's Rules - and more Study Guides, Projects, Research Computer Science in PDF only on Docsity!
In this project we are going to use numerical integration combined with linear algebra to solve integral equations. The scripts and functions you will be asked to produce are all quite short; but they are tricky to write and debug. You will also learn more about passing functions and their arguments to other functions in Matlab.
Fredholm integral equations
A Fredholm integral equation has the form^1
g(s) =
0
K(s, t)f (t) dt − μf (s). (1)
The function K(k, t) is called the kernel of the equation. Given a functions g(s) we want to determine the function f (t). If μ is zero, (1) is a Fredholm equation of the first kind; otherwise it is an equation of the second kind. Equations of the first kind have many important applications: computer aided tomography and image deblurring among others. But they are generally hard to solve. Here we will consider equations of the second kind, which are usually more tractable.
A strategy for computing an approximation
Only in very simple cases is it possible to find an analytic solution of (1). Instead we will compute an approximation to the solution by replacing the integral with an integration formula. The process can be illustrated by the by the trapezoidal. We begin by dividing the interval [0, 1] into n equally spaced points, s 1 = t 1 , s 2 = t 2 ,... , sn = tn. Set gi = g(si),
If we knew the values of f exactly we would have
gi =
0
K(si, t)f (t) dt − μf (si). (2)
Since we do not, we try to determine approximations fi to f (ti), that satisfy
gi = CT^10 K(si, t)f (t) − μfi,
where CT is the composite trapezoidal rule on the points t 1 ,... , tn. Let h = 1/(n−1) be the distance between consecutive ti. Then replacing the integral in (2) by the composite trapezoidal rule, we get the equation
gi =
h 2 (K(si, t 1 )f 1 + 2K(si, t 2 )f 2 + · · · + 2K(siti, tn− 1 )fn− 1 + K(si, tn)fn) − μfi. (^1) The notation used here is slightly nonstandard.
This equation must hold for i = 1,... , n. This gives us n linear equations, which may be solved for the fi. The matrix form for the system system is easily exhibited. Let kij = K(si, tj ). Then we have for n = 5
g 1 g 2 g 3 g 4 g 5
h 2
k 11 − 2 μ/h 2 k 12 2 k 13 2 k 14 k 15 k 21 2 k 22 − 2 μ/h 2 k 23 2 k 24 k 25 k 31 2 k 32 2 k 33 − 2 μ/h 2 k 34 k 35 k 41 2 k 42 2 k 43 2 k 44 − 2 μ/h k 45 k 51 2 k 52 2 k 53 2 k 54 k 55 − 2 μ/h
f 1 f 2 f 3 f 4 f 5
Thus our algorithm is to evaluate the gi, form the matrix in (3), and use the Matlab operator \ to solve the system for the fi.
Generating test cases
One of the difficulties in debugging any scientific program is generating test cases with known solutions. This is particularly so of our problem, since in will in general be impossible to integrate (2) analytically. Of course, we could use numerical integration to generate the gi, but that would create errors of its own. We are going to use numerical integration, but in such a way that it is exact. Specif- ically, our kernel will be Kp(s, t) =
1 − (s − t)^2
)p (4)
and we will let ft be a low order polynomial; e.g., 1 + t^2. Since for fixed s, Kp(s, t) is a polynomial of degree 2p in t. Hence if the degree of f is q, then a Gaussian quadrature whose order satisfies 2 − 1 ≥ 2 p + q will integrate (2) exactly. The first step is to code the following function.
function intgrl = gauss24(a, b, fun, varargin) % GAUSS24 Gauss--Legendre integration with 24 points. % INTGRL = GAUSS24(A, B, FUN, VARARGIN) returns a 24 point % Gauss-Legendre approximation to the integral over the % interval [A,B] of the function FUN. FUN is a string % containing the name of the function to be evaluated. It’s % calling sequence (within GAUSS24) is % % VAL = FEVAL(FUN, T, VARARGIN(:)); % % FUN may return a vector of values, in which case GAUSS % will return a vector of integrals.
Setting up the equations.
We have now generated a vector g that was evaluated at the points in the vector s for the function f. We must now generate the matrix in (3), but for general n. Write the following function.
function K = Ktrap(n, p, mu) % K = KTRAP(N, P, MU) sets up the matrix of the system % (3) in the Project 3 writup. P is the order of the kernel, % MU is as in (3).
In our experiments, we want to compare the trapezoidal rule with Simpson’s rule. Code the following function.
function K = Ksimp(n, p, mu) % K = KSIMP(N, P, N) is the analogue of KTRAP but with % the compound Simpson’s rule replacing the trapezoidal rule. % An error return is made if N is even.
You will find your life easier if you first write out the matrix as in (2).
The bottom line
In writing your code you should try to avoid for loops. You will be much helped by the Matlab. operator that implements elementwise operations. For example, the loop
for i=1:n f(i) = 1 + 2*x(i)^2; end
can be replaced by
f = 1 + 2*x.^2;
The fewer loops in your code, the faster it will run and the better your grade will be. Run your routines with the script in Figure 1. It may help you in debugging to know that the plot should be two straight lines with negative slopes. Turn listings of your functions, the plot, and the answer to the following question. What are numbers mysterynumtrap and mysterynumsimp and how do they relate to the errors in the composite trapezoidal and Simpson’s rules?
k = 20; mu = -1;
% Compute the trapezoidal approximations and plot the error.
nrm = []; for n = 21:10: s = linspace(0, 1, n)’; f = fun(s); g = ggen(s, ’fun’, mu, k); K = Ktrap(n, k, mu); ft = K\g; nrm = [nrm, norm(ft - f)/norm(f)]; end loglog(21:10:101,nrm, ’-’) hold nrm [m,k] = size(nrm); mysterynumtrap = log10(nrm(k)/nrm(1))/log10(101/21)
% Compute the Simpson approximations and plot the error.
nrm = []; for n = 21:10: s = linspace(0, 1, n)’; f = fun(s); g = ggen(s, ’fun’, mu, k); K = Ksimp(n, k, mu); ft = K\g; nrm = [nrm, norm(ft - f)/norm(f)]; end loglog(21:10:101,nrm,’.-’) nrm [m,k] = size(nrm); mysterynumsimp = log10(nrm(k)/nrm(1))/log10(101/21) hold
Figure 1: Script for Project 3
