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}$

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle N \approx 10^5} , so we sample the function at Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle 10^5} points

Trapezoidal rule

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \frac{c}{12} \, \left( \frac{1}{N} \right)^2 < 10^{-10}}

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle N \approx 10^5} , so we sample the function at Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle 10^5+1} points

Hence the midpoint and Trapezoidal rules are lousy, but for a concave-up function (i.e. Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle 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 Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \left[ a,b \right]} (let $h={\frac {b-a}{2}}$ be the midpoint between Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle a} and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle b} ):

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \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:

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \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 Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle f^{(4)}} is either constant or does not change much, we can assume Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle f''(\mu) \approx f''(\xi)} :

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \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:

• Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \epsilon > 0}
• "flexibility coefficient" Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle k} (usually 10)
• interval Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle [a,b]}

1. Compute Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \Delta} (do simpson's rule on Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle [a,b]} and on half-intervals)
2. if Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \left| \Delta \right| < k \, \epsilon} , then stop
3. Otherwise, subdivide Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \left[ a,b \right]} into half-intervals with Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \frac{\epsilon}{2}} accuracy
4. Go to 1 (recursive)