Numerical Differentiation and Integration

Numerical Differentiation

Problem Defined

In many practical applications, we may want the derivative of a function $f$, without evaluating the true expression for $f^{\prime }(x)$. This may occur because:

• The true expression for $f^{\prime }(x)$ is extremely complicated, or expensive to compute!
• We only have points for $f$ - $\{(x_i,f(x_i))\}$, so we don’t even have an expression for $f$ to compute a derivative with!

How can we find such derivatives?

Derivations Using Taylor Approximations

One way we can derive an iterative solution for derivatives is by using their Taylor Approximation!

For example, recall that

$$f^{\prime }(x)=\underset{h\to 0}{\lim }\frac{f(x+h)-f(x)}{h}$$

If $h>0$ and is very small, we approximately have that

$$f^{\prime }(x)\approx \frac{f(x+h)-f(x)}{h}$$

By the Taylor Remainder Theorem of $f(x+h)$ about $x$, we have that for some $z\in (x,x+h)$,

$$\begin{array}{rl}f(x+h)& =f(x)+f^{\prime }(x)\cdot h+\frac{f^{\prime \prime }(x)}{2!}{h}^2+\frac{f^{\prime \prime \prime }(x)}{3!}{h}^3+\frac{f^{\prime \prime \prime }(z)}{4!}{h}^4\\ f^{\prime }(x)& =\frac{f(x+h)-f(x)}{h}-\frac{f^{\prime \prime }(x)}{2!}h-\frac{f^{\prime \prime \prime }(x)}{3!}{h}^2-\frac{f^{\prime \prime \prime }(z)}{4!}{h}^3\\ f^{\prime }(x)& =\frac{f(x+h)-f(x)}{h}+\text{Er(h)}\end{array}$$

Where $\text{Er(h)}$ is the error value on $h$, which converges to 0 as h becomes small. This has a $1^{st}$ order error, and approximates our derivative from the left-hand side, known as the forward difference!

Orders of Error

For any error term $\text{Er(h)}$, we say that it is a $p^{th}$ order method ($O(h^P)$), where $p$ is the leading power term (lowest exponent) of $h$ in $\text{Er(h)}$, as the lowest exponent will contribute the most to the error.

The higher the order of error, the better!

We similarly have a backward difference,

$$f^{\prime }(x)=\frac{f(x)-f(x-h)}{h}+\text{Err(h)}$$

which essentially approximates our derivative from the right-hand side (with $1^{st}$ order of error), and a central difference,

$$f^{\prime }(x)=\frac{f(x+h)-f(x-h)}{2h}+\text{Er}(h)$$

which approximates our derivative using both sides. We can show that this error term has a $2^{nd}$ order of error, making it better than the previous two!

Central Difference for $f^{\prime \prime }$

We also can find a central difference for $f^{\prime \prime }$ as

$$f^{\prime \prime }(x)=\frac{f(x+h)-2f(x)+f(x-h)}{h^2}+\text{Er}(h)$$

Which similarly has an error term for $O(h^2)$.

In general, we have for any $k\in \mathbb{R}$, that there exists $z\in (x,x+kh)$

$$f(x+kh)=f(x)+f^{\prime }(x)\frac{(kh)}{1!}+f^{\prime \prime }(x)\frac{(kh{)}^2}{2!}+f^{\prime \prime \prime }(x)\frac{(kh{)}^3}{3!}+f^{\prime \prime \prime \prime }(z)\frac{(kh{)}^4}{4!}$$

We can use combinations of these Taylor expansions to approximate our derivatives!

Derivations Using Undetermined Coefficients

Sometimes, simple Taylor expansions are not enough - in these cases, we may still be able to approximate our derivative using the method of undetermined coefficients! In such a problem, we are given the following values

$$\{f(x+k_{i}h){\}}_{i=1}^n$$

where $k_i\in \mathbb{R}$, $i\ne j\to k_i\ne k_j$, and want to find a linear combination of these values to approximate some derivative $f^r(x)$. We then use their Taylor expansions to solve for the coefficients.

To see how this is done, consider the following example.

Say we want to approximate $f^{\prime \prime }(x)$ using a linear combination of $f(x-h),f(x),f(x+h)$. Then, we can approximate it as

