topic_06_numerical_differentiation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Numerical differentiation\n",
    "\n",
    "In this lesson, we will estimate the first and second derivatives of functions using different methods. By the end of the lesson, you will be able to:\n",
    "\n",
    "1.\tfind first derivatives using backward, forward, and central difference methods\n",
    "2.\tfind second derivatives using central difference method.\n",
    "\n",
    "In engineering, we frequently want to know the rate of change of a quantity with respect to time, or distance, or other quantity. For example, velocity is the rate of change of displacement with respect to time. Differentiation is determining the rate of change of a quantity. \n",
    "For relatively straightforward equations, we can determine the exact derivative of the quantity - you did this when doing the Newton's method problems. It is common that the exact derivative is not possible to find. Examples would include when there is no equation describing the behaviour of the quantity in question - it may be in the form of an experimental data set. In those cases, numerical differentiation is used to find an approximation to the derivative. We will examine several methods for estimating the first derivative and one method for estimating the second derivative.\n",
    "\n",
    "## First derivative methods\n",
    "\n",
    "### Forward difference\n",
    "Consider a function $f(x)$. By expanding the function in a Taylor's series, it can be shown that the first derivative of the function at $x=a$, $f'(a)$ is approximated by\n",
    "\n",
    "$$f'(a)=\\frac{f(a+h)-f(a)}{h}$$\n",
    " \n",
    "where $h$ is a small increment.\n",
    "\n",
    "The smaller $h$ is the more accurate the approximation $f'(a)$ is. The error of the approximation is on the order of $h$.\n",
    "\n",
    "Consider a function $f(x)=x^2+3x+4$ . Let's approximate the derivative of this function at $x=3$  using the forward difference method.\n",
    "\n",
    "Copy the script and two functions below into a new Python script file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x\tForward\n",
      "3.0\t9.0010\n"
     ]
    }
   ],
   "source": [
    "def f(x):\n",
    "    y = x**2+3*x+4\n",
    "    return y\n",
    "\n",
    "def forward_diff(x, h):\n",
    "    y = (f(x+h) - f(x)) / h\n",
    "    return y\n",
    "\n",
    "# Script to numerically differentiate a function using forward differencing\n",
    "x = 3\n",
    "h = 0.001\n",
    "\n",
    "forward = forward_diff(x,h)\n",
    "print('x\\tForward')\n",
    "print(f'{x:2.1f}\\t{forward:4.4f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise\n",
    "Modify the script above to include the exact derivative of the function at $x=3$. Include the exact derivative in the table.\n",
    "We will now look at two other methods of estimating the first derivative of a function.\n",
    "\n",
    "### Backward difference\n",
    "Consider a function $f(x)$. By expanding the function in a Taylor's series, it can be shown that the first derivative of the function at $x=a$, $f'(a)$ is approximated by\n",
    "\n",
    "$$f'(a)=\\frac{f(a)-f(a-h)}{h}$$\n",
    " \n",
    "where $h$ is a small increment.\n",
    "\n",
    "The smaller $h$ is the more accurate the approximation $f'(a)$ is. The error of the approximation is on the order of $h$ like the forward difference method.\n",
    "\n",
    "## Exercise\n",
    "Modify your script above to approximate the derivative using the backward difference method. Add the back difference result to the table output to the shell.\n",
    "\n",
    "### Central difference\n",
    "Consider a function $f(x)$. By expanding the function in a Taylor's series, it can be shown that the first derivative of the function at $x=a$, $f'(a)$ is approximated by\n",
    "\n",
    "$$f'(a)=\\frac{f(a+h)-f(a-h)}{2h}$$\n",
    " \n",
    "where $h$ is a small increment.\n",
    "\n",
    "The smaller $h$ is the more accurate the approximation $f'(a)$ is. The error of the approximation for this method is on the order of $h^2$. So it is more accurate than both the forward and backward methods. This is indicated in the plot below, which shows the numerical differentiation of the function at $x$ using each method.\n",
    "<img src=\"./img/diff_plot.png\" alt=\"Comparing numerical differentiation methods\" align='center' style=\"width: 600px;\"/>\n",
    "## Exercise\n",
    "Modify your script above to approximate the derivative using the central difference method. Add the central difference result to the table output to the shell.\n",
    "## Problem 1\n",
    "Consider the function\n",
    "$$f(x)=\\sqrt{1+x^2}e^{\\sqrt{1+x^2}}$$\n",
    " \n",
    "The exact derivative is given by\n",
    "$$f'(x)=xe^{\\sqrt{1+x^2}}\\left(1+\\frac{1}{\\sqrt{1+x^2}}\\right)$$\n",
    "Find $f'(1.3)$ using the backward, forward, and central difference methods. Find the % error for each method. Output your results in the form of a table. Take an interval size $h=0.1$.\n",
    "\n",
    "Examine and discuss what happens when the interval size is increased and decreased.\n",
    "\n",
    "## Problem 2\n",
    "The distribution of the velocity $u$ across the cross-section of a straight pipe for steady flow is given by Poiseuille's equation:\n",
    "\n",
    "$$u(r)=\\frac{k_s}{4\\mu}\\left(r^2-a^2\\right)$$\n",
    " \n",
    "where $k_s=-200$ Pa/m is the pressure gradient; $\\mu=0.0035$ Pas is the viscosity; $a=15$ mm is the radius.\n",
    "The shear stress in the fluid is given by\n",
    "$$\\tau=\\mu\\frac{du}{dr}$$\n",
    "Find the shear stress on the inner wall by the forward, backward, and central divided difference methods. Compare those estimates with the exact solution. Let $h=0.001$.\n",
    "\n",
    "## Problem 3\n",
    "A radar station is tracking the motion of an aircraft. The time $t$, recorded distance to the aircraft $r$, and the angle $\\theta$ during a period of 24 s is given below:\n",
    "\n",
    "Time (s):   0       4       8       12       16      20       24\n",
    "\n",
    "$r$  (km):     18.803  18.861  18.946  19.042   19.148  19.260   19.376\n",
    "  \n",
    "$\\theta$  (rad):    0.7854  0.7792  0.7701  0.7594   0.7477  0.7350   0.7215\n",
    "\n",
    "The angle of instantaneous velocity can be calculated by:\n",
    "$$v=\\sqrt{\\left(\\frac{dr}{dt}\\right)^2+\\left(r\\frac{d\\theta}{dt}\\right)^2}$$\n",
    "Find the velocity at 0 using forward difference method, 12 s using central difference, and 24 s using the backward difference method. Take $h=0.01$. Note: you will need to use the interplation function in python to find $r$ and $\\theta$ at $t-h$ and $t+h$.\n",
    "\n",
    "## Second derivative methods\n",
    "Consider a function $f(x)$. By expanding the function in a Taylor's series, it can be shown that the first derivative of the function at $x=a$, $f''(a)$ is approximated by\n",
    "\n",
    "$$f''(a)=\\frac{f(a+h)-2f(a)+f(a-h)}{h^2}$$\n",
    " \n",
    "where $h$ is a small increment.\n",
    "\n",
    "## Exercise\n",
    "Approximate the second derivative of the following function at $x=1.5$.\n",
    "\n",
    "$$f(x)=x^4+3x^3+5x^2-4x+10$$\n",
    " \n",
    "Compare the approximation to the exact derivative at $x=1.5$. Let $h=0.001$.\n",
    "\n",
    "## Problem 4\n",
    "A 1 m long uniform beam is simply supported at both ends and is subjected to a load. The deflection of the beam is given by the following differential equation:\n",
    "\n",
    "$$\\frac{d^2y}{dx^2}=-\\frac{M}{EI}$$\n",
    "\n",
    "where $y$ is the deflection, $x$ is the distance along the beam, $M$ is the bending moment, and $EI$ is the flexural rigidity of the beam.\n",
    "\n",
    "The following experimental data is obtained from measuring the deflection of the beam along its length.\n",
    "\n",
    "$x$  (m)   0   0.2     0.4     0.6       0.8     1.0\n",
    "  \n",
    "$y$  (m)   0   0.0778  0.1068  0.0838    0.0397  0\n",
    "  \n",
    "Find the bending moment at 0.4 m and 0.6 m along the beam.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}