User:John Randall/Quadrature

From J Wiki
Jump to: navigation, search

Let f be a continuous function on [a,b]. Let n be an integer, h=(b-a)/n, x_i=a+ih, f_i=f(x_i). Then we can estimate \int_a^b f(x)\,dx as the linear combination \sum c_i f_i for some choice of coefficients c_i.

In the case n=6 we have, for example \sum c_if_i=\frac{h}2(f_0+2f_1+2f_2+2f_3+2f_4+2f_5+f_6)\mbox{ (trapezoidal rule)} \sum c_if_i=\frac{h}3(f_0+4f_1+2f_2+4f_3+2f_4+4f_5+f_6)\mbox{ (Simpson's rule)} Our aim is to explain where the magic coefficients come from.

Suppose the points \{x_i\} are divided into p contiguous panes, each with (m+1) points, so that n=pm. For the trapezoidal rule, we have p=6, m=1, and for Simpson's rule we have p=3, m=2. Consider a single pane with (m+1) points, and suppose a=0, b=m, h=1. The general result will follow from the linearity properties of the integral. Let L_i(x) be the Lagrange basis polynomials of degree m satisfying L_i(j)=\delta_{ij} for 0\leq i,j\leq m. Then f(x)\simeq\sum L_i(x) f_i \int_a^b f(x)\,dx \simeq \sum \left(\int_a^b L_i(x)\,dx\right)f_i=\sum c_i f_i

require 'numeric'
V=:^/~ @: i.@:#
L=:clean@(%. V)"1
basis=:L @: =@i.@>:
int=:(0&p.. p. <:@:#)"1
   c 1
0.5 0.5
   c 2
0.333333 1.33333 0.333333

For a single pane, these give

\int_{x_0}^{x_1} f(x)\,dx\simeq \frac{h}2(f_0+f_1)\mbox{ (trapezoidal rule)} \int_{x_0}^{x_2} f(x)\,dx\simeq \frac{h}3(f_0+4f_1+f_2)\mbox{ (Simpson's rule)} The result for multiple panes is simple additivity of the single pane result:

tz=:#~ -.@(*./\.)@(0&=)  NB. remove trailing zeros
ppr=:+//.@(*/)           NB. polynomial product
coeff=:(c@:]) pprtz 0 = ] | i.@:*
   6 coeff 1
0.5 1 1 1 1 1 0.5
   3 coeff 2
0.333333 1.33333 0.666667 1.33333 0.666667 1.33333 0.333333

The same method works for generating formulas with more interpolating points (although there are good reasons for stopping at Simpson's rule). If we take m=3, we get:

   2 coeff 3
0.375 1.125 1.125 0.75 1.125 1.125 0.375

\sum c_if_i=\frac{3h}8(f_0+3f_1+3f_2+2f_3+3f_4+3f_5+f_6)\ \ (Simpson's\ \frac38\ rule)