$$f^{\prime \prime }(x)=c_{1}f(x-h)+c_{2}f(x)+c_{3}f(x+h)=L(h)$$

We expand each using their Taylor Expansions (performed above) to obtain

$$\begin{array}{rl}L(h)& =c_1\left[f(x)-f^{\prime }(x)h+\frac{f^{\prime \prime }(x)}{2!}{h}^2-\frac{f^{\prime \prime \prime }(x)}{3!}{h}^3+\frac{f^{\prime \prime \prime \prime }(z_1)}{4!}{h}^4\right]\\ & +c_{2}f(x)\\ & +c_3\left[f(x)+f^{\prime }(x)h+\frac{f^{\prime \prime }(x)}{2!}{h}^2+\frac{f^{\prime \prime \prime }(x)}{3!}{h}^3+\frac{f^{\prime \prime \prime \prime }(z_2)}{4!}{h}^4\right]\\ & =(c_1+c_2+c_3)f(x)+(-c_1+c_3)h{f}^{\prime }(x)+(c_1+c_3)\frac{h^2}{2}{f}^{\prime \prime }(x)\\ & +(-c_1+c_3)\frac{h^3}{3!}{f}^{\prime \prime \prime }(x)+(c_1{f}^{\prime \prime \prime \prime }(z_1)+c_3{f}^{\prime \prime \prime \prime }(z_2))\frac{h^4}{4!}\end{array}$$

We now attempt to find $c_1,c_2,c_3$ necessary to find our answer. As we want our second derivative, we want to drop the coefficients for our $f$ and $f^{\prime }$ terms to 0, but make the coefficient for the $f^{\prime \prime }$ term 1!

$$\begin{array}{rl}{c}_1+c_2+c_3& =0\\ -c_{1}h+c_{3}h& =0\\ c_1\frac{h^2}{2}+c_3\frac{h^2}{2}& =1\end{array}$$

We solve this as

$$c_1=c_3=\frac{1}{h^2}{c}_2=-\frac{2}{h^2}$$

Subbing this back into our equation, we get

$$f^{\prime \prime }(x)=\frac{f(x-h)-2f(x)+f(x+h)}{h^2}+O(h^2)$$

Which we can use to approximate our $f^{\prime \prime }(x)$ value!

Rates of Convergence

Given an approximation $L(h)$ of $T$, the error is given as

$$E(h)=|T-L(h)|$$

Suppose the error is approximately equal to $C{h}^p$, where $C$ is a constant, indicating that it is a $p$ order method.

If we take the logarithm of this error, we obtain

$$\log (E(h))\approx \log C+p\log h$$

We can use this fact to to find our order of error (rate of convergence)! So, given our error on $h_1,h_2$,

$$E(h_1)=C{h}_1^{p}E(h_2)=C{h}_2^p$$

We can find our order of error as

$$\begin{array}{r}\frac{E(h_1)}{E(h_2)}={\left(\frac{h_1}{h_2}\right)}^p\\ \log \frac{E(h_1)}{E(h_2)}=\log {\left(\frac{h_1}{h_2}\right)}^p\\ p=\frac{\log E(h_1)/\log E(h_2)}{\log h_1/\log h_2}\end{array}$$

Richardson’s Extrapolation

Richardson’s Extrapolation describes a general method to boost the accuracy of our approximations.

Consider the $p^{th}$ order approximation

$$T=L(h)+e_0{h}^p+e_1{h}^{p+1}+\dots$$

Where $T$ is our target, $L(h)$ is our algorithm, and every $e_i{h}^{p+i}$ is part of the error term.

Note that if we wanted to change the order of our approximation, we’d need to eliminate the $e_0{h}^p$ term.

Consider what happens if we replace $h$ with $2h$.

$$T=L(2h)+e_0{2}^p{h}^p+e_1{2}^{p+1}{h}^{p+1}+\dots$$

Then, we can use this equation to elminate our $e_0{h}^p$ error term. Multiplying our first equation by $2^p$ and subtracting our second equation from it, we get

