AMSC/CMSC 460 Project 3 Due 3/27In this project we are going to use numerical integration combined with linear algebrato solve integral equations. The scripts and functions you will be asked to produce areall quite short; but they are tricky to write and debug. You will also learn more aboutpassing functions and their arguments to other functions in Matlab.Fredholm integral equationsA Fredholm integral equation has the form1g(s) =Z10K(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) wewant to determine the function f(t). If µ is zero, (1) is a Fredholm equation of thefirst kind; otherwise it is an equation of the second kind. Equations of the first kindhave many important applications: computer aided tomography and image deblurringamong others. But they are generally hard to solve. Here we will consider equations ofthe s ec ond kind, which are usually more tractable.A strategy for computing an approximationOnly in very simple cases is it poss ible to find an analytic solution of (1). Instead we willcompute an approximation to the solution by replacing the integral with an integrationformula. The process can be illustrated by the by the trapezoidal.We begin by dividing the interval [0, 1] into n equally spaced points, s1= t1, s2=t2, . . . , sn= tn. Setgi= g(si),If we knew the values of f exactly we would havegi=Z10K(si, t)f (t) dt − µf(si). (2)Since we do not, we try to determine approximations fito f (ti), that satisfygi= CT10K(si, t)f (t) − µfi,where CT is the composite trapezoidal rule on the points t1, . . . , tn.Let h = 1/(n−1) be the distance between consecutive ti. Then replacing the integralin (2) by the composite trapezoidal rule, we get the equationgi=h2(K(si, t1)f1+ 2K(si, t2)f2+ · · · + 2K(siti, tn−1)fn−1+ K(si, tn)fn) − µfi.1The notation used here is slightly nonstandard.1AMSC/CMSC 460 Project 3 Due 3/27This equation must hold for i = 1, . . . , n. This gives us n linear equations, which maybe solved for the fi.The matrix form for the system system is easily exhibited. Let kij= K(si, tj). Thenwe have for n = 5g1g2g3g4g5=h2k11− 2µ/h 2k122k132k14k15k212k22− 2µ/h 2k232k24k25k312k322k33− 2µ/h 2k34k35k412k422k432k44− 2µ/h k45k512k522k532k54k55− 2µ/hf1f2f3f4f5(3)Thus our algorithm is to evaluate the gi, form the matrix in (3), and use the Matlaboperator \ to solve the system for the fi.Generating test casesOne of the difficulties in debugging any scientific program is generating test cases withknown solutions. This is particularly so of our problem, since in will in general beimpossible to integrate (2) analytically. Of course, we could use numerical integrationto ge nerate 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 beKp(s, t) =1 − (s − t)2p(4)and we will let ftbe a low order polynomial; e.g., 1 + t2. Since for fixed s, Kp(s, t) is apolynomial of degree 2p in t. Hence if the degree of f is q, then a Gaussian quadraturewhose order ` satisfies 2` − 1 ≥ 2p + 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 GAUSS24% will return a vector of integrals.2AMSC/CMSC 460 Project 3 Due 3/27Note that the specifications are in Matlab format, which requires that variables andfunctions be written in capital letters for emphasis. They should be reduced to lowercase in your code. The ordinates and abscissas for the quadrature will be provided.The argument varargin is one solution to a problem that arises frequently in sci-entific computation. The utility function gauss24 must evaluate a function. Its nameis provided by the argument fun, and gauss24 will generate the values t at which it isto be evaluated. However, fun requires more information — for example, the order p ofthe kernel. One solution would be to include the additional information in the callingsequence of gauss24; e.g.,intgrl = gauss24(0, 1, ’myfunction’, p);But then gauss24 would not work with another function that did not require p orrequired more than one additional argument. Matlab’s solution is to collect all the finalarguments after ’myfunction’ into a cell array and pass it to gauss24. When gauss24invokes feval to evaluate fun, Matlab unpacks varargin and places its objects at theend of the calling sequence of fun. For further details, use the help command.To evaluate the giin (2) you must write two further functions:function y = kerf(t, f, s, p)% KERF(T, F, S, P) is called by gauss24 to evaluate% the integrand of the integral equation (1) (in% the writeup for Project 3) for the scalar T and% the vector S. F is the name of the function in% the integrand and P is the order of the kernel (4).% The calling sequence for F is%% FEVAL(F, T);andfunction g = ggen(s, f, mu, p)% GGEN(S, F, MU, P) generates the vector G defined by (2) in the% writeup for Project 3. S is the vector of S(i), F is the name% of the function f in the kernel with calling sequence%% FEVAL(F, T);%% MU is as in equation (2) and P is the order of the kernel.The body of each function can be written in a single line of code.3AMSC/CMSC 460 Project 3 Due 3/27Setting up the equations.We have now generated a vector g that was evaluated at the points in the vector s forthe function f. We must now generate the matrix in (3), but for ge neral n. Write thefollowing 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