Wichita State University Logo

Math 451: Computational Mathematics

3.7 Differentiation

3.7.1

One way to approximate the derivative of a univariant function , that is a function of one variable, at a point $x$ in the domain of $f$ is the definition of derivative and the difference quotient

$$ \begin{align*} f'(x) &= \displaystyle\lim_{h\rightarrow 0} \dfrac{f(x+h)-f(x)}{h} \\ \\ f'(x) &\approx \dfrac{f(x+h)-f(x)}{h} \end{align*} $$
where $h$ is a small number. If $h>0$ then the forward difference approximation is given by

$$ f'(x) \approx \dfrac{f(x+h)-f(x)}{h} $$
Forward Difference Approximation Figure 1


The backward difference approximation is given by

$$ f'(x) \approx \dfrac{f(x)-f(x-h)}{h} $$
Backward Difference Approximation Figure 2
The centered difference approximation is given by

$$ f'(x) \approx \dfrac{f(x+h)-f(x-h)}{2h} $$
Centered Difference Approximation Figure 3
The centered difference formula uses Mean Value Theorem from single variable, differential calculus.

The Mean Value Theorem

If $f:[a,b]\rightarrow\mathbb{R}$ is continuous on the closed interval $[a,b]$ and differentiable on the interior $(a,b)$, then there exists at least one point $c\in(a,b)$ so that

$$ f'(c) = \dfrac{f(b)-f(a)}{b-a} $$

In our example $c \approx \frac{a+b}{2} = \frac{x_0-h + x_0+h}{2} = x0$.

3.7.2 Accuracy of the Forward and Backward Difference Methods

If our function has at least two continuous derivatives ($f\in C^2[a,b]$), then we can determine the accuracy of the difference methods using Taylor's Theorem. For the forward difference method we have

$$ f(x+h) = f(x) + hf'(x) + \frac{h^2}{2}\,f''(\xi),\qquad\xi\in[x,x+h] $$
Solving for $f'(x)$ one obtains

$$ f'(x) = \dfrac{f(x+h)-f(x)}{h} - \frac{h}{2}\,f''(\xi) $$
So we have a truncation error or discretization error of

$$ \left|\,f'(x)-\dfrac{f(x+h)-f(x)}{h}\,\right| = \frac{h}{2}\,f''(\xi) $$
Similarly

$$ \begin{align*} f(x-h) &= f(x) - hf'(x) + \frac{h^2}{2}\,f''(\xi),\qquad\xi\in[x-h,x] \\ \\ f'(x) &= \dfrac{f(x)-f(x-h)}{h} + \frac{h}{2}\,f''(\xi) \\ \\ &\left|\,f'(x) - \dfrac{f(x)-f(x-h)}{h}\,\right| = \frac{h}{2}\,f''(\xi) \end{align*} $$
In either case the discretization error is $\frac{h}{2}\,f''(\xi)$, or $O(h)$.

However roundoff error also plays contributes to the finite difference quotient. For forward difference

$$ \dfrac{f(x+h)(1+\delta_1)-f(x)(1+\delta_2)}{h} = \dfrac{f(x+h)-f(x)}{h} + \dfrac{\delta_1f(x+h)-\delta_2f(x)}{h} $$



$$ \begin{align*} \left|\,\dfrac{f(x+h)(1+\delta_1)-f(x)(1+\delta_2)}{h} - \dfrac{f(x+h)-f(x)}{h}\,\right| &= \left|\,\dfrac{\delta_1f(x+h)-\delta_2f(x)}{h}\,\right| \\ \\ &\le \dfrac{\epsilon\left(|f(x)| + |f(x+h)| \right)}{h} \\ \\ &= O\left(\frac{\epsilon}{h}\right) \end{align*} $$
since both $\delta_1$ and $\delta_2$ are less that machine precision $\epsilon$. If we use a value for $h$ that is too small, the roundoff error will be large. If we us a value for $h$ that is too large, the truncation error will be large. The best accuracy we can hope for is for the two values for error are approximately equal.

$$ \begin{align*} h &\approx \frac{\epsilon}{h} \\ h^2 &\approx \epsilon \\ h &\approx \sqrt{\epsilon} \end{align*} $$
Thus the truncation error becomes $O(h) = O\left(\sqrt{\epsilon}\right) \approx \sqrt{\epsilon}$, and the roundoff error $O\left(\frac{\epsilon}{h}\right) = O\left(\frac{\epsilon}{\sqrt{\epsilon}}\right) \approx \sqrt{\epsilon}$.