$$\begin{array}{rl}(2^p-1)T& =2^{p}L(h)-L(2h)\\ & +2^p(e_0{h}^p+e_1{h}^{p+1}+\dots )-2^p(e_0{h}^p+2e_1{h}^{p+1}+4e_2{h}^{p+2}+\dots )\\ T& =\frac{2^{p}L(h)-L(2h)}{2^p-1}+{\hat{e}}_1{h}^{p+1}+{\hat{e}}_2{h}^{p+2}+\dots \end{array}$$

Giving us a new algorithm that increases our order of error!

Numerical Integration

In other practical applications, we may want the integral of a function $f$, ${\int }_a^{b}f(x)dx$. This may occur because

• We may not know $f(x)$ in order to take the integral.
• We may know $f(x)$, but a closed form expression of the integral may not be available.
• We may be able to compute the integral, but it is difficult to compute.
• We may only have $f$ at a few nodes!

To do these approximations, we’ll perform what are known as quadrature rules, rules that let us approximate integrals.

Integration via Interpolation

Possibly one of the more intuitive approaches to approximations is to first interpolate $f$ given its points, and then integrate this interpolation.

Assume we’re given points $\{(x_i,f(x_i){\}}_{i=0}^n$, where $\forall 0\le i\le n$, $x_i\in [a,b]$. Recall that for these points, we can find our Lagrange Polynomial as

$$f(x)\approx P_n(x)=\sum _{i=0}^{n}f(x_i)l_i(x)$$

Where $l_i(x)$ are the bais polynomials. Then, as this approximates $f$, we may also be able to approximate its integral using $P_n(x)$!

$${\int }_a^{b}f(x)dx\approx {\int }_a^b{P}_n(x)dx=\sum _{i=0}^{n}f(x_i){\int }_a^b{l}_i(x)dx$$

This can alternatively be given as

$${\int }_a^{b}f(x)dx\approx \sum _{i=0}^n{A}_{i}f(x_i)A_i={\int }_a^b{l}_i(x)dx$$

Similar to our Lagrange Interpolation, this has an advantage in that we can save on a lot of computation for different functions on the same nodes!

Note that this is essentially taking a “weighted sum”, as the integrals define a value determining how much each $f(x_i)$ contributes to the result.

We discuss various ways we can use this below!

Closed Newton-Cotes Quadrature

Suppose we have our equipspaced interpolation nodes on $[a,b]$, given by

$$x_i=a+ih\forall 0\le i\le n$$

Where $h=\frac{b-a}{n}$. Note that by this, $x_0=a$, and $x_n=b$.

Then, using our integral approximation from before, we’ll obtain approximation

$${\int }_a^{b}f(x)dx=\sum _{i=0}^n{A}_{i}f(x_i)+Er(h)$$

Theorem: Error

There exists a $ϵ\in (a,b)$ such that

$$Er(h)=\{\begin{array}{ll}\frac{h^{n+3}{f}^{n+2}(ϵ)}{(n+2)!}{\int }_0^n{t}^2(t-1)(t-2)\dots (t-n)dt& \text{n even and }f\in C^{n+2}[a,b]\\ \frac{h^{n+2}{f}^{n+1}(ϵ)}{(n+1)!}{\int }_0^{n}t(t-1)\dots (t-n)dt& \text{n odd and }f\in C^{n+1}[a,b]\end{array}$$

Note that this error is 0, meaning this quadrature is exact if $f(x)$ is a polynomial such that

• It’s degree $\mathrm{deg}(f)\le n+1$ if $n$ is even
• It’s degree $\mathrm{deg}(f)\le n$ if $n$ is odd

As in these cases, the derivative in the error term drops to 0.

Some important quadratures from this rule are given below.

• When $n=1$, we have the Trapezoidal rule, giving us $$\begin{array}{rlr}& x_0=a& x_1=b\\ & {\int }_a^{b}f(x)dx=\left(\frac{f(a)+f(b)}{2}\right)(b-a)-\frac{h^3}{12}{f}^{\prime \prime }(ϵ)& \end{array}$$
• When $n=2$, we have the Simpson’s rule, giving us $$\begin{array}{rl}& x_0=a{x}_1=\frac{a+b}{2}{x}_2=b\\ & {\int }_a^{b}f(x)dx=\frac{f(a)+4f(\frac{a+b}{2})+f(b)}{6}(b-a)-\frac{h^5}{90}{f}^4(ϵ)\end{array}$$

Open Newton-Cotes Quadrature

Suppose we have a set of nodes $\{x_i{\}}_{i=0}^n\subset (a,b)$, such that $\forall 0\le i\le n$,

$$x_i=a+(i+1)h$$

Where $h=\frac{b-a}{n+2}$. Note that by this definition, $x_0=a+h$, and $x_n=b-h$.

Then, for some $ϵ\in (a,b)$, we obtain

$${\int }_a^{b}f(x)dx=\{\begin{array}{ll}{\sum }_{i=0}^n{A}_{i}f(x_i)+\frac{h^{n+3}{f}^{n+2}(ϵ)}{(n+2)!}{\int }_{-1}^{n+1}{t}^2(t-1)(t-2)\dots (t-n)dt& \text{n even}\\ {\sum }_{i=0}^n{A}_{i}f(x_i)+\frac{h^{n+2}{f}^{n+1}(ϵ)}{(n+1)!}{\int }_{-1}^{n+1}{t}^2(t-1)(t-2)\dots (t-n)dt& \text{n odd}\end{array}$$

Note that we have exactness if $f(x)$ is a polynomial under the same conditions as the closed-form algorithm.

Some important quadratures from this rule are given below.

• When $n=0$, we have the Midpoint rule, giving us $$\begin{array}{rl}& x_0=\frac{a+b}{2}\\ & {\int }_a^{b}f(x)dx=f\left(\frac{a+b}{2}\right)(b-a)+\frac{h^3}{3}{f}^{\prime \prime }(ϵ)\end{array}$$

  This takes the integral of a constant line, whose value is equal to $f$‘s midpoint.

Composite Quadrature Rules

In each of these quadrature rules, because our errors depend on $h$, as $n\to \infty$, $h\to 0$, giving us lower errors! So, we should interpolate with higher-degree polynomials, right?

Wrong! As we use higher-degree polynomials, we run into another issue - Runge’s Phenomena! So, to minimize our errors, we should use other methods (as described below).

One such method is by applying a composite quadrature rule, which partitions $[a,b]$ into smaller intervals and applies lower-order quadrature rules on each. So, given $\{x_i{\}}_{j=0}^m$, we’ll have partition

$$x_0<x_1<x_2<\cdots <x_{m-1}<x_m$$

And can apply our quadrature rules on each $[x_j,x_{j+1}]$ subinterval with size $x_{j+1}-x_j$, whose composite will form our function integral.

$${\int }_a^{b}f(x)dx={\int }_{x_0}^{x_1}f(x)dx+{\int }_{x_1}^{x_2}f(x)dx+\cdots +{\int }_{x_{m-1}}^{x_m}f(x)dx$$

Each subquadrature is given as

$${\int }_{x_j}^{x_{j+1}}f(x)dx\approx \sum _{i=0}^{n_j}{A}_i^{j}f(x_i^j)$$

Where $A_i$ are the quadrature weights on $[x_j,x_{j+1}]$, and $\{x_i^j{\}}_{i=0}^{n_j}$ are the quadrature nodes in our subinterval. For each subquadrature, we commonly will apply our midpoint rule, trapezoidal rule, or simpson’s rule.

Suppose we applied our midpoint rule. We can prove that our new error is of order 2!

Proof

For simplicity, assume that the size of our subintervals is constant. Then, if we applied our midpoint rule, for any subinterval, we get error term

$$Err=\frac{h^3}{3}{f}^{\prime \prime }(ϵ)$$

For some $ϵ$ in our subinterval. Then, our total composite error is given as

$$Err=\sum _{j=0}^{m-1}E{r}_j=\frac{h^3}{3}\sum _{j=0}^{m-1}{f}^{\prime \prime }({ϵ}_j)$$

Defining $M=\frac{1}{m}{\sum }_{j=0}^{m-1}{f}^{\prime \prime }({ϵ}_j)$, the average of these node values, we know that this average is bounded by $f^{\prime \prime }(x)$‘s maximum and minimum along the entire interval! So, we can apply intermediate value theorem to find a $c$ such that $f^{\prime \prime }(c)=M$.

This gives us error

$$Er=\frac{m{h}^3}{3}{f}^{\prime \prime }(c)$$

For some $c$ in our interval $[a,b]$. Applying our assumption that $h=\frac{(b-a)}{m}$, we get

$$Er=(b-a)\frac{h^2}{3}{f}^{\prime \prime }(c)$$

Telling us that our error term is a 2nd order method!

We can similary find that the composition trapezoidal rule is a 2nd order method, and our composition simpson’s rule is a 4th order method.

Romberg Integration

Similarly to how we have Richardson’s Extrapolation for boosting accuracy, we can also apply something known as Romberg Integration to boost the accuracy of composite integration.

Applying a similar process to Richardson’s Extrapolation, we obtain

$$T=\frac{2^{p}Q(h)-Q(2h)}{2^p-1}+\widehat{e_1}{h}^{p+1}+\widehat{e_2}{h}^{p+2}$$

Giving us a new quadrature which is of a higher order error!

However, we need to be a bit careful with this - replacing $h$ by $2h$ means we’re halving the number of subintervals we have for this second quadrature! This is only possible if $m$ is even, as otherwise, we won’t be able to evenly cover our entire domain.

Doing this also means we have to evaluate $f$ at different points, which are actually where the partition boundaries once were!

Gaussian Quadratures

Note that in the prior quadrature rules, we formed approximations of the form

$${\int }_a^{b}f(x)dx\approx \sum _{i=0}^n{A}_{i}f(x_i)$$

Where our quadrature nodes are given.

How do we pick the quadrature nodes and quadrature weights such that our approximate is exact for the maximal degre of polynomials?

Example: Exact Quadratures

The quadrature rule

$${\int }_{-1}^{1}f(x)dx\approx f\left(-\frac{1}{\sqrt{3}}\right)+f\left(\frac{1}{\sqrt{3}}\right)$$

Is exact on $P_3[-1,1]$.

So, suppose we have $(n+1)$ quadrature nodes and $(n+1)$ quadrature weights. This gives us $(2n+2)$ unknowns! Thus, it may be possible to exactly represent polynomials in ${\prod }_{2n+1}$.

Note that any polynomial $p\in {\prod }_n[a,b]$ can be expressed as

$$p(x)=\sum _{i=0}^{n}p(x_i)l_i(x)l_i=\prod _{j=0}^n\frac{(x-x_j)}{(x_i-x_j)}j\ne i,0\le i\le n$$

So, our rule is exact on ${\prod }_n[a,b]$ if it is equal to the integral of this!

$$\begin{array}{r}\sum _{i=0}^n{A}_{i}p(x_i)={\int }_a^{b}p(x)dx=\sum _{i=0}^{n}p(x_i){\int }_a^b{l}_i(x)dx\\ ⟹A_i={\int }_a^b{l}_i(x)dx\end{array}$$

Below, we state a theorem to help us determine what nodes and weights we should be choosing.

Theorem: Properties of Orthogonal Functions

If $f$ is a non-zero function on $[a,b]$ such that

1. $f\in C[a,b]$
2. $f$ is orthogonal to ${\prod }_n[a,b]$, or in other words $${\int }_a^{b}f(x)p(x)dx=0p\in \prod _n[a,b]$$

Then, $f$ will change sign at least $(n+1)$ times in $[a,b]$ - it has at least $n$ real distinct zeroes in $[a,b]$.

Theorem: Exactness of Quadrature Rules

Let $q$ be a polynomial satisfying the following properties:

1. $q(x)\not\equiv 0$ and is of the degree $(n+1)$.
2. $q(x)$ is orthogonal to ${\prod }_n[a,b]$.
3. $q(x)$ has real, simple roots, and all roots are in $[a,b]$.

Then, the quadrature rule ${\sum }_{i=0}^n{A}_{i}f(x_i)$ with these roots as the quadrature nodes and the quadrature weights

$$A_i={\int }_a^b{l}_i(x)l_i=\prod _{j=0}^n\frac{(x-x_j)}{(x_i-x_j)}j\ne i$$

Is exact on ${\prod }_{2n+1}[a,b]$.

Proof

Let $p(x)={\prod }_{2n+1}[a,b]$. By property of polynomials, we can factorize $p(x)$ as

$$p(x)=q(x)\cdot \tilde{p}(x)+r(x)$$

Where $q(x)$ is a polynomial of most $(n+1)$ satisfying the above theorem, and $\tilde{p}(x)$ and $r(x)$ are polynomials of at most $n$.

At the roots of $q(x)$ (which are guaranteed by the above theorem), we have

$$q(x_i)=0⟹p(x_i)=r(x_i)$$

Taking the integral of $p(x)$, we have

$$\begin{array}{rl}{\int }_a^{b}p(x)dx& ={\int }_a^{b}q(x)\tilde{p}(x)dx+{\int }_a^{b}r(x)dx\\ & =0+{\int }_a^{b}r(x)\\ & =\sum _{i=0}^{n}r(x_i){\int }_a^b{l}_i(x)dx\\ & =\sum _{i=0}^{n}p(x_i){\int }_a^b{l}_i(x)dx\end{array}$$

So, if we can find such a polynomial $q(x)$, we can find exactness!

How can we find such a $q(x)$?

Recall the Legendre Polynomials on $[-1,1]$, which are orthogonal

$${\varphi }_0(x)=1{\varphi }_1(x)=x{\varphi }_2(x)=\frac{3x^2-1}{2}$$

These polynomials have the following properties:

1. The degree of each polynomial ${\varphi }_k(x)$ is $k$.
2. ${\varphi }_{k+1}(x)$ is orthogonal to ${\prod }_k[-1,1]$.
3. By an earlier theorem, ${\varphi }_k(x)$ has $k$ distinct, real roots on $[-1,1]$.

So, on $[-1,1]$, we can choose $q(x)$ to be the Legendre polynomial on ${\varphi }_{n+1}(x)$!

We can, in fact, find our earlier example using ${\varphi }_2(x)$!

Example: Exact Quadratures

Recall

$${\int }_{-1}^{1}f(x)dx\approx f\left(-\frac{1}{\sqrt{3}}\right)+f\left(\frac{1}{\sqrt{3}}\right)$$

This is exact on ${\prod }_3[-1,1]$, so $2n+1=3⟹n=1$, so we choose $q(x)={\varphi }_{n+1}(x)={\varphi }_2(x)$.

We have roots at $x_0=-1/\sqrt{3}$ and $x_1=1/\sqrt{3}$. We choose these to be our quadrature nodes. We then have

$$\begin{array}{r}{A}_0={\int }_{-1}^1{l}_0(x)dx={\int }_{-1}^1\frac{x-x_1}{x_0-x_1}dx=1\\ A_1={\int }_{-1}^1{l}_1(x)dx={\int }_{-1}^1\frac{x-x_0}{x_1-x_0}dx=1\end{array}$$

Thus, our exact quadrature rule is $A_{0}f(x_0)+A_{1}f(x_1)$, giving us the above!

What if our interval isn’t $[-1,1]$? Then, we can apply a change of variables!

If $x\in [a,b]$, then we can set

$$t=\frac{2x-a-b}{b-a}⟺x=\frac{1}{2}[(b-a)t+a+b]$$

transforming our domain into the interval $[-1,1]$! This gives us integral

$$\begin{array}{rl}{\int }_a^{b}f(x)dx& ={\int }_{-1}^{1}f\left(\frac{1}{2}[(b-a)t+a+b]\right)\frac{dx}{dt}dt\\ & ={\int }_{-1}^{1}f\left(\frac{1}{2}[(b-a)t+a+b]\right)\frac{(b-a)}{2}dt\\ & ={\int }_{-1}^1\tilde{f}(t)dt\\ \tilde{f}(t)& =\frac{(b-a)}{2}f\left(\frac{1}{2}[(b-a)t+a+b]\right)\end{array}$$

We can use this to find the roots $t_i$ of ${\varphi }_{n+1}(t)$ in $(-1,1)$, and convert them to the quadrature nodes $x_i$!