MATH 417 Lecture 10

« previous | Thursday, February 13, 2014 | next »

Section 4.3: Numerical Integration

If we have a function $f(x)$ that is impossible (or really hard) to integrate exactly on an interval $[a,b]$, we can estimate the integral based on the values of $f$ at $n+1$ points $x_0,x_1,\dots ,x_n$.

1. Choose a partition of the interval $[a,b]$ → ${\left\{I_j\right\}}_{j=1}^N$
2. Pick ${\xi }_j\in I_j$
3. Form the riemann sum $\sum _{j=1}^{N}f({\xi }_j)\left|I_j\right|$.

We usually define partition intervals by a diameter $h=\underset{j}{\max }\left|I_j\right|$. The Riemann integral, if it exists, is the limit of the sum as $h\to 0$.

In practice, we approximate $f(x)$ with a polynomial $p_n(x)$, and then integrate the polynomial:

$\underset{a}{\overset{b}{\int }}f(x)\mathrm{d}x=\underset{a}{\overset{b}{\int }}{p}_n(x)\mathrm{d}x+\underset{a}{\overset{b}{\int }}\frac{f^{(n+1)}(\xi (x))}{(n+1)!}(x-x_0)\dots (x-x_n)\mathrm{d}x$

Our error is going to be 0 if $f$ is a polynomial of degree $n$.

We can get an upper bound on this error by approximating $x-x_0\le b-a$. Hence

$|\underset{a}{\overset{b}{\int }}\frac{f^{(n+1)}(\xi (x))}{(n+1)!}(x-x_0)\dots (x-x_n)\mathrm{d}x|\le \frac{C}{(n+1)!}{|b-a|}^{n+1}$

For equidistant points, this bound can be better estimated by

$|\underset{a}{\overset{b}{\int }}\frac{f^{(n+1)}(\xi (x))}{(n+1)!}(x-x_0)\dots (x-x_n)\mathrm{d}x|\le \frac{C}{(n+1)!}{h}^n$

1. Construct a numerical rule which is exact for degree at most $n$
2. Find the degree of accuracy of the rule

Quadrature Rule

Given $x_0,x_1,\dots ,x_n$,

$\underset{a}{\overset{b}{\int }}f(x)\mathrm{d}x\approx \underset{a}{\overset{b}{\int }}{p}_n(x)\mathrm{d}x=\underset{a}{\overset{b}{\int }}\sum _{j=0}^{n}f(x_j)\prod _{i=0;i\ne j}^n\frac{x-x_i}{x_j-x_i}\mathrm{d}x=\sum _{j=0}^{n}f(x_j)\underset{a}{\overset{b}{\int }}{p}_{n,j}(x)\mathrm{d}x$

Thus

$\underset{a}{\overset{b}{\int }}f(x)\mathrm{d}x\approx Q_n(f)=\sum _{j=0}^n{A}_{j}f(x_j)$, where $A_j=\underset{a}{\overset{b}{\int }}{p}_{n,j}(x)\mathrm{d}x=\underset{a}{\overset{b}{\int }}\frac{(x-x_0)\dots (x-x_n)}{(x_j-x_0)\dots (x_j-x_n)}\mathrm{d}x$

In a simple case where $x_0=a$ and $x_1=b$, we get the trapezoid rule.

The degree of accuracy is between $n$ and $2n+1$.

To find $A_i$'s, we can set up a linear system instead of integrating the Lagrange polynomials

$\begin{array}{rl}{b}_0=\underset{a}{\overset{b}{\int }}1\mathrm{d}x& =A_0\cdot 1+A_1\cdot 1+\dots +A_n\cdot 1\\ b_1=\underset{a}{\overset{b}{\int }}x\mathrm{d}x& =A_0\cdot x_0+A_1\cdot x_1+\dots +A_n\cdot x_n\\ b_2=\underset{a}{\overset{b}{\int }}{x}^2\mathrm{d}x& =A_0\cdot x_0^2+A_1\cdot x_1^2+\dots +A_n\cdot x_n^2\\ ⋮& =⋮\\ b_n=\underset{a}{\overset{b}{\int }}{x}^n\mathrm{d}x& =A_0\cdot x_0^n+A_1\cdot x_1^n+\dots A_n\cdot x_n^n\end{array}$

Change of Variables

Suppose we have a rule for $f(x)$ when $x\in [-1,1]$. We can scale the rule to fit an arbitrary interval $y\in [a,b]$ as follows:

$-1\le x\le 1$

$a\le y(x)\le b$

$y(x)=\alpha x+\beta =\frac{b-a}{2}x+\frac{a+b}{2}$

$f(x)=g(y(x))$, where $g$ is a polynomial of the same degree:

Performing a substitution in the integral $\underset{-1}{\overset{1}{\int }}g(x)\mathrm{d}x$ gives $\underset{-1}{\overset{1}{\int }}g(y(x))y^{\prime }(x)\mathrm{d}x$

$y^{\prime }(x)=\frac{b-a}{2}$

This integral can now be approximated by $B_{0}g(y(-1))+B_{1}g(y(0))+B_{2}g(y(1))$, where

$B_j=\frac{b-a}{2}{A}_j$

$\begin{array}{rl}f(x)& =1+x+x^2\\ g(y)& =1+\frac{y+1}{2}+{\left(\frac{y+1}{2}\right)}^2\end{array}$

If $\underset{a}{\overset{b}{\int }}f(x)\mathrm{d}x=\sum _{j=0}^n{A}_{j}f(x_j)$ for $f(x)=1,x,x^2,\dots ,x^4$, then it is also exact for $p(x)$ of degree $\le n$.

Take $\underset{a}{\overset{b}{\int }}f(x)\mathrm{d}x=\underset{a}{\overset{b}{\int }}\sum _{k=0}^n{c}_k{x}^k=\sum _{k=0}^n{c}_k\underset{a}{\overset{b}{\int }}{x}^k\mathrm{d}x=\sum _{k=0}^n{c}_k\sum _{j=0}^n{A}_j{x}_j^k=\sum _{k=0}^n\sum _{j=0}^n{A}_j{c}_k{x}_j^k=\sum _{j=0}^n{A}_j\sum _{k=0}^n{c}_k{x}_j^k=\sum _{j=0}^n{A}_{j}f(x_j)$