1、An Introduction to Programming and Numerical Methods in MATLAB-S.R. Otto & J.P. Denier-PAGE（234-254） （原文）Chapter 7 : Numerical Integration7.1 Introduction In this chapter we shall discuss techniques whereby functions can be integrated; these are quite classical. We will give the derivation of the ap

2、proximations and will mention the likely failings of the techniques. Two versions of the code will be given for each of the main methods, so we can start to appreciate the power of MATLAB. It is of course possible to write the code as if it were Fortran or C, but that would waste the power of the pa

3、ckage. Before we start we note: if we consider two points we can fit a straight line through them; with three we can fit a quadratic and with four we can fit a cubic (this was discussed in more detail in Chapter 5). We start by breaking down the integration regime into small intervals and approximat

4、ing the area below the curve by slices. This method is tantamount to counting the squares and dealing with the parts of squares at the tops of the columns in sophisticated ways.In this case the area under the curve is approximately7 4 “the area of the boxes” (the fourth row up has a few part squares

5、 as does the fifth row). We obviously need a scheme which is slightly more robust.7.2 Integration Using Straight LinesIn this chapter our objective is to calculate the value of the integral=x=abfxdx.We shall assume we can calculate (easily) the value of f (x) for all values of x between a and b. Thi

6、s means we are dealing with a function rather than a set of data values. We shall explore the latter case in due course but at this juncture we shall use a new routine integrand.m:% % integrand.m % input a set of values (x) % output function values f(x)%function f = integrand(x) % Here we use f = si

7、n(x2) as a % sample function. f = sin(x.2); We shall now take a while to derive the method we are going to use to integrate the function on the grid of points. We shall use N points and as such we use the code step = (b-a)/(N-1); x = a:step:b; f = integrand(x); In order to derive the form for the in

8、tegration we introduce the nomenclature that the points we have defined above are (xj,fj), where j runs from 1 to N. Let us consider the consecutive points (xj,fj ) and (xj+1,fj+1) and approximate the curve between them by a straight line. We use the formula for the line through the points (x1,y1) a

9、nd (x2,y2) which isY-Y1Y2-Y1 = X-X1X2-X1 or y=y1 + Y2-Y1X2-X1 (x-x1). In our case this givesfL(x)=fj+fj+1-fjh(x-xj),where we have introduced h = xj+1 xj and fL(x) the formula for the straight line (This is just a direct application of Newtons Forward Differences, (5.1).) Let us now perform the analy

10、tical integration of the function fL(x) between xj and xj+1 to determine the area Aj,j+1:Aj,j+1=x=xjxj+1fLxdx=x=xjxj+1fj+1-fjhx-xjdx+fjdx. We introduce the linear transformation X= x xj, where dX = dx and when x = xj, x = 0 and x= xj+1 correspond to x= h. The integral becomesAj,j+1=X=0hfj+1-fjhX+fjd

11、X, which can be integrated to yieldAj,j+1=h2fj+fj+1. This is the area of a trapezium with vertices at (xj, 0), (xj, fj), (xj+1, fj+1) and (xj+1, 0). The area of a trapezium is the mean of the length of the two parallel sides times the perpendicular distance between them; that is (fj+fj+1)/2 times (x

12、j+1-xj). This method is unsurprisingly called the trapezium rule. In order to determine the total area from x= a to x = b we sum all the partsArea=j=1j=N-1Aj,j+1=j=1j=N-1h2fj+fj+1. In order to perform this summation it is instructive to write out the seriesArea= h2f1+f2+h2f2+f3+h2f3+f4+h2fN-2+fN-1+h

13、2fN-1+fN, so thatArea= h2f1+2f2+2f3+2fN-1+2fNx=abfxdx, Example 7.1 We shall calculate the integral of the function fx = x3sinx between zero and one. We shall use N points N = 10; x = linspace(0,1,N); h = x(2)-x(1); f = x.3.*sin(x); g = h*(sum(f)-f(1)/2-f(N)/2); Here we have constructed a grid of poi

14、nts running from zero to one and then worked out the gap between successive points (that is h). We now construct the function f and work out the expression for the trapezium rule, which is the sum of the values of f minus half of the end values. This gives the value of the integral as g.7.2.1 Errors

15、 in the Trapezium Method The number of points required for a calculation depends on how exactly you need to know the answer (in general). The error in this scheme is encountered because we approximate the curve between the points (xj,fj) and (xj+1,fj+1) by a straight line. This can be reduced by usi