Therefore the best accuracy we can hope for using forward or backward difference is $\sqrt{\epsilon}$ using a step size of $h = \sqrt{\epsilon}$.

3.7.3 Accuracy of the Centered Difference Method

If our function $f$ has at least three continuous derivative $f\in C^3[a,b]$, then we again use Taylor's formula

$$ \begin{align*} f(x+h) &= f(x) + hf'(x) + \frac{h^2}{2}\,f''(x) + \frac{h^3}{6}\,f'''(\xi),\qquad\xi\in[x,x+h] \\ \\ f(x-h) &= f(x) - hf'(x) + \frac{h^2}{2}\,f''(x) - \frac{h^3}{6}\,f'''(\eta),\qquad\eta\in[x,x+h] \\ \\ f(x+h)-f(x-h) &= 2hf'(x) + \frac{h^3}{6}\left(\,f'''(\xi)+f'''(\eta)\,\right) \\ \\ f'(x) &= \dfrac{f(x+h)-f(x-h)}{2h} - \frac{h^2}{12}\left(\,f'''(\xi)+f'''(\eta)\,\right) \\ \end{align*} $$
Our discretization error becomes

$$ \left|\,f'(x)-\dfrac{f(x+h)-f(x-h)}{2h}\,\right| = \frac{h^2}{12}\left(\,f'''(\xi)+f'''(\eta)\,\right) = O(h^2) $$
Balancing our truncation error with roundoff error yields

$$ \begin{align*} h^2 &\approx \frac{\epsilon}{h} \\ h^3 &\approx \epsilon \\ h &\approx \sqrt[3]{\epsilon} \end{align*} $$
Hence the truncation error becomes $O(h^2) \approx \epsilon^{2/3}$, and the roundoff error $O\left(\frac{\epsilon}{h}\right) = O\left(\frac{\epsilon}{\sqrt[3]{\epsilon}}\right) \approx \epsilon^{2/3}$.

With the centered difference we attain an accuracy of approximately $\epsilon^{2/3}$.

In [2]: 
format short
eps^(2/3)

ans =

   3.6669e-11

3.7.4 Richardson Extrapolation

Suppose that our function has four continuous derivatives so that $f\in C^4[a,b]$. The we can use Taylor's Theorem to derive

$$ \begin{align*} f(x+h) &= f(x) + hf'(x) + \frac{h^2}{2}\,f''(x) + \frac{h^3}{6}\,f'''(x) + \frac{h^4}{24}f^{(4)}(x) + \frac{h^5}{120}f^{(5)}(\xi) \\ \\ f(x-h) &= f(x) - hf'(x) + \frac{h^2}{2}\,f''(x) - \frac{h^3}{6}\,f'''(x) + \frac{h^4}{24}f^{(4)}(x) + \frac{h^5}{120}f^{(5)}(\eta) \end{align*} $$
where $\xi\in[x,x+h]$ and $\eta\in[x-h,x]$. Giving us

$$ f(x+h)-f(x-h) = 2h\,f'(x) + \frac{h^3}{3}\,f'''(x) + O\left(h^5\right) $$
Solving for $f'(x)$ yields

$$ f'(x) = \dfrac{f(x+h)-f(x-h)}{2h} - \frac{h^2}{6}\,f'''(x) + O\left(h^4\right) $$
We can define a new function

$$ f'(x) := \phi(h) = \dfrac{f(x+h)-f(x-h)}{2h} - \frac{h^2}{6}\,f'''(x) + O\left(h^4\right) $$
and compute our function $\phi$ at $\frac{h}{2}$, or compute $f$ at two more points $x-\frac{h}{2}$ and $x+\frac{h}{2}$.

$$ \begin{align*} f'(x) = \phi\left(\frac{h}{2}\right) &= \dfrac{f\left(x+h/2\right)-f\left(x-h/2\right)}{2(h/2)} - \frac{(h/2)^2}{6}\,f'''(x) + O\left(h^4\right) \\ \\ &= \dfrac{f\left(x+h/2\right)-f\left(x-h/2\right)}{h} - \frac{h^2}{24}\,f'''(x) + O\left(h^4\right) \end{align*} $$
The idea here is to combine the expressions for $\phi(h)$ and $\phi\left(\frac{h}{2}\right)$ so that the $f'''(x)$ terms add to zero. We need to values $A$ and $B$ so that

