Quadrature CS3220 - Summer 2008 Jonathan Kaldor
What is Quadrature? • Numerical evaluation / approximation of definite integrals • ∫ ab f(x) d x
What is Quadrature? • Recall from calculus that we can compute the integral of a function f(x) by computing its antiderivative F(x) [ the function such that F’(x) = f(x) ] • Fundamental Theorem of Calculus then says that ∫ ab f(x) dx = F(b) - F(a)
What is Quadrature? • This is of course extremely powerful and useful, but it has a downside: computing the antiderivative symbolically can be difficult • Computers (and calculators) have gotten better at it • Still, there are plenty of integrals out there that exist but are difficult / impossible to evaluate this way
What is Quadrature? • Consider the integral ∫ ab e x2 dx • This integral exists, but we cannot find an elementary function F(x) such that F’(x) = e x2 • As a result, if we want to evaluate this integral, we must do so numerically
What is an Integral? • So we already know that a definite integral is just the antiderivative evaluated at the endpoints • It also corresponds to a more fundamental notion: the (signed) area under a curve
What is an Integral?
What is an Integral?
What is an Integral? • So we can approximate the integral by computing an approximation of the (signed) area under the curve • We would like to compute this approximation using only evaluations of f(x) • We would also like to be able to reduce the amount of error by computing f(x) more times
Applications • By now, the usual suspects • Structural Analysis (Finite Element Analysis) • Physics (Electricity and Magnetism) • Physical Simulation (Forces and Energy) • Also: rendering (have to evaluate some really nasty integrals efficiently)
What is an Integral? • Can use an extremely important property of integrals to refine the interval over which we are integrating: ∫ ab f(x) d x = ∫ ac f(x) d x + ∫ cb f(x) d x • This allows us to break up the interval into pieces and approximate each piece separately (called composite rules )
Midpoint Approximation • Easiest approximation: approximate area under curve using a rectangle whose width is equal to the interval width and height is equal to the function evaluated at some point • If we evaluate it at the midpoint of our interval, we get the midpoint approximation
Midpoint Approximation • Area of rectangle is just base times height, so we get our estimated integral as M = (b-a) f((a+b)/2) • Note that if f((a+b)/2) is less than zero, we end up with negative area, which is what we want
Midpoint Approximation
Midpoint Approximation • What error do we have in our approximation? • One way to consider error is to look at the smallest degree polynomial we cannot integrate exactly • Why polynomials? We can compute their exact answer analytically
Midpoint Approximation • For midpoint, we can exactly integrate constant functions and linear functions, but not quadratic or above • The order of this method is thus 2 • Why does this work for linear functions?
Midpoint Approximation • Another way of deriving midpoint is to make the assumption that the function is constant over the entire interval. • What happens when we assume the function is linear over the interval? • Easiest approach: assume it is linear and passes through the endpoints of the interval
Trapezoid Rule
Trapezoid Rule • Definition of Trapezoid Rule is then T = (b-a) (f(a) + f(b)) / 2 • By definition, this integrates linear functions exactly, but not quadratics (order 2)
Extending it Further • We’ve approximated our function using constant pieces and linear pieces... can certainly use higher order approximations as well, which require more samples. Easiest method is to use evenly spaced sample points • Quadratic pieces (requires three points): S = (b-a)/6 (f(a) + 4 f((a+b)/2) + f(b))
Simpsons Rule • This is known as Simpsons Rule. If we break our integral into [a,c] and [c,b], and apply Simpsons Rule on each of the two halves with midpoints d and e, we get (b-a)/12 (f(a) + 4f(d) + 2f(c) + 4f(e) + f(b)) which is a composite rule. Note the 1 4 2 4 2 ... 4 2 1 repetition of function evaluations
Simpsons Rule • Clearly, this integrates quadratic functions exactly. It turns out that this also integrates cubic functions exactly as well (but not quartics or higher) • We can of course use higher degree approximations, but they quickly become unwieldy (as we’ve seen, high degree polynomials have some undesirable properties)
Simpsons Rule • We can also express Simpsons Rule in terms of our midpoint and trapezoid rules: S = (b-a)/6 (f(a) + 4 f((a+b)/2) + f(b)) = 2/3 M + 1/3 T
Error in Quadrature Rules • You might expect that the error in our approximation decreases as we improve our approximation from constant functions to linear, then quadratic, etc • Simpsons rule does have less error than either Midpoint or Trapezoid, but Midpoint tends to have less error than Trapezoid, somewhat surprisingly (about one half the error)
Error in Quadrature Rules • Your book observes that Midpoint seems to have half the error of Trapezoid and that the errors are opposite signs, but doesn’t explain why • Follows from an analysis of the Taylor expansion of the function around the midpoint (a+b)/2
Error in Quadrature Rules • f(x) = f(m) + f’(m)(x-m) + f’’(m)/2 (x-m) 2 +... • We can integrate this function and evaluate it over [a,b]: • ∫ f(x) = f(m) x + f’(m)/2 (x-m) 2 + f’’(m)/6 (x-m) 3 + ... = f(m) b - f(m) a + f’’(m)/6 (b-a) 3 + ... = (b-a)f(m) + f’’(m)/6 (b-a) 3 + ... = M + E 2 + ...
Error in Quadrature Rules • So the actual integral is our midpoint evaluation, plus some terms that depend on the second (and fourth, and sixth, and...) derivative. • This allows us to formalize why midpoint exactly integrates constant and linear functions (those functions have no second- and-higher derivatives, so midpoint is exact)
Error in Quadrature Rules • We can also use the Taylor series expansion to determine the error in the Trapezoid rule. • Unsurprisingly, the terms involving the first derivative drop off as well (since we can integrate constant and linear functions exactly), but the terms involving the second derivative do not • End up with ∫ f(x) = T - 2 E, where E is the same term as in Midpoint
Error in Quadrature Rules • Thus, we have a mathematical explanation for a.) why midpoint is the same order as trapezoid, and b.) why midpoint has about half the error, and opposite sign, of the trapezoid rule • We can then use both rules to estimate this error: E ≈ (T - M) / 3
Error in Quadrature Rules • We can also try and use both Midpoint and Trapezoid rules to eliminate the dominant error term • Since trapezoid has twice the dominant error E and opposite sign, we want to add midpoint and trapezoid together, scaling midpoint by 2 and dividing the entire term by 3
Error in Quadrature Rules • ∫ f = M + E + ... ∫ f = T - 2E + ... • 2 ∫ f = 2M + 2E + ... ∫ f = T - 2E + ... • Add two equations together: 3 ∫ f = 2M + T + ... ∫ f = 2/3 M + 1/3 T + ... • Note: this is just Simpsons rule!
Composite Quadrature • As we have noted, we can further refine the interval [a,b] to improve our approximation, computing the integral over [a,c] and [c,b] for some choice of c • If we divide the interval in half, Trapezoid and Simpson’s Rules have the benefit of having their refined points of evaluation be a superset of the original points
Composite Quadrature
Composite Quadrature
Composite Quadrature
Composite Quadrature • Note that this is not true for midpoint - we need to recompute all midpoint values when we refine
Adaptive Quadrature • The other benefit is that by using two quadrature rules together, we can estimate the amount of error in our approximation • Use it as a subdivision rule to divide our interval if the error is above some tolerance ε • Only refine intervals that need additional terms - avoid refining smooth regions whose integral is well-approximated
Adaptive Quadrature • For instance, using both midpoint and trapezoid together, we get an estimate of the leading term error • If it is above some tolerance, refine intervals and recompute. Ends up focusing effort where it is needed
Adaptive Quadrature
Quadrature in MATLAB • Unsurprisingly, v = quad(f, a, b, tol) computes the quadrature of the function handle f over the closed interval [a,b] using an adaptive Simpson’s rule with tolerance tol • Also has several other quadrature methods, with varying performance on certain types of functions (quadl, quadgk)
Recommend
More recommend
