Quadrature CS3220 - Summer 2008 Jonathan Kaldor What is - - PowerPoint PPT Presentation

Quadrature

CS3220 - Summer 2008 Jonathan Kaldor

What is Quadrature?

• Numerical evaluation / approximation of

definite integrals

• ∫ab f(x) dx

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 ex2 dx

• This integral exists, but we cannot find an

elementary function F(x) such that F’(x) = ex2

• 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) dx = ∫ac f(x) dx + ∫cb f(x) dx

• 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

• ver 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
• ur 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
• f 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 + E2 + ...

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

• f functions (quadl, quadgk)