$$ A\phi\left(\frac{h}{2}\right) + B\phi(h) = A\,\dfrac{f\left(x+h/2\right)-f\left(x-h/2\right)}{2h} + B\,\dfrac{f(x+h)-f(x-h)}{2h} + O\left(h^4\right) $$
Thus

$$ \begin{align*} A\left(- \frac{h^2}{24}\,f'''(x)\right) + B\left(- \frac{h^2}{6}\,f'''(x)\right) &= 0 \\ \\ \left(-\frac{A}{24} - \frac{B}{6}\right)h^2\,f'''(x) &= 0 \\ \\ \frac{A}{24} + \frac{B}{6} &= 0 \\ \\ A = 4 &\qquad B = -1 \end{align*} $$
Combining results

$$ \begin{align*} 4\phi\left(\frac{h}{2}\right) - \phi(h) &= 4f'(x) - f'(x) = 3f'(x) \\ \\ &= 4\,\dfrac{f\left(x+h/2\right)-f\left(x-h/2\right)}{h} - \dfrac{f(x+h)-f(x-h)}{2h} + O\left(h^4\right) \end{align*} $$
Hence

$$ f'(x) = \frac{4}{3}\,\dfrac{f\left(x+h/2\right)-f\left(x-h/2\right)}{h} - \frac{1}{3}\dfrac{f(x+h)-f(x-h)}{2h} + O\left(h^4\right) $$
This gives us a truncation error $O\left(h^4\right)$ and balancing the discretization error with roundoff error yields

$$ \begin{align*} h^4 &\approx \frac{\epsilon}{h} \\ h^5 &\approx \epsilon \\ h &\approx \epsilon^{1/5} \end{align*} $$
Thus the accuracy of our numerical derivative using Richardson Extrapolation with a step size of $\epsilon^{1/5}$ is $\epsilon^{4/5}$.

In [4]: 
eps^(4/5)

ans =

   3.0002e-13

It is possible to perform more extrapolation steps with more terms of the Taylor approximation to get more accurate results. However we compute the value of our function at four points instead of two.

3.7.5 Structures

It is possible for one to define their own data structure in MATLAB. We encountered several data structures already

  • An array is a aggregate data structure that groups elements of the same type in contiguous memory and accesses the elements of the array by their index.

  • A built-in MATLAB data types are structures maintained by the MATLAB environment and language.

A structure is an aggregate data type that groups elements of different data types. These elements are called fields of the structure.

Elements of an array are accessed by their index or indices using round brackets 

>> A = [ 1 2 3; 4 5 6; 7 8 9 ]

A =

     1     2     3
     4     5     6
     7     8     9

>> A(2,3)

ans =

     6
>>


However fields of a structure are accessed by their field name . 

>> student.lastName = 'Rubino'

student = 

  struct with fields:

    lastName: 'Rubino'

>>


The string 'lastName' is the field name. The string 'Rubino' is the field value . 

>> student.firstName = 'Anthony'

student = 

  struct with fields:

     lastName: 'Rubino'
    firstName: 'Anthony'

>> student.age = 19

student = 

  struct with fields:

     lastName: 'Rubino'
    firstName: 'Anthony'
          age: 19

>>


The variable student is a MATLAB structure.

This MATLAB structure has three (3) fields identified by their field names

  1. lastName
  2. firstName
  3. age

Each of the fields of the structure has an associated field value

  1. student.lastName = 'Rubino'
  2. student. firstName = 'Anthony'
  3. student.age = 19

In order to access a specific field value one must use both the label of the structure student and the fieldname 'firstName' 

>> student.firstName

ans =

    'Anthony'

>>

3.7.6 Example of Structures

For example, we want to write a function derivative that computes the numerical derivative of a function at an array of points.

function approx = derivative(f, x, varargin)

  1. The first argument is a function handle that references a function to evaluate.
  2. The second argument is an array of double precision floating point numbers in the domain of function f . The derivative of the function f is to be approximated at these points.
  3. The third argument varargin is a $1\times n$ cell array of optional arguments.


    The output is an array of double precision floating point number with the same size as input array x that should contain the values of the numerical approximations of the derivative of f at each value in x .

We have four methods above for computing a numerical derivatives and they are very similar. We require that a default method be used to compute the derivative and suppose that our derivative function utilizes other methods we write to perform the actual computations.

We also need to determine an appropriate step size h . It is cumbersome however to pass so many arguments to all of our functions. We need create a structure to contain all of the parameters for computing a numerical derivative and pass a single structure variable instead.

We can write local functions 'forward', 'backward', 'center', and 'richardson' to compute numerical derivatives. Each of these functions has a single input argument data that contains all of the above information. Each of these functions returns the array of numerical approximations to the derivative(s). 

function approx = forward(data)
    ...
end
In [8]: 
x = pi/4;
h = 10.^(-(1:16));

exact  = cos(pi/4);
approx = derivative(@sin, pi/4, 'StepSize', h)

loglog(h,abs(exact - approx))

approx =

  Columns 1 through 10

       0.6706      0.70356      0.70675      0.70707       0.7071      0.70711      0.70711      0.70711      0.70711      0.70711

  Columns 11 through 16

      0.70711       0.7071      0.70832      0.71054      0.77716       1.1102

3.7.7 Exercises


1. Download the partially completed function derivate.m .


This partially completed function will already compute numerical approximations to the derivatives at input x using forward difference and backward difference. Finish this function by completing the local functions 

function approx = center(data)


and 

function approx = richardson(data)


Use vectorization. Notice that your code should be able to handle arrays of inputs x and arrays of inputs h .

In [9]: 
DifferenceTable(@sin, @cos, pi/4)

=====================================

             Method = forward difference
------------------------------------------------------------
    h            (f(x+h)-f(x))/h         Error
------------------------------------------------------------
 1.0e-01        0.6706029729039897      3.7e-02
 1.0e-02        0.7035594916892096      3.5e-03
 1.0e-03        0.7067531099743674      3.5e-04
 1.0e-04        0.7070714246693033      3.5e-05
 1.0e-05        0.7071032456451577      3.5e-06
 1.0e-06        0.7071064277441863      3.5e-07
 1.0e-07        0.7071067464892167      3.5e-08
 1.0e-08        0.7071067842367995      3.1e-09
 1.0e-09        0.7071068175434903      3.6e-08
 1.0e-10        0.7071077057219100      9.2e-07
 1.0e-11        0.7071121466140083      5.4e-06
 1.0e-12        0.7071010443837622      5.7e-06
 1.0e-13        0.7083222897108499      1.2e-03
 1.0e-14        0.7105427357601002      3.4e-03
 1.0e-15        0.7771561172376095      7.0e-02
 1.0e-16        1.1102230246251565      4.0e-01
------------------------------------------------------------

2. Download the function DifferenceTable.m .


This function will compute the values of a function and its exact derivative at a single point for several step sizes and compare the results. You should use this function to verify that optimum step size for both centered difference and Richardson Extrapolation.

  • Create at least two polynomial functions using polyval and two periodic trigonometric functions.

  • Create the output and plots for each of these four functions. This will be a total of three (3) tables and plots for each function; one each for 'forward', 'center', and 'richardson'.

  • Include your tables, graphs, and explanation of your results in a pdf file. You can create a word document and then use the "Export" function to create the pdf to submit with your completed derivative.m .

  • Did your numerical experiments validate our computations for the optimum step size for each method? Was the optimum step size the same for all four functions?

Submit your m-file and your pdf for the assignment.

Creative Commons Logo - White


Your use of this self-initiated mediated course material is subject to our Creative Commons License .


Creative Commons Attribution-NonCommercial-ShareAlike 4.0

Creative Commons Logo - Black
Attribution
You must give appropriate credit, provide a link to the license, and indicate if changes were made. You may do so in any reasonable manner, but not in any way that suggests the licensor endorses you or your use.

Creative Commons Logo - Black
Noncommercial
You may not use the material for commercial purposes.

Creative Commons Logo - Black
Share Alike
You are free to share, copy and redistribute the material in any medium or format. If you adapt, remix, transform, or build upon the material, you must distribute your contributions under the same license as the original.