
	Trending ▼ ICSE CBSE 10th ISC CBSE 12th CTET GATE UGC NET Vestibulares ResFinder

integral solution

19 pages, 0 questions, 0 questions with responses, 0 total responses,

	gangu	+Fave Message
	Profile Timeline Uploads

Home > gangu >

Formatting page ...

D. Keffer, ChE 505 ,University of Tennessee, August, 1999 Advanced Analytical Techniques for the Solution of Single- and Multi-Dimensional Integral Equations David Keffer Department of Chemical Engineering University of Tennessee, Knoxville August 1999 Table of Contents 1. 2. 3. 4. 5. 6. 7. Definition of integral equations Ordinary versus partial integral equations Linearity versus nonlinearity of integral equations Single integral equations vs. systems of integral equations Categorization of integral equations Kernels of integral equations Analytical solutions to integral equations Example 1. Volterra equations of the second kind (generalized solution) Example 2. Volterra equations of the first kind (generalized solution) Example 3. Fredholm equations of the second kind (generalized solution) Example 4. Volterra equations of the second kind (specific example) 8. Numerical solutions to linear integral equations 9. Numerical solutions to nonlinear integral equations 10. Numerical solutions to systems of integral equations 11. Numerical solutions to higher-order linear integral equations 12. Numerical solutions to higher-order nonlinear integral equations 13. Numerical solutions to multivariate integral equations 14. Applications of integral equations Example #1: The Ornstein-Zernike closure Example #2: The Yvon-Born-Green Hierarchy 15. Sources for more information on integral equations i 1 1 1 2 2 3 3 3 6 6 7 10 14 18 19 23 25 26 27 29 31 D. Keffer, ChE 505 ,University of Tennessee, August, 1999 Introduction This hand-out is a primer on integral equations. On the analytical side, it assumes that the reader knows calculus, differential equations, and the basic operations of linear algebra. On the numerical side, it assumes that the reader knows the trapezoidal rule for numerical integration and the Newton-Raphson method of root-finding in nonlinear algebraic equations. 1. Definition of Integral equations What is an integral equation? An integral equation (IE) is an equation in which an unknown function appears within an integral, just as a differential equation is an equation in which an unknown function appears within a derivative. Just as the solution to a differential equation is a function, so too is the solution to an integral equation a function. 2. Ordinary versus partial integral equations When dealing with differential equations, we encounter ordinary differential equations (ODEs) and partial differential equations (PDEs). Are there analogous ordinary integral equations (OIEs) and partial integral equations (PIEs)? Yes. There are integral equations in which the integration is carried out with respect to a single variable (OIEs) and there are integral equations in which the integration is carried out with respect to multiple variables (PIEs). However, in the mathematics of integral equations, they are not referred to as OIEs and PIEs. The jargon is different. However, for our purposes, we might as well call them OIEs and PIEs, since we understand that terminology. An example of an ordinary integral equation: x ( x ) = f ( x ) + N (x , y) ( y)dy (1) a An example of a partial integral equation: x 2 x1 (x 1 , x 2 ) = f ( x1 , x 2 ) + N ( x1 , x 2 , y1 , y 2 ) ( y1 , y 2 )dy1dy 2 (2) a 2 a1 Hopefully you see the analogy between ODEs and PDEs with OIEs and PIEs. 3. Linearity versus nonlinearity of integral equations Do integral equations come in linear and nonlinear flavors? Yes. An integral equation is called linear if linear operations are performed in it upon the unknown function, that is, if it has the form x A( x ) ( x ) + B(x ) + N ( x, y) ( y)dy = 0 (3) a 1 D. Keffer, ChE 505 ,University of Tennessee, August, 1999 Here is an example of a nonlinear integral equation x (x ) = f ( x ) + F [ x, y, ( y)]dy (4) a where F [ x , y, ( y)] is a nonlinear function of ( y) , for example F [ x , y, ( y)] = sin[ ( y)] (5) 4. Single integral equations vs. systems of integral equations Do integral equations appear in systems as well as individually. You bet they do. All the time. 5. Categorization of integral equations We have seen that when dealing with algebraic equations, we used only the linear/nonlinear categorization. When dealing with ODEs, we could expand the categorization to include first order/second order, homogeneous/nonhomogeneous, separable/nonseparable, etc. When dealing with PDEs, we categorized them as elliptic, hyperbolic, and parabolic. How do mathematicians categorize integral equations? Mathematicians categorize ordinary integral equations in much the same way that ODEs are classified by small sub-groups, for which a particular solution method exists. Thus, the categorization does not apply to all ODEs (the way that the linear/nonlinear categorization does). Linear OIEs are divided into two basic types: Volterra integral equations, in which the upper limit of integration is variable, x ( x ) = f ( x ) + N (x , y) ( y)dy (6) a and Fredholm integral equations, in which the limits of integration are fixed, b (x ) = f (x ) + N (x , y) ( y)dy (7) a Both Volterra and Fredholm integral equations can be subdivided into two groups: equations of the first and second kind. In equations (6) and (7), we have written Volterra and Fredholm integral equations of the second kind. Volterra and Fredholm integral equations of the first kind have the form, respectively, x f ( x ) = N( x , y) ( y)dy (8) a and Fredholm integral equations, in which the limits of integration are fixed, 2 D. Keffer, ChE 505 ,University of Tennessee, August, 1999 b f ( x ) = N( x , y) ( y)dy (9) a These categorizations apply to both linear and nonlinear OIEs. 6. Kernels of integral equations What is a kernel? In equations (6) to (9), the function N ( x, y) is called the kernel of the integral equation. Every integral equation has a kernel. Kernels are important because they are at the heart of the solution to integral equations. 7. Analytical solutions to integral equations Example 1. We cannot possibly present analytical solutions to al of the types of integral equations. However we can present a few examples. Let us consider the Volterra equation of the second kind. x ( x ) = f ( x ) + N (x , y) ( y)dy (6) a Certain assumptions have to be made about N ( x, y) being defined and bounded over the domain, as well as f ( x ) being a Riemann-integrable function in the domain. Also, ( x ) is an arbitrary Riemann-integrable function in the domain. The first step in the solution of this integral equation is a transformation called the Dirichlet transformation. In the Dirichlet transformation, we state that the order of integration doesn t matter. b q b M ( x , y)dx dy = M ( x, y)dy dx a p a p q (10) Let us define a function M x ( x, y) = N (x , y) N( y, s) (s) (11) Substituting equation (11) into (10), we have xx x M x (s, y)ds dy = M x (s, y)dy ds a a a a (12) xx y N( x , y) N( y, s) (s )ds dy = N( x, y) N( y, s )dy (s)ds a a a s (13) x x We can now proceed to the solution of the Volterra equation. If there exists an integrable function ( x ) which satisfies the Volterra equation (6), then this function also satisfies the so-called iterated equation, obtained by substituting the right hand side of equation (6) under the integral sign of equation (6) obtaining. 3 D. Keffer, ChE 505 ,University of Tennessee, August, 1999 y (x ) = f (x ) + N ( x, y) f ( x ) + N ( y, s) (s )ds dy a a x (14) For the moment, let us proceed under the assumption that this substitution is valid. When we have completed the transformation, we will verify the assumption. We now apply the Dirichlet Transformation to equation (14) obtaining the iterated equation x x (x ) = f ( x ) + N( x , y)f ( y)dy + 2 N1 ( x, y) ( y)dy a (15) a where N1 ( x , y) is called the iterated kernel and is written as x N1 ( x, y) = N( x, s) N(s, y)ds (16) y Repeating the above transformation, we obtain a two-fold iterated equation x x (x ) = f (x ) + [N( x, y) + N1 (x , y)]f ( y)dy + N 2 (x , y) ( y)dy 3 a (17) a where N 2 ( x , y) is called the two-fold iterated kernel and is written as x N 2 ( x , y) = N (x , s) N1 (s, y)ds (18) y After the nth repetition of the transformation we have: n 1 (x ) = f ( x ) + N( x , y) + i N i (x , y) f ( y)dy + n +1 N n ( x, y) ( y)dy i =1 a a x x (19) where N n ( x, y) is called the n-fold iterated kernel and is written as x N n (x , y) = N( x , s ) N n 1 (s, y)ds (20) y Intricate proofs have been performed which show that if N ( x, y) was defined and bounded over the domain all of the i-fold iterated kernels are also defined and bounded over the domain of integration [Pogorzelski]. The standard procedure is to rewrite equation (19) as x (x ) = f (x ) + ( x , y, )f ( y)dy (21) a 4 D. Keffer, ChE 505 ,University of Tennessee, August, 1999 where ( x, y, ) is called the resolvent kernel and is written as ( x , y, ) = N (x , y) + i N i ( x, y) (20) i =1 Now this looks like we cheated in the substitution because, with this resolvent kernel, equation (21) ought to be x (x ) = f (x ) + ( x , y, )f ( y)dy + n +1 a x N ( x , y)[ ( y) f ( y)]dy (22) a but we have omitted the last term entirely. Therefore, we must show that it is important, which, fortunately, mathematicians have already done [Pogorzelski, p. 11]. To check that our solution in equation (21) satisfies the Volterra equation (6), substitute (21) into (6). We need to do this because we made an assumption moving from equation (13) to equation (14) that has not yet been proved. y (x ) = f (x ) + N ( x, y) f ( y) + ( y, s, )f (s )ds dy a a (23) y (x ) = f (x ) + N ( x, y)f ( y)dy + N (x , y) ( y, s, )f (s)ds dy a a a (24) x x x 2 We change the order of integration for y and s x (x ) = f (x ) + N( x, y)f ( y)dy + N( x, y) ( y, s, )dy f (s )ds a a s x x 2 (25) We transpose the dummy variables y and s x (x ) = f (x ) + N( x, y)f ( y)dy + N( x, s) (s, y, )ds f ( y)dy a a y (26) x x (x ) = f (x ) + N( x , y) + N (x , s) (s, y, )ds f ( y)dy a y (27) x x 2 The term in brackets in equation (27) can be rewritten, using the definition of the resolvent kernel in equation (20) x N (x , y) + N ( x, s) (s, y, )ds = (x , y, ) y So that we can rewrite equation (27) as 5 (28) D. Keffer, ChE 505 ,University of Tennessee, August, 1999 x (x ) = f ( x ) + ( x, y, )f ( y)dy (29) a which is what we have in equation (21), which we substituted into the integral equation for checking. So, our check is complete. A Volterra equation of the second kind has one and only one bounded solution, given by the formula in equation (29). This solution is convergent for all values of . This solution technique is also valid for partial Volterra equations of the second kind, x 2 x1 (x 1 , x 2 ) = f ( x1 , x 2 ) + N ( x1 , x 2 , y1 , y 2 ) ( y1 , y 2 )dy1dy 2 (2) a 2 a1 Example 2. Let us consider the Volterra equation of the first kind x f ( x ) = N( x , y) ( y)dy (8) a Using Leibnitz s rule for the differentiation of integrals: ( ) ( ) 2 d2 dF( x, ) d d F( x, )dx = ) d dx + F( 2 , ) d 2 F( 1 , ) d 1 d 1 ( ) 1( (30) we can differentiate both sides of equation (8), effectively transforming the Volterra equation of the first kind into a Volterra equation of the second kind: f ( x ) + N ( x , y) ( y)dy N (x , x ) N(x, x ) a x (x ) = (31) which is a Volterra equation of the second kind. If N(x, x ) = 0 , we can differentiate again to obtain a transformation with N (x , x ) in the denominator. Example 3. Fredholm s integral equation of the second kind. Fredholm s integral equation of the second kind has the form b (x ) = f (x ) + N (x , y) ( y)dy (7) a The solution to this equation is b (x ) = f (x ) + ( x , y, )f ( y)dy (32) a 6 D. Keffer, ChE 505 ,University of Tennessee, August, 1999 where ( x, y, ) is called the resolvent kernel and is written as ( x , y, ) = i 1 N i ( x, y) (33) i =1 where N1 ( x, y) = N( x , y) . Example 4. Let s work a Volterra equation of the second kind with numbers. x (x ) = e x e x y ( y)dy (34) 0 so f ( x ) = e , = 1 , N ( x , y) = e x x y , and a = 0 . x x x y y y N1 ( x, y) = N( x , s) N (s, y)ds = e x s e s y ds = e x y ds = e x y (x y ) x x N 2 ( x , y) = N (x , s) N1 (s, y)ds = e e x x s s y y N 2 ( x , y) = e y x y (s y )ds = e x y (s y )ds y x y y2 2 x y x yx + y = e yx + 2 2 2 2 2 2 2 x x x s2 s2 y2 y2 N 3 ( x , y) = N ( x, s) N 2 (s, y)ds = e x s e s y ys + ds = e x y ys + ds 2 2 2 2 y y y x 3 yx 2 y 2 x y 3 yy 2 y 2 y x 3 yx 2 y 2 x y 3 = e x y N 3 ( x , y) = e x y + + + 2 2 6 2 2 2 2 6 6 6 Let s truncate the resolvent kernel at the third iterated kernel because things are getting ugly fast. ( x , y, ) = N (x , y) + i N i ( x, y) (20) i =1 x2 x 3 yx 2 y 2 x y 3 y2 ( x, y, ) = e x y e x y (x y ) + e x y yx + e x y + 2 2 2 2 6 6 x2 y 2 x 3 yx 2 y 2 x y 3 ( x, y, ) = e x y 1 x + y + yx + + + 2 6 2 2 2 6 Now, let s create the solution using both the first, second, and third approximations to the resolvent kernel 7 D. Keffer, ChE 505 ,University of Tennessee, August, 1999 x y2 1 (x ) = e e [1 x + y] e dy = e e [1 x + y] dy = e e y xy + 2 0 0 0 x x2 02 x2 x 1 (x ) = e 1 x xx + + 0 x 0 + = e x 1 x + 2 2 2 0 x x x 2 (x) = e e x x x y x y 0 y x x x x x y2 x2 y2 y x2 x yx + e dy = e 1 1 x + y + yx + dy 1 x + y + 2 2 2 2 0 x x2 x3 y 2 x y3 y2 x 2 2 ( x ) = e 1 y xy + + + = e x 1 x + y 2 6 2 6 0 2 2 x x x2 y 2 x 3 yx 2 y 2 x y 3 y 3 ( x ) = e x e x y 1 x + y + yx + + + e dy 2 2 6 2 2 6 0 x2 x3 x4 3 ( x ) = e x 1 x + + 2 6 24 We can pretty much guess what the infinite solution looks like: ( x )i i =0 i! ( x ) = e x 1.4 n=2 1.2 n>3 0.8 n=3 phi 1 0.6 0.4 0.2 0 0 0.1 0.2 0.3 0.4 0.5 x 0.6 0.7 0.8 0.9 1 n=1 We can plot the analytical solution, varying the number of iterated kernels that we include in the resolvent kernel. The analytical solution is phi(x) = 1.0. On the following page is the MATLAB code used to generate this plot. 8 D. Keffer, ChE 505 ,University of Tennessee, August, 1999 % % plot % if (j==1) plot (x,phi(:,j),'k-'), xlabel( elseif (j==2) plot (x,phi(:,j),'r-'), xlabel( elseif (j==3) plot (x,phi(:,j),'b-'), xlabel( elseif (j==4) plot (x,phi(:,j),'g-'), xlabel( elseif (j==5) plot (x,phi(:,j),'m-'), xlabel( elseif (j==6) plot (x,phi(:,j),'k:'), xlabel( elseif (j==7) plot (x,phi(:,j),'r:'), xlabel( elseif (j==8) plot (x,phi(:,j),'b:'), xlabel( elseif (j==9) plot (x,phi(:,j),'g:'), xlabel( elseif (j==10) plot (x,phi(:,j),'m:'), xlabel( else plot (x,phi(:,j),'k-'), xlabel( end hold on fprintf (1,'FINISHED CALCULATING %4i \n', j, ncutoff(j)); end hold off phi function volterra2_analytic % % Solution to a Volterra Equation of the Second Kind % % The Volterra Equations have four parameters % % the lower limit of integration, a % the constant prefactor outside the integral, lam % the function outside the integral, f, given in a function at the bottom of this file % the kernel, kernelorig, given in a function at the bottom of this file % % This code gives teh analytic solution to a=0, lam=-1, f=e=^x and N=e^(x-y) % % Author: David Keffer, University of Tennessee, Department of Chemical Engineering % clf time_01 = cputime; a = 0; lam = -1; % % set up a grid of x values % xstart = 0.0; xf = 1.0; dx = 0.01; nx =(xf-xstart)/dx + 1; x = zeros(nx,1); for i = 1:1:nx x(i) = xstart + (i-1)*dx; end % % let's run different runs where we truncate the resolvent % kernel at different points just to see the effect % nruns = 10; ncutoff=[1;2;3;4;5;6;7;8;9;10]'; phi = zeros(nx,nruns); % for j = 1:1:nruns % % calculate and plot purely analytic phi % for i = 1:1:nx xi = x(i); for k = 0:1:ncutoff(j) phi(i,j) = phi(i,j) + (-xi)^k/factorial(k); end phi(i,j) = phi(i,j)*f(xi); end function y = f(x) y = exp(x); function y = factorial(x) fact = 1; for i = 1:1:x fact = fact*i; end y = fact; 9 'x' ), ylabel ( 'phi' ) 'x' ), ylabel ( 'phi' ) 'x' ), ylabel ( 'phi' ) 'x' ), ylabel ( 'phi' ) 'x' ), ylabel ( 'phi' ) 'x' ), ylabel ( 'phi' ) 'x' ), ylabel ( 'phi' ) 'x' ), ylabel ( 'phi' ) 'x' ), ylabel ( 'phi' ) 'x' ), ylabel ( 'phi' ) 'x' ), ylabel ( 'phi' ) THE %2i RUN where ncutoff = D. Keffer, ChE 505 ,University of Tennessee, August, 1999 8. Numerical solutions to linear integral equations One can imagine different types of numerical solutions to a linear integral equation like the Volterra equation. Frequently we are going to encounter integral series, the which we have no intention of analytically evaluating either because we don t know how to evaluate the integrals analytically or it s not worth our time. Faced with this problem, we could code the analytical solution above, using numerical methods to numerically evaluate all of the integrals required for the kernel, the iterated kernels, and finally for phi. For example, we could use the trapezoidal rule to evaluate all of the iterated kernels in equation (20) (up to however many terms we want to use) and use the trapezoidal rule again to evaluate the solution for phi as given in equation (21). No two ways about it, THIS IS COMPUTATIONALLY INTENSIVE. It results in the evaluation of each iterated integral n2 times where n is the number of points along our axis. So if we want to evaluate phi at 100 x-values, using 5 iterated kernels, we will be forced to evaluate 50,000 integrals! This sucks. If the equation is nonlinear, we may be forced to do this. A second, much more elegant method is to take advantage of the linearity of the equation. We can use a numerical method for the evaluation of the integrals but plug this method directly into the integral equation. Consider Volterra s equation of the second kind x (x ) = f (x ) + N (x , y) ( y)dy (6) a and consider the trapezoidal rule for n intervals: b g( y)dy = a b a f (a ) n f (b) + f (y j) + n 2 2 j= 2 (35) If we substitute equation (35) in for the integral in equation (36) we have N (x i , a ) (a ) i 1 N( x i , x i ) ( x i ) (x i ) = f (x i ) + h + N( x i , y j ) ( y j ) + 2 2 j= 2 h i 1 N (x i , a ) N(x i , x i ) (a ) h N( x i , y j ) ( y j ) + 1 h ( x i ) = f ( x i ) 2 2 j= 2 (37) (38) for i = 1, equation (38) becomes, where x 1 = a (a ) = f (a ) (39) and for i = 2, we have h N ( x 2 , x1 ) N( x 2 , x 2 ) ( x1 ) + 1 h ( x 2 ) = f ( x 2 ) 2 2 (40) and for i = 3, we have h N( x 3 , x1 ) N( x 3 , x 3 ) ( x1 ) hN( x 3 , x 2 ) ( x 2 ) + 1 h (x 3 ) = f (x 3 ) 2 2 10 (41) D. Keffer, ChE 505 ,University of Tennessee, August, 1999 Now we have n equations of type equation (37), where x 1 = a and x n +1 = b . If you look at these equations, you will see that they form a set of linear algebraic equations! Whoa boy! We know how to solve that. The unknowns n unknowns are the function evaluated at each of our n points. Our problem is going to look like this: A = b (42) where = [ ( x1 = a ) ( x 2 ), ( x 3 ) K ( x i ) L ( x n ), ( x n +1 )] (43) A off ,i, j = hN ( x i , y j ) (44) N( x i , x i ) A diag ,i = h 1 2 (45) T 1 0 0 0.5A A diag , 2 0 off , 2 ,1 0.5A off ,3,1 A off ,3, 2 A diag ,3 A= 0.5A off ,i,1 A off ,i, 2 A off ,i ,3 A off ,n , 2 A off ,n , 3 0.5A off ,n ,1 0.5A off ,n +1,1 A off ,n +1, 2 A off ,n +1,3 0 0 0 A diag ,i A off ,n ,i A off ,n +1,i 0 0 0 0 A diag ,n A off ,n +1,n A diag ,n+1 0 0 0 0 0 (46) and where f ( x 1 = a ) f (x ) 2 f (x 3 ) b= f (xi ) f (x n ) f (x ) n +1 (47) Now, the accuracy of the solution depends on the step-size. As an illustration, we solve the same problem that we solved analytically in the last section. This time we solve it with this method for various values of the step-size. 11 D. Keffer, ChE 505 ,University of Tennessee, August, 1999 1 .0 0 5 n=16 & n=32 n=8 1 phi 0 .9 9 5 n=4 0 .9 9 0 .9 8 5 0 .9 8 n=2 0 .9 7 5 0 0 .1 0 .2 0 .3 0 .4 0 .5 x 0 .6 0 .7 0 .8 0 .9 1 Clearly as the number of intervals, n, increases, the accuracy of our solution to the integral equation also increases. In a result that should not surprise us, Linz has shown that should one use better methods for numerical integration than the trapezoidal rule, one will obtain more accurate results. 12 D. Keffer, ChE 505 ,University of Tennessee, August, 1999 amat(i,i) = 1.0; else amat(i,i) = 1.0 - fac*0.5*kernelorig(xi,xi) ; end bnum(i) = f(xi); for j = 1:1:npoints yj = x(j); if (j < i) if (j == 1) amat(i,j) = - fac*0.5*kernelorig(xi,yj); else amat(i,j) = -fac*kernelorig(xi,yj); end end end end amatinv = inv(amat); phi = amatinv*bnum % plot if (jj==1) plot (x,phi(:),'k-'), xlabel( 'x' ), ylabel ( 'phi' elseif (jj==2) plot (x,phi(:),'r-'), xlabel( 'x' ), ylabel ( 'phi' elseif (jj==3) plot (x,phi(:),'b-'), xlabel( 'x' ), ylabel ( 'phi' elseif (jj==4) plot (x,phi(:),'g-'), xlabel( 'x' ), ylabel ( 'phi' elseif (jj==5) plot (x,phi(:),'m-'), xlabel( 'x' ), ylabel ( 'phi' elseif (jj==6) plot (x,phi(:),'k:'), xlabel( 'x' ), ylabel ( 'phi' elseif (jj==7) plot (x,phi(:),'r:'), xlabel( 'x' ), ylabel ( 'phi' elseif (jj==8) plot (x,phi(:),'b:'), xlabel( 'x' ), ylabel ( 'phi' elseif (jj==9) plot (x,phi(:),'g:'), xlabel( 'x' ), ylabel ( 'phi' elseif (jj==10) plot (x,phi(:),'m:'), xlabel( 'x' ), ylabel ( 'phi' else plot (x,phi(:),'k-'), xlabel( 'x' ), ylabel ( 'phi' end hold on fprintf (1,'FINISHED CALCULATING THE %2i RUN where nintervals = %10f \n', jj, ninterval); end hold off function volterra2_analytic % % Solution to a Volterra Equation of the Second Kind % % The Volterra Equations have four parameters % the lower limit of integration, a % the constant prefactor outside the integral, lam % the function outside the integral, f, given in a function at the bottom of this file % the kernel, kernelorig, given in a function at the bottom of this file % % This code gives the numeric solution to a=0, lam=-1, f=e=^x and N=e^(x-y) % Author: David Keffer, University of Tennessee, % Department of Chemical Engineering clf time_01 = cputime; a = 0; lam = -1; % % set up a grid of x values % xstart = 0.0; xf = 1.0; % % let's run different runs where we change the step size % just to see the effect % nruns = 5; nintervals=[2;4;8;16;32;64;128;256;512;1024]'; npoints = nintervals+1; for jj = 1:1:nruns % % set up the x grid % ninterval = nintervals(jj); npoints = ninterval+1; phi = zeros(npoints,1)'; dx = (xf-xstart)/ninterval; x = zeros(npoints,1); for i = 1:1:npoints x(i) = xstart + (i-1)*dx; end % % calculate and plot purely analytic phi % amat = zeros(npoints,npoints); bnum = zeros(npoints,1); fac = lam*dx; for i = 1:1:npoints xi = x(i); if (i == 1) function y = f(x) y = exp(x); function y = kernelorig(x,y) y = exp(x-y); 13 ) ) ) ) ) ) ) ) ) ) ) D. Keffer, ChE 505 ,University of Tennessee, August, 1999 9. Numerical solutions to nonlinear integral equations There are many ways to skin a cat. There are at least as many ways to solve nonlinear integral equations. I will speak briefly on two methods: (1) an extension to the linear method of Section 8, and (2) the method of successive approximations. Consider a general nonlinear integral equation x (x ) = f ( x ) + F [ x, y, ( y)]dy (4) a where F [ x , y, ( y)] is a nonlinear function of ( y) , for example F [ x , y, ( y)] = sin[ ( y)] (5) First, in Section 8, we showed how to reduce a linear integral equation to a system of linear algebraic equations, which we then know how to solve. It should come as no surprise that we can reduce a nonlinear integral equation into a system of nonlinear algebraic equations, which we also know how to solve, although it frequently is a pain to do so. Still, if we go through the trouble to code it up once, then we have it for time immemorial, so long as there continue to be FORTRAN compilers. Since FORTRAN will surely outlast us, we are alright. We can discretize the x-axis into n intervals. Thus we obtain n+1 unknowns ( x i ) for i = 1 to n+1. If we are smart, we choose to discretize the y-axis (the dummy axis of integration) in the exact same way that we discretized the x-axis. Then, we can write something like, again using Trapezoidal rule, (x 1 ) = f ( x 1 = a ) (48) (x 2 ) = f (x 2 ) + h [F [ x 2 , x 1 , ( x1 )] + F [ x 2 , x 2 , ( x 2 )]] 2 (49) (x 3 ) = f ( x 3 ) + h [F [ x 3 , x 1 , ( x 1 )] + 2F [ x 3 , x 2 , (x 2 )] + F [x 3 , x 3 , (x 3 )]] 2 (50) j= i 1 h (x i ) = f ( x i ) + F [x i , x 1 , ( x 1 )] + 2 F [ x i , x j , ( x j )] + F [ x i , x i , ( x i )] 2 j=2 (51) This is a system of nonlinear algebraic equations. We could solve this using any technique in our arsenal for solving systems of nonlinear algebraic equations. This case is particularly amenable to the Newton-Raphson method since all of the partial derivatives in the Jacobian are going to have the same functional form. Remember the derivatives are taken with respect to ( x i ) and not with respect to x i . The second method for solving a nonlinear integral equation is the method of successive approximations. Again, referring to Pogorzelski [page 192], we see that equation (4) can be solved using the recursive relation x i ( x ) = f ( x ) + F [ x , y, i 1 ( y)]dy (52) a 14 D. Keffer, ChE 505 ,University of Tennessee, August, 1999 You repeat this procedure until i +1 ( x ) i ( x ) is acceptably small. The only additional information needed for this method is the starting point 0 (x) = f (x) (53) Now, certain qualifications have to be made about when this method works. The solution has to satisfy the Lipschitz condition and various other criteria have to be met. Well, at this point, I wave my hands and leave this to the mathematicians. We are engineers and we are here to use the equations. If we try and use it and it doesn t converge, then we go to the other method of solutions! Mathematicians please forgive us ignorant and sacrilegious engineers. Now, if the truth be told, 0 ( x ) can be any arbitrary continuous function in our range of interest. We set it to f(x) just for the heck of it. You could just as easily set it to unity. Let s work an example and solve it three ways, using (1) an extension to the linear method of Section 8, and (2) the method of successive approximations with 0 ( x ) = f ( x ) and (3) the method of successive approximations with 0 ( x ) = 1.0 . x (x ) = e x sin[ ( y)]dy (54) 0 Here we plot the solution to equation (54) using successive approximations. 1.9 n=1 n=2 and 3 1.8 1.7 1.6 phi 1.5 1.4 1.3 1.2 1.1 1 0 0.1 0.2 0.3 0.4 0.5 x 0.6 0.7 0.8 0.9 In this figure, n is the number of iterations performed. In this case, 0 ( x ) = f ( x ) = e . x 15 1 D. Keffer, ChE 505 ,University of Tennessee, August, 1999 1.9 n=1 n=2 and 3 1.8 1.7 1.6 phi 1.5 1.4 1.3 1.2 1.1 1 0 0.1 0.2 0.3 0.4 0.5 x 0.6 0.7 0.8 0.9 1 In this figure, n is the number of iterations performed. In this case, 0 ( x ) = 1.0 . From this we can see that, for this particular equation, the approximation converges pretty dang quickly. It is left up to the student to verify these results using the other method, namely that of using a method for solving a system of nonlinear algebraic equations. Below is the MATLAB code used to implement the successive approximations solution method. 16 D. Keffer, ChE 505 ,University of Tennessee, August, 1999 kernphi1 = kern(phi(1,jjm1)); integral_mid = 0.0; phi(1,jj) = f(x(1)); for i = 2:1:npoints integral_mid = integral_mid + kern(phi(i,jjm1)); integral = dxh*( kernphi1 + 2.0*integral_mid kern(phi(i,jjm1)) ); phi(i,jj) = f(x(i)) + lam*integral; end % % plot % if ( jjm1 == ncutoff(jrun) ) jrun = jrun + 1; if (jj==1) plot (x,phi(:,jj),'k-'), xlabel( 'x' ), ylabel ( 'phi' ) elseif (jj==2) plot (x,phi(:,jj),'r-'), xlabel( 'x' ), ylabel ( 'phi' ) elseif (jj==3) plot (x,phi(:,jj),'b-'), xlabel( 'x' ), ylabel ( 'phi' ) elseif (jj==4) plot (x,phi(:,jj),'g-'), xlabel( 'x' ), ylabel ( 'phi' ) elseif (jj==5) plot (x,phi(:,jj),'m-'), xlabel( 'x' ), ylabel ( 'phi' ) elseif (jj==6) plot (x,phi(:,jj),'k:'), xlabel( 'x' ), ylabel ( 'phi' ) elseif (jj==7) plot (x,phi(:,jj),'r:'), xlabel( 'x' ), ylabel ( 'phi' ) elseif (jj==8) plot (x,phi(:,jj),'b:'), xlabel( 'x' ), ylabel ( 'phi' ) elseif (jj==9) plot (x,phi(:,jj),'g:'), xlabel( 'x' ), ylabel ( 'phi' ) elseif (jj==10) plot (x,phi(:,jj),'m:'), xlabel( 'x' ), ylabel ( 'phi' ) else plot (x,phi(:,jj),'k-'), xlabel( 'x' ), ylabel ( 'phi' ) end hold on fprintf (1,'THE %2i RUN where nintervals = %4i and iteration = %4i, phi(%5.2f) = %5.2f \n', jjm1, ninterval, ncutoff(jrun1), x(npoints), phi(npoints,ncutoff(jrun)) ); end end hold off function volterra2_nonlinear % % Solution to a Nonlinear Volterra Equation of the Second Kind % % The Volterra Equations have four parameters % % the lower limit of integration, a % the constant prefactor outside the integral, lam % the function outside the integral, f, given in a function at the bottom of this file % the kernel, kernelorig, given in a function at the bottom of this file % This code gives the numeric solution to a=0, lam=-1, f=e=^x and K=sin(phi) % % Author: David Keffer, University of Tennessee, Department of Chemical Engineering % clf a = 0; lam = -1; % % set up a grid of x values % xstart = 0.0; xf = 1.0; % % USING SUCCESSIVE APPROXIMATIONS % ninterval = 20; npoints = ninterval+1; dx = (xf-xstart)/ninterval; x = zeros(npoints,1); for i = 1:1:npoints x(i) = xstart + (i-1)*dx; end nruns = 3; ncutoff=[1;2;3;4;5;6;7;8;9;10]'; phi = zeros(npoints,ncutoff(nruns))'; % % initialize phi(:,1) % for i = 1:1:npoints % phi(i,1) = f(x(i)); phi(i,1) = 1.0; end % % start iteration loop % jrun = 1; for jj = 2:1:ncutoff(nruns)+1 jjm1 = jj-1; dxh = 0.5*dx; function y = f(x) y = exp(x); function y = kern(x,y) y = sin(x); 17 D. Keffer, ChE 505 ,University of Tennessee, August, 1999 10. Numerical solutions to systems of integral equations We know that the formalism for solving systems of equations is often precisely the same as that for solving a single equation. We have seen this again and again. With algebraic equations, we can use the NewtonRaphson method to solve a single non-linear algebraic equation. Similary, we can use the multivariate NewtonRaphson method to solve a system of non-linear algebraic equations. With ODEs, we use Runge-Kutta to solve 1 ODE and multivariate Runge-Kutta to solve a system of ODEs. With PDEs, the same trend holds. With IEs the same trend still holds. For linear integral equations, we reduced a single IE to a system of n linear algebraic equations, where n was the number of intervals used in our discretization of the x-axis. For a system of linear integral equations, of the form: m i ( x ) = f i ( x ) + i N i , j (x , y) j ( y) dy a j=1 x (55) where i ranges from 1 to m we then have a system of mn equations, where n is again the number of intervals used in our discretization of the x-axis, and m is the number of integral equations. I stop here because the derivation of the method is totally analogous to what was given above. For nonlinear equations, we can use successive approximation on a system of nonlinear IEs just as easily as we used it on a single nonlinear IE. We have to simultaneously successively approximate all m i ( x ) based on the previous m, but so what? That s a piece of cake. The extension to systems requires no new intellectual effort whatsoever. 18

Formatting page ...

gangu chat

Horizontal lines at:
Guest Horizontal lines at:
AutoRM Data:
Box geometries: