MATH 417 Lecture 13

« previous | Tuesday, February 25, 2014 | next »

Gaussian Rule

for $n$ points $x_{1},x_{2},\ldots ,x_{n}$, our degree of accuracy is $2n-1=2(n-1)+1$.

Legendre Polynomials

${\begin{aligned}P_{n}(x)&={\frac {\mathrm {d} ^{n}}{\mathrm {d} x^{n}}}\left((1-x^{2})^{n}\right)\\&=(-1)^{n}{\frac {(2n)!}{n!}}\,\left(x-x_{1}\right)\,\left(x-x_{2}\right)\,\dots \,\left(x-x_{n}\right)\end{aligned}}$

Our goal is to find roots where $P_{n}(x)=0$

Recall that the inner product space $L^{2}[-1,1]$ is $\left\langle f,g\right\rangle =\int _{-1}^{1}f(x)\,g(x)\,\mathrm {d} x$

$\int _{-1}^{1}P_{n}(x)\,q(x)\,\mathrm {d} x=0$ for any $q(x)$, where $\mathrm {deg} {(q)}\leq n-1$.

Theorem

Theorem. [Gaussian Rule]. The degree of accuracy of the Gaussian Rule is $2n-1$ if and only if for all polynomials $p(x)$ of degree $\leq 2n-1$

$\int _{-1}^{1}P(x)\,\mathrm {d} x=\sum _{i=1}^{n}A_{i}\,p(x_{i})$

Proof. Take $p(x)$ such that the degree of $p$ is at most $n-1$. Then $p(x)$ is identical the Lagrange polynomial $\sum _{i=1}^{n}\sum _{i=1}^{n}p(x_{i})\,L_{n,i}(x)=\sum _{i=1}^{n}p(x_{i})\,\left(\prod _{j=1;j\neq i}^{n}{\frac {(x-x_{j})}{(x_{i}-x_{j})}}\right)$

Since the polynomials are equal, then their integrals are equal:

$\int _{-1}^{1}p(x)\,\mathrm {d} x=\sum _{i=1}^{n}p(x_{i})\,\int _{-1}^{1}L_{n,i}(x)\,\mathrm {d} x=\sum _{i=1}^{n}P(x_{i})\,A_{i}$

Now take $p(x)$ such that the degree of $p$ is at most $2n-1$. Then $p(x)=Q(x)\,P_{n}(x)+R(x)$, where $\mathrm {deg} {(Q)},\mathrm {deg} {(R)}\leq n-1$.

$\int _{-1}^{1}p(x)\,\mathrm {d} x={\cancel {\int _{-1}^{1}Q(x)\,P_{n}(x)\,\mathrm {d} x}}+\int _{-1}^{1}R(x)=\int _{-1}^{1}R(x)\,\mathrm {d} x=\sum _{i=1}^{n}A_{i}\,p(x_{i})$

Observe that the integral is equilavent to zero because $\mathrm {deg} {(Q)}\leq n-1$, and integrating by parts $n$ times gives ${\mbox{boundary terms}}+\int _{-1}^{1}\underbrace {{\frac {\mathrm {d} ^{n}Q}{\mathrm {d} x^{n}}}(x)} _{=0}\,(1-x^{2})^{n}$ Evaluating the boundary terms $(1-x^{2})^{n}\,(1+x^{2})^{n}$ at the points $-1$ and $1$ both give zero, so we are left with $0+0=0$.

quod erat demonstrandum

Now on an arbitrary interval, we have ${\frac {2}{b-a}}\int _{a}^{b}P(x)\,\mathrm {d} x=\sum _{i=1}^{n}A_{i}\,P(x_{i})$ simply by change of variables.

Example

$P_{5}(x)=-60x(120-560x^{2}+504x^{4})$ has 5 roots $x_{1}<x_{2}<x_{3}<x_{4}<x_{5}$. We know that $x_{3}=0$, $x_{2}=-x_{4}$, $x_{1}=-x_{5}$.

Solving the equation gives

Error Estimation

Remember Hermite polynomials?

$f(x)-H_{2n-1}(x)={\frac {f^{(2n)}(\xi (x))}{(2n)!}}\,(x-x_{1})^{2}\,(x-x_{2})^{2}\,\dots \,(x-x_{n})^{2}$

Our error is then

$E(x)=\int _{a}^{b}f(x)\,\mathrm {d} x-\int _{a}^{b}H_{2n-1}(x)\,\mathrm {d} x=\int _{a}^{b}{\frac {f^{(2n)}(\xi (x))}{(2n)!}}\,(x-x_{1})^{2}\,(x-x_{2})^{2}\,\dots \,(x-x_{n})^{2}$

By the intermediate value theorem,

$\int _{a}^{b}{\frac {f^{(2n)}(\xi (x))}{(2n)!}}\,(x-x_{1})^{2}\,(x-x_{2})^{2}\,\dots \,(x-x_{n})^{2}={\frac {f^{(2n)}(\zeta )}{(2n)!}}\int _{a}^{b}(x-x_{1})^{2}\,(x-x_{2})^{2}\,\dots \,(x-x_{n})^{2}$

If we let $h:=\max _{i}{(x-x_{i})}$, then

$\left|E(x)\right|\leq {\frac {\left|f^{(2n)}(\zeta )\right|}{(2n)!}}h^{2n+1}$

Composite Rule

Recall that the compasite simpson's rule error was

$-{\frac {f{(4)}(\zeta )}{90}}\,h^{4}$

In comparison, the composite gaussian rule error is at most

$\max _{a\leq x\leq b}\left|f^{(2n)}(x)\right|\,{\frac {1}{(2n)!}}\,h^{2n}$

Example

For example, if $\left|f^{(4)}\right|\leq 1$ and $\epsilon =10^{-10}$, then how many intervals $N$ do we need?

Gaussian

${\frac {c}{(2n)!}}\,\left({\frac {1}{N}}\right)^{2n}<10^{-10}$

If we use the $n=5$ rule as used in the previous example then $N=9$, so we sample the function at $n\cdot N=45$ points.

Simpson

${\frac {c}{90}}\,\left({\frac {1}{N}}\right)^{4}<10^{-10}$

$N\approx 10^{\frac {10}{4}}=10^{2.5}\approx 300$, so we sample the function at $2\cdot 300=600$ points

Midpoint

${\frac {c}{24}}\,\left({\frac {1}{N}}\right)^{2}<10^{-10}$

$N\approx 10^{5}$, so we sample the function at $10^{5}$ points

Trapezoidal rule

${\frac {c}{12}}\,\left({\frac {1}{N}}\right)^{2}<10^{-10}$

$N\approx 10^{5}$, so we sample the function at $10^{5}+1$ points

Hence the midpoint and Trapezoidal rules are lousy, but for a concave-up function (i.e. $f''>0$), the trapezoidal rule always overestimates, and the midpoint rule always underestimates. Therefore, we can get a bound on the value of the exact integral

Section 4.6: Adaptive Integration

If we have a function that "behaves nicely" (i.e. is flat or behaves like a cubic function) on certain intervals and does something interesting on other intervals, we can use several different rules on each interval, but adaptive integratino does this sort of thing "automatically"

Suppose we use Simpson's rule over an entire integral $\left[a,b\right]$ (let $h={\frac {b-a}{2}}$ be the midpoint between $a$ and $b$):

$\int _{a}^{b}-{\frac {b-a}{2}}\left(f(a)+4f(a+h)+f(b)=-{\frac {f^{(4)}(\mu )}{90}}h^{5}\right)$

Now what if we split the interval at the midpoint and compute two Simpson's rules:

$\int _{a}^{a+h}f(x)\,\mathrm {d} x+\int _{a+h}^{b}f(x)\,\mathrm {d} x-\left[S\left(a,{\frac {a+b}{2}}\right)+S\left({\frac {a+b}{2}},b\right)\right]=\ldots ={\frac {1}{16}}\left(-{\frac {f^{(4)}(\xi )}{90}}\,h^{5}\right)$

If we assume $f^{(4)}$ is either constant or does not change much, we can assume $f''(\mu )\approx f''(\xi )$:

$\Delta =S\left(a,{\frac {a+b}{2}}\right)+S\left({\frac {a+b}{2}},b\right)-S(a,b)\approx 15\cdot \left(-{\frac {1}{16}}\,{\frac {f^{(4)}(\xi )}{90}}h^{5}\right)=15\cdot {\mbox{error of}}\ S\left(a,{\frac {a+b}{2}}\right)+S\left({\frac {a+b}{2}},b\right)$

Algorithm 4.3

Given:

• $\epsilon >0$
• "flexibility coefficient" $k$ (usually 10)
• interval $[a,b]$

1. Compute $\Delta$ (do simpson's rule on $[a,b]$ and on half-intervals)
2. if $\left|\Delta \right|<k\,\epsilon$, then stop
3. Otherwise, subdivide $\left[a,b\right]$ into half-intervals with ${\frac {\epsilon }{2}}$ accuracy
4. Go to 1 (recursive)