Class 13. Numerical IntegrationSimple Monte Carlo Integration (NRiC §7.6)• Can use RNGs to estimate integrals.• Suppose we pick n random points x1, ..., xNuniformly in a multi-D volume V .• Basic theorem of Monte Carlo integration:ZVf dV ' V <f> ± Vs<f2> − <f>2N,where<f> ≡1NNXi=1f(xi) and <f2> ≡1NNXi=1f2(xi).• The ± term is a 1-σ error estimate, not a rigorous bound.• Previous formula works fine if V is simple.• What if we want to integrate a function g over a region W that is not easy t o samplerandomly?• Solution: find a simple volume V that encloses W and define a new function f(x),x ∈ V , such that:f(x) =(g( x) for all x ∈ W ,0 otherwise.E.g., suppo se we want to integrate g(x, y) over the shaded area inside area A below:b(x)Area ATo integrat e, take random samples over the whole rectangle, setf(xi, yi) =(g( xi, yi) yi≤ b(xi),0 otherwise,and computeZshaded areag( x, y) dx dy 'ANXif(xi, yi).1– Nifty example: π can be estimated by integratingp(x, y) =(1 x2+ y2≤ 1,0 otherwise,over a 2 × 2 square:π =Z1−1Z1−1p(x, y) dx dy'4NXip(xi, yi).– See NRiC for another worked example.• Optimization strategy: make V as close as possible to W , since zero values of f willincrease the rel ative error estimate.• Principal disadvantage: accuracy increases only as square root of N.• Fancier routines exist for faster convergence: NRiC §7.7–7.8 .• Monte Carlo techniques used in a variety of other contexts: anywhere statistical sam-pling is useful. E.g.,– Predicting motion of bodies with short Lyapunov times if starting positions andvelocities poorly known.– Determining model fit significance by testing the model against many sets ofrandom synthetic dat a with the same mean and variance.Numerical Integration (Quadrature)• NRiC §4.• Already seen Monte Carlo integration.• Can cast problem a s a differentia l equation (DE):I =Zbaf(x) dxis equivalent to solving for I ≡ y(b) the DE dy/dx = f(x) with the boundary condition(BC) y(a) = 0 .– Will learn about ODE solution methods next class.2Trapezoidal and Simpson’s rules• Have abscissas xi= x0+ ih, i = 0, 1, ..., N + 1.• A function f(x) has known values f (xi) = fi.• Want to integrate f (x) between endpoints a and b.• Trapezoidal rule (2-point closed formula):Zx2x1f(x) dx = h12f1+12f2+ O(h3f00),i.e., the area of a trapezoid of base h and vertex height s f1and f2.• Simpson’s r ule(3-point closed formula):Zx3x1f(x) dx = h13f1+43f2+13f3+ O(h5f(4)).Extended trapezoidal rule• If we a pply the trapezoidal rule N − 1 times and add the results, we get:ZxNx1f(x) dx = h12f1+ f2+ f3+ ... + fN −1+12fN+ O"(b − a)3f00N2#.• Big a dvantage is it builds on previous work:– Coarsest step: average f at endpoints a and b.– Next refinement: add value at midpoint to average.– Next: add values a t14and34points.– And so on. This is implemented as trapzd() in NRiC.More sophistication• Usually don’t know N in adva nce, so iterate to a desired accuracy: qtrap().• Higher-order method by cleverly adding refinements to cancel error terms: qsimp().• Generalization to order 2k (Richardson’s deferred approach to the limit): qromb().– Uses extrapolation methods to set h → 0.• For improper integrals, generally need open formulae (not evaluated at endpoints).• For multi-D, use nested 1-D