Download Numerical Integration: Approximating the Area Under a Curve and more Study notes Computer Science in PDF only on Docsity!
1 Numerical Integration
Integrating a function is crucial to many problems in science and engineering, and yet it is often the case that finding the analytic form of an integral is difficult. If what one has is just data that defines a function (data that might come from conducting an experiment, for example), then one doesn’t even have the actual function expressed in an analytic form. The solution to this is numerical integration, also called quadrature. In essence, if we can compute the value of a function given values of the argu- ments to that function, then we can compute the integral of that function, at least to some degree of numerical accuracy. It is important to mention that numerical integration is not the same thing as finding the area under the curve. If the function is always positive, then the two are the same, but one must remember that for functions that go negative the integral produces a “positive area” value and a “negative area” value that are summed to get the integral. This is an obvious thing and not conceptually difficult. The problem with numerical integration, unlike doing a calculus exercise, is that with a black box for the function one might not be aware that the function takes on negative values and thus might overlook the distinction.
1.1 The Simplest Method for Numerical Integration
Consider a function known to be positive over an interval xL to xR. (This means that we can speak of the integral and the area under the curve in- terchangeably.) In elementary calculus, one of the first things proven about integrals is that the integral is the limit of the Riemann sum, and the Rie- mann sum is the sum of rectangles that approximate the area under the curve. For example, the area S 0 = f (xL)(xR − xL) is the area of a rectangle of height f (xL) and width xR − xL. If the function is a constant y = f (xL), then we have that f (xL) = f (xR), the height of the function is a constant, and the area under the curve is exactly S 0. Most functions are not constants, though, so we have to continue our approximation. Let’s define x 0 = xL, x 1 = (xR − xL)/2, x 2 = xR, so we have broken the interval up into two intervals with left endpoints x 0 and x 1. Consider the sum S 2 = f (x 0 )(x 1 − x 0 ) + f (x 1 )(x 2 − x 1 ). Actually, it’s simpler to notice that Δ 2 = x 1 − x 0 = x 2 − x 1 and thus that our approximation is
S 2 = Δ 2 · (f (x 0 ) + f (x 1 )). We would hope that this is a better approximation to the area under the curve, because now we are using two points instead of just one to approximate the function. (Draw the picture, which I would do except that’s much more work.) Going further, let Δ 4 = (xR − xL)/4, define
x 0 = xL x 1 = xL + Δ 4 x 2 = xL + 2 · Δ 4 x 3 = xL + 3 · Δ 4 x 4 = xL + 4 · Δ 4 = xR
and let’s use the four rectangles with left endpoints x 0 , x 1 , x 2 , x 3 as our ap- proximations, so the integral will be approximated by
S 4 = (f (x 0 ) + f (x 1 ) + f (x 2 ) + f (x 3 )) · Δ 4
= (f (x 0 ) + f (x 1 ) + f (x 2 ) + f (x 3 )) ·
xR − xL 4
By now you are getting the picture. Namely, we let
ΔN = (xR − xL)/N,
compute the sum
SN =
(N − 1
i=
f (xL + i · ΔN )
ΔN
for increasing values of N. Note that we could just as easily have chosen to use the right endpoints of the subintervals instead of the left endpoints. As the subintervals get shorter and shorter, the difference between the two choices will disappear, so we lose no generality by simply sticking with the left endpoints. In truth, what we are doing with this approach numerically is exactly what is done in proving that the integral exists and is equal to the area under the curve for positive functions.
What this does is turn the previous sum
SN =
(N − 1
i=
f (xL + i · ΔN )
ΔN
into a better approximation
SN =
(N − 1
i=
f (xL + i · ΔN ) + f (xL + (i + 1) · ΔN ) 2
ΔN
If we were to use this on our linear function, we would get an exact answer every time, because the trapezoids we are summing up fix the “curve” exactly. In general, the trapezoid approximation will be better than the rectangular approximation because we will be including the little triangle on top of the rectangle to our computation of area.
1.4 Quadratic Approximations
A final comment is that we can go one step further and use quadratic poly- nomials for approximating the area under a curve. The trapezoid method is essentially looking at the function f (x) over an interval xi to xi+1 and approximating the function with the straight line between f (xi) and f (xi+1). The theorem that we use is the one that says that n + 1 points determine a unique polynomial of degree n. That is, two points uniquely determine a straight line, three points uniquely determine a quadratic, four points a cubic, and so forth. This is true unless there is some degeneracy involved; for example, if the “two” points are actually the same point, then we don’t get a line, and if three points actually lie on a straight line, then we don’t get a quadratic. But let us proceed as if we were working in the general case. Let’s assume we have three points (x 0 , y 0 ), (x 1 , y 1 ), and (x 2 , y 2 ). We can fit a quadratic to these three points by brute force solution of the equations for the coefficients. If our quadratic is f (x) = ax^2 + bx + c for unknowns a, b, and c, then we have ax^20 + bx 0 + c = y 0 ax^21 + bx 1 + c = y 1 ax^22 + bx 2 + c = y 2
We have three linear equations in three unknowns, and with some effort we can solve these for a, b, and c. Or we can pull the rabbit out of the hat and use Lagrange interpolation. The function
f (x) = y 0
(x − x 1 )(x − x 2 ) (x 0 − x 1 )(x 0 − x 2 )
+y 1
(x − x 0 )(x − x 2 ) (x 1 − x 0 )(x 1 − x 2 )
+y 2
(x − x 0 )(x − x 1 ) (x 2 − x 0 )(x 2 − x 1 )
is a quadratic polynomial that takes on the value y 0 at x 0 , y 1 at x 1 , and y 2 at x 2. If we believe the theorem that three points determine a unique quadratic polynomial, then f (x) as defined immediately above must be that quadratic. And if we look at the algebra, we find that f (x) = ax^2 + bx + c for
a =
y 0 (x 0 − x 1 )(x 0 − x 2 )
y 1 (x 1 − x 0 )(x 1 − x 2 )
y 2 (x 2 − x 0 )(x 2 − x 1 )
b =
−y 0 (x 1 + x 2 ) (x 0 − x 1 )(x 0 − x 2 )
−y 1 (x 0 + x 2 ) (x 1 − x 0 )(x 1 − x 2 )
−y 2 (x 0 + x 1 ) (x 2 − x 0 )(x 2 − x 1 )
c =
y 0 (x 1 x 2 ) (x 0 − x 1 )(x 0 − x 2 )
y 1 (x 0 x 2 ) (x 1 − x 0 )(x 1 − x 2 )
y 2 (x 0 x 1 ) (x 2 − x 0 )(x 2 − x 1 ) To do a quadratic approximation for the area under the curve of a func- tion, we start with x 0 = xL, x 1 = x 0 +ΔN , and x 2 = x 0 +2ΔN as three points that determine a quadratic. We note that x 1 − x 0 = ΔN and x 2 − x 0 = 2ΔN , so with these equally spaced points we can simplify the formulas further.
a =
y 0 − 2 y 1 + y 2 2Δ^2 N
b = −
y 0 (x 1 + x 2 ) − 2 y 1 (x 0 + x 2 ) + y 2 (x 0 + x 1 ) 2Δ^2 N
c =
y 0 (x 1 x 2 ) − 2 y 1 (x 1 x 2 ) + y 2 (x 0 x 1 ) 2Δ^2 N