16、ng a quadratic instead over the interval spanned by the requisite three points. We could have used the fundamental code which did not make use of the combinations of the terms: integral = 0; for i = 2:N integral = integral+(f(i)+f(i-1)/2*h; end The advantage of this form of the code is it is far eas

17、ier to convert to one in which intermediate results are available. We note that the command /2*h divides by two and multiplies by h rather than dividing by 2h: this would be accomplished via parentheses, that is /(2*h). Before we proceed we should consider the error involved in approximating the are

18、a by a set of trapezia. The error is the unshaded part below the curve. To estimate the extent of this we again use “Taylor Series” and note thatfx=fxj+fxj+1-fxjhx-xj+x-xjx-xj+12d2fdx2|x=E where E xj, xj+1. This can be integrated to givex=xjxj+1fxdx=h2fxj+fxj+1+h36d2fdx2|x=E Consequently the error i

19、n using just the rst term is proportional to h3 and fE. Notice that if the second derivative is zero over the range the error is actually zero. Unsurprisingly this corresponds to fx being a straight line.7.3 Integration Using Quadratics We need to construct a curve which passes through the points (x

20、j-1,fj-1 ), (xj,fj) and (xj=1,fj+1). Let us consider a value of j = 1: this is purely to reduce the verbosity of our expressions. The quadratic passing through these three points can be written asfq=f0+f0x-x0h+2f0x-x0x-x12h2,using a truncation of Newtons forward difference formula, (5.1).In passing

21、we mention this series expansion is similar to a Taylor series expansion but instead of the terms being differentials they are differences. In order to derive the formula we now need to integrate the function fqx from x=x0 to x=x2. This is quite straightforward: however at this point we will exploit

22、 the symbolic capabilities of MATLABsyms x f0 f1 f2 h x0q = f0+(f1-f0)/h*(x-x0)+(f2-2*f1+f0)/(2*h2)*(x-x0)*(x-(x0+h);iq = int(q,x0,x0+2*h);simplify(iq)ans =1/3*h*(f0+4*f1+f2)We shall dissect this portion of code so you can appreciate what is happening. In the first line we assign the variables x,f0,

23、f1,f2,h and x0 to be symbolic. This means that MATLAB does not have to know the value of the variables and treats them as mathematical objects. In the second line we set up the function q (which is just the quadratic fqx). This is integrated in the third line (between x0 and x2=x0+2*h. (This produce

24、s a verbose answer.) Finally we simplify our answer.Hence we have the answerx=x0x=x2fqxdx=h3f0+4f1+f2.Although we know there is only one quadratic through a given set of points, it is instructive to re-derive this in another way. We shall start with the general quadratic.fqx=a+bx-x0+cx-x0x-x1.The re

25、quirement that the quadratic goes through the point (x0,f0) gives us that a=f0. Using the point (x1,f1) gives us b=f1-f0x1-x0.Finally using (x2,f2) we havef2=f0+f1-f0x1-x0x2-x0+cx2-x0x2-x1.which can be manipulated to give:c=1x2-x1f2-f0x2-x1-f1-f0x1-x0.This gives us the general form of the quadratic

26、through three points (which can be used for Task 7.10). We now return to the form for regularly spaced points so xj=h and integrating from x=x0 to x=x2=x0+2h:x=x0x=x0+2hfqxdx=x=x0x=x0+2hf0+f1-f0hx-x0+1hf2-f02h-f1-f0hx-x0x-x1dxNow using the substitution X=x-x0, so that x=x0 corresponds to X=0 and x=x

27、2=x0+2h corresponds to X=2h the expression x-x1=X+x0-x1=X-h and dx=dX. Hence we have=X=02hf0X+f1-f0hX+12h2f2-f0-2f1-f0XX-hdX=f0X+f1-f0hX22+12h2f0-2f1+f2X33-hX2202h=2hf0+f1-f0h2h2+12h2f0-2f1+f28h33-2h3=h2f0+2f1-f0+13f0-2f1+f2=h3f0+4f1+f2We now need to divide the range of integration into the appropri

28、ate number of subintervals. Notice the number of intervals needs to be even (and hence the number of points needs to be odd). The total integral is approximated byx=abfxdxh3f0+4f1+f2+h3f2+4f3+f4+h3f4+4f5+f6+h3fN-4+4fN-3+fN-2+h3fN-2+4fN-1+fN.This can be simplified to givex=abfxdxh3f0+4f1+2f2+4f3+2f4+

29、2fN-2+4fN-1+fN.This is called Simpsons 1/3 rule.We shall now construct a code to determine the integral. At this point we could use a code which used a conventional approach, but we shall try to use a version which exploits the power of MATLAB.% Simpsons 1/3 rule.%x = 0.0:0.1:1.0;h = x(2)-x(1);N = l

30、ength(x);if mod(N,2) = 0disp(Routine needs an odd number of points)break % Ensure the number of points is oddendrodd = 1:2:N;reven = 2:2:(N-1);weights(rodd) = 2; weights(1) = 1; weights(N) = 1;weights(reven) = 4;f = sin(x.2);integral = h/3*sum(weights.*f);format long edisp(integral)This calculates t

31、he value of the expression01sinx2dx.The value which MATLAB comes out with is 3.2209e-01, where we have changed the way in which these numbers are displayed by using the command format long e. This answer can be checked using MATLABs symbolic capabilities:syms x ff = sin(x2);f1 = int(f,0,1)f1 =1/2*Fr

32、esnelS(2(1/2)/pi(1/2)*2(1/2)*pi(1/2)This value can be compared with that attained using the symbolic toolbox, using the command vpa (variable precision arithmetic):vpa(f1)ans =.(Note we can change the number of digits using the command digits(10). We further note that this is still a symbolic object

33、: in order to obtain a value (which can be used for plotting, for example) we use the command double. This is a case of MATLAB being too clever: it has solved the integral and written it as a Fresnel integral (for further details see Abramowitz and Stegun Handbook of Mathematical Functions). Our sim

34、ple integration using ten points does quite well: in fact the error is proportional to h4, which can be shown using the same technique as we used for the trapezium rule on page 230.In the above code we have introduced the term weights which is applied to the terms before they are added together to o

35、btain the integral. The MATLAB command mod allows the user to determine the remainder when the rst argument is divided by the second. In this case when considering mod 2 we are checking for parity (that is whether N is even or odd).In the code for Simpsons one third rule we have re-introduced anothe

36、r MATLAB command, namely break. This stops the code, or more exactly it exits the current level: for example if it is used in a nested loop it will terminate the current level and return to the previous level. This manual entry for break is instructive here help breakBREAK Terminate execution of WHI

37、LE or FOR loop.BREAK terminates the execution of FOR and WHILE loops.In nested loops, BREAK exits from the innermost loop only.We now progress to consider cubic approximations to the function in the hope that this will provide even more accurate answers.Again we can use a technique which involves su

38、mming the separate intervals, which you may prefer:integral = 0;for j = 1:2:N-2integral = integral + h*(f(j)+4*f(j+1)+f(j+2)/3;end7.4 Integration Using Cubic PolynomialsAs in the previous sections we need to define an approximating curve and in order to do so we need four points (x0, f0), (x1, f1),

39、(x2, f2) and (x3, f3). Using the same form as above we can write the cubic equation asfcx=f0+f0x-x0h+2f0x-x0x-x12h2+3f0x-x0x-x1x-x26h3.In addition to the terms defined earlier we have introduced the third-order forward difference. This can be defined recursively using 2, so that3f0=2f0=2f1-2f0=f3-2f

40、2+f1-f2-2f2+f0=f3-3f2+3f1-f0,using the Newton forward difference for 2f0. We now need to integrate from x=x0 to x=x3, and we again exploit the symbolic tool box for this, using the codesyms x f0 f1 f2 f3 h x0x1 = x0+h; x2 = x0+2*h; x3 = x0+3*h;t1 = (f1-f0)/h*(x-x0);t2 = (f2-2*f1+f0)/(2*h2)*(x-x0)*(x

41、-x1);t3 = (f3-3*f2+3*f1-f0)/(6*h3)*(x-x0)*(x-x1)*(x-x2);q = f0+t1+t2+t3;q1 = int(q,x,x0,x3);simplify(q1)This gives the answer 3/8*h*(f0+3*f1+3*f2+f3). Thus we havex=x0x3fcxdx=3h8f0+3f1+3f2+f3.You should note it is possible to do all these integrals by hand, but we wish to demonstrate the power and t

42、he utility of this particular toolbox. We now need to combine all the subintervals, sox=abfxdx3h8f0+3f1+3f2+f3+3h8f3+3f4+3f5+f6+3k8fN-3+3fN-2+3fN-1+fN,which can be simplified to givex=abfxdx3h8f0+3f1+3f2+2f3+3f4+3f5+2fN-3+3fN-2+3fN-1+fN.This is referred to as Simpsons 3/8rule. We shall now give a MA

43、TLAB program based on the same structure as the previous one. Notice this time the number of points needs to be divisible by three. % % simpsons 3/8 rule. % N = 10; x = linspace(0,1,N); h = x(2)-x(1); ms = Number of intervals should be divisible by three;if mod(N-1,3) = 0disp(ms)breakendm = (N-1)/3;

44、rdiff = 3*(1:(m-1)+1;weights = 3*ones(1,N);weights(1) = 1; weights(N) = 1;weights(rdiff) = 2;f = sin(x.2);integral = 3*h/8*sum(weights.*f);disp(integral)The answer given by this is 0.964. In fact this is not quite as good as the previous method. The error is again proportional to h4 but the constant

45、 of proportionality is larger. We could also write the integral in separate regions, that is without combination, sointegral = 0;for j = 1:3:N-3integral = integral + h*(f(j)+3*f(j+1) .+3*f(j+2)+f(j+3)*3/8;endEach of these methods has a restriction on the numbers of points one can use. However these

46、can be circumvented by using an amalgamation of the two schemes. There are other methods available for integrating functions and some of these will be met in due course. We could continue this process, especially with the symbolic toolbox at our disposal to perform the algebra. For instance the form

47、ula obtained by using a quartic over five points isx=x0x4fxdx2h457f0+32f1+12f2+32f3+7f4.Although you might think the higher the order the polynomial the more accurately you would know the answer, there are problems which are intrinsic to using high-order polynomials. In fact the optimum method is to

48、 use a combination of the two Simpson methods. If the number of points supplied is even we use the one third rule for the first N 3 points (which is necessarily an odd number of points and consequently an even number of intervals) and then we use the three-eighths rule on the remaining points.7.5 In

49、tegrating Using MATLAB CommandsAs with many examples in this text we can also use standard MATLAB commands, for instance quad and quad8 (see help quad). These commands use similar techniques to those above but with the advantage of automation. For instance they exploit grids which can adapt. This me

50、ans that in regions which are harder to integrate (perhaps with more variation in the function) the scheme adds extra points.The syntax is:tol = 1e-4 1e-5;a = 0; b = pi;trace = 1;q = quad(sin,a,b,tol,trace)These quantities are respectively the extent of the domain x a, b, the tolerances (relative an

51、d absolute) and whether the user wants to see a trace or not (setting trace as non-zero shows the evolution of the calculation).Example 7.2 We show the integration of the function j1x for x = 0 to x = 10. This is actually a Bessel function which occurs as one of the solutions to the differential equ

52、ation and is revisited in Chapter 8.x2d2ydx2+xdydx+x2-1y=0.Fortunately MATLAB has a routine which evaluates Bessel functions but we need to write our own routine to make sure that it is available for quad:function value=ourbess(x)value = besselj(1,x);and then we simply need the code:q = quad(ourbess

53、,0,10,1e-5 1e-5,1);This gives a value of 1.673. The integral of the function j1x is -j0x and hence the value we seek is -j010+j00 1.; so as we can see the integration scheme does very well.7.6 Specific Examples of IntegralsWe shall now describe how we can deal with other problems which arise when ev

54、aluating numerical integrals.7.6.1 Infinite Integrals and Removable Singularities Using the various methods we can evaluate integrals from a to b, but if one or more of these values is infinite we need a different treatment. This will vary depending on the form of the integral or more exactly the in

55、tegrand.Example 7.3 For example let us consider the integralI=0e-x2dx,In this case the integrand e-x2 decays very quickly so we can adequately determine the value of the integral usingIX=0Xe-x2dx,for a suitable value of X (in this case X = 10 is more than sufficient). We can evaluate the integral fo

56、r a few values of X until IX is constant (to within a defined tolerance). We could develop an algorithm to decide what value of X to use, but in general we will use the above method. We can also exploit the symmetry of problems, for instance-e-x2dxdx=20e-x2dx20xe-x2dx,provided X 1. We can also use t

57、ransformations to rescale the regions of integration.In fact in this example we could have used the MATLAB command erf which gives the error functionEq=2q=0xe-q2dq,We can use erf(Inf) which gives us the value unity.We now consider how we might perform an integral for which the integrand is singular at an end-point of the range. For instance:q=011q1/2dq=2q1/2q=01=2,Example 7.4 We