Lecture3conv.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction to Programming in Python, Mathematics, Statistics, and Modeling\n",
    "## Lecture 3:\n",
    "This lecture contains the contents: \n",
    "* Recap: Linear Algebra, Taylor Series, self-study\n",
    "* Introduction to numerical apprixomations\n",
    "    * Numerical differentiation\n",
    "    * Numerical integration\n",
    "    * Numerical errors (Truncation error & Rounding error)\n",
    "* Simple Ordinary Differential Equation as example\n",
    "* Simple Partial Differential Equation as example\n",
    "\n",
    "This lecture will cover a re-cap of previous lecture on linear algebra and part of the self-study material on Taylor Series. Then numerical differentiation is introduced. This is applied to a simple example and it is shown how different numerical differentiation techniques affect the resulting derivative. This is followed by numerical integration, the rules are explained and an example is done. Then it is illustrated how numerical differentiation and function approximation (i.e. Taylor Series) are linked and how they apply to Ordinary Differential Equations (ODE) and to Partial Differential Equations (PDE's). Then we do an example of an ODE for which we know the anlytical solution and we compare the performance of two of our numerical schemes (here we learn the difference between an implicit and explicit scheme). Finally, we will be solving a simple PDE example (Laplace equation)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img style=\"float: left;\" src=\"lecture_03/Slide3.PNG\" width=\"100%\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img style=\"float: left;\" src=\"lecture_03/Slide4.PNG\" width=\"100%\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img style=\"float: left;\" src=\"least_squares/Slide3.PNG\" width=\"100%\">"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Short exercise to check Curve fitting exercise and reading files:\n",
    "# data is stored in file: \"data_record_ex.txt\"\n",
    "# location at which the data is recorded is stored in file: \"data_location_ex.txt\"\n",
    "\n",
    "# Load the data into the Jupyter Notebook:\n",
    "data = ???\n",
    "x_coord = ???\n",
    "\n",
    "# Using Least Squares Approximation, fit a curve to the data, particularly a constant, cosine, and sine function:\n",
    "# i.e. for: f(x) = A + B*cos(x) + C*sin(x) \n",
    "\n",
    "# Construct linear system:\n",
    "num_coef = ???\n",
    "num_points = ???\n",
    "\n",
    "A = np.ones((num_points, num_coef))\n",
    "\n",
    "for ii in range(num_points):\n",
    "    A[ii, 0] = ???\n",
    "    A[ii, 1] = ???\n",
    "    A[ii, 2] = ??? \n",
    "\n",
    "rhs = np.zeros((num_points,))\n",
    "rhs[:] = ???\n",
    "\n",
    "# Solve linear system:\n",
    "solution = np.linalg.???(???)\n",
    "\n",
    "# Create polynomial fit:\n",
    "cosines_fit = np.zeros((num_points,))\n",
    "for ii in range(num_points):\n",
    "    cosines_fit[ii] = ???\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot commands:\n",
    "plt.figure(num=None, figsize=(15, 5), dpi=80, facecolor='w', edgecolor='k')\n",
    "\n",
    "font = {'family': 'serif',\n",
    "        'color':  'darkred',\n",
    "        'weight': 'normal',\n",
    "        'size': 16,\n",
    "        }\n",
    "\n",
    "plt.plot(x_coord, data, 'or', markersize=3)\n",
    "plt.plot(x_coord, cosines_fit, 'b', markersize=3)\n",
    "plt.title('Curve fitting', fontdict=font)\n",
    "plt.xlabel('x', fontdict=font)\n",
    "plt.ylabel('data', fontdict=font)\n",
    "plt.gca().legend(('data', 'curve fit'))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Recap: Taylor Series "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Any function $ f(x) $ with $ n + 1 $ derivatives can always be expressed into a Taylor series around point x:\n",
    "\n",
    "$ f(x + \\Delta x) = f(x) + \\Delta x f'(x) + \\frac{\\Delta x^2}{2}f''(x) + \\frac{\\Delta x^3}{6}f'''(x) + \\ldots + \\frac{\\Delta x^n}{n!}f^{(n)}(x) + \\mathcal{O}(\\Delta x^{n+1})$ \n",
    "\n",
    "Note that in the above example the Taylor expansion is of \\$ f(x + \\Delta x) \\$ around \\$ x \\$, the general formula for the Taylor series is:\n",
    "\n",
    "$ f(x) = \\sum_{n=0}^{\\inf} \\frac{d^{n}f(c)}{dx^n} \\frac{(x - c)^n}{n!} $\n",
    "\n",
    "Below is an example of the Taylor series expansion for the function $ f(x) = sin(x) $, for different orders of approximation (1, 3, 5, 7, 9, 11, and 13), around $ 0 $, i.e. $ c = 0 $ (Maclaurin Series):"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example: Isothermal (fluid) compressibility"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img style=\"float: left;\" src=\"lecture_03/Slide6.PNG\" width=\"100%\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img style=\"float: left;\" src=\"lecture_03/Slide7.PNG\" width=\"100%\">"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from math import sin, cos, pi, e\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "\n",
    "# Simple example (play with compressibility, the smaller the better the approximation):\n",
    "num_p = 101\n",
    "p0 = 0\n",
    "pmax = 100\n",
    "p_vec = np.linspace(p0, pmax, num_p)\n",
    "V_0 = 1\n",
    "beta_t = 1e-2\n",
    "\n",
    "comp_act = np.zeros((num_p,))\n",
    "comp_tay = np.zeros((num_p,))\n",
    "\n",
    "### write your code here ###\n",
    "\n",
    "# Plot result and differences:\n",
    "plt.figure(num=None, figsize=(15, 5), dpi=80, facecolor='w', edgecolor='k')\n",
    "\n",
    "font = {'family': 'serif',\n",
    "        'color':  'darkred',\n",
    "        'weight': 'normal',\n",
    "        'size': 16,\n",
    "        }\n",
    "\n",
    "plt.plot(p_vec, comp_act, 'red', linewidth=2)\n",
    "plt.plot(p_vec, comp_tay, 'blue', linewidth=2)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numerical approximation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img style=\"float: left;\" src=\"lecture_04/Slide2.PNG\" width=\"100%\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## One way of looking at a numerical differentiation: Newton (or Leibniz) definition of df/dx"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img style=\"float: left;\" src=\"lecture_04/Slide3.PNG\" width=\"100%\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Another way of looking at numerical differentiation: Application of Taylor Series\n",
    "(note: $ h $ is stepsize here and is the same as $ \\Delta x $)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img style=\"float: left;\" src=\"lecture_05/Slide29.PNG\" width=\"100%\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Usually we determine the error as the **normed** difference between the **numerical approximation** and a **reference solution** (preferably the analytical solution). \n",
    "\n",
    "When the analytical function or derivatives are unknown, we still can atleast find the **order of approximation**. As can be seen in the above example, the central difference scheme is clearly a better approximation than the forward and backward schemes. When decreasing the \\$ \\Delta x \\$, the error goes to zero at a higher rate for the central scheme than for the other two. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example: Take the derivative of sin(x) & fractional flow curve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Short example of taking numerical deriviative:\n",
    "import numpy as np\n",
    "from math import sin, cos, pi\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "\n",
    "# Define which function to use:\n",
    "fraction_flow = False\n",
    "\n",
    "def my_function(x):\n",
    "    return sin(x)\n",
    "\n",
    "def my_func_der(x):\n",
    "### write your code here ###\n",
    "\n",
    "# Some paramters:\n",
    "x_left = 0\n",
    "x_right = 2 * pi\n",
    "dx = pi / 5.0\n",
    "\n",
    "num_points_num = int(np.ceil( (x_right - x_left) / dx )) + 1\n",
    "x_numeric = np.linspace(x_left, x_right, num_points_num)\n",
    "\n",
    "# Allocate memory for derivatives:\n",
    "function_ana = np.zeros((num_points_num,))\n",
    "derivative_num_forw = np.zeros((num_points_num,))\n",
    "derivative_num_centr = np.zeros((num_points_num,))\n",
    "derivative_ana_for_error = np.zeros((num_points_num,)) \n",
    "\n",
    "# Please note that the range is from 1 to (N - 2), this is to simplify indexing etc.,\n",
    "# to be precise you would need to run from 0 to (N - 1) and set rules for when ii == 0 and ii == (N - 1)\n",
    "for ii in range(1, num_points_num - 1):\n",
    "    # Calculate numerical function:\n",
    "    function_ana[ii] = my_function(x_numeric[ii])\n",
    "\n",
    "    # Calculate analytical derivative at location of numerical derivative for error calculation:\n",
    "    derivative_ana_for_error[ii] = my_func_der(x_numeric[ii])\n",
    "\n",
    "    # Calculate forward derivative:\n",
    "    derivative_num_forw[ii] = ???\n",
    "                \n",
    "    # Calculate central derivative:\n",
    "    derivative_num_centr[ii] = ???\n",
    "    \n",
    "\n",
    "# Compute error for three derivatives (NOTE: only include inside interval, excluding boundaries):\n",
    "error_forw = np.linalg.norm(derivative_num_forw[1:-1] - derivative_ana_for_error[1:-1])\n",
    "error_centr = np.linalg.norm(derivative_num_centr[1:-1] - derivative_ana_for_error[1:-1])\n",
    "\n",
    "# Create analytical solution and derivatives:\n",
    "num_point_ana = 1001\n",
    "x_analytic = np.linspace(x_left, x_right, num_point_ana)\n",
    "funcion_ana = np.zeros((num_point_ana,))\n",
    "derivative_ana = np.zeros((num_point_ana,))\n",
    "\n",
    "for ii in range(num_point_ana): \n",
    "    funcion_ana[ii] = my_function(x_analytic[ii])\n",
    "    derivative_ana[ii] = my_func_der(x_analytic[ii])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot data:\n",
    "plt.figure(num=None, figsize=(15, 5), dpi=80, facecolor='w', edgecolor='k')\n",
    "\n",
    "font = {'family': 'serif',\n",
    "        'color':  'darkred',\n",
    "        'weight': 'normal',\n",
    "        'size': 16,\n",
    "        }\n",
    "\n",
    "plt.subplot(121)\n",
    "plt.plot(x_analytic, funcion_ana, 'r', linewidth=2)\n",
    "plt.plot(x_numeric, function_ana, 'ob', markersize=3)\n",
    "plt.title('Function', fontdict=font)\n",
    "plt.xlabel('x', fontdict=font)\n",
    "plt.ylabel('function', fontdict=font)\n",
    "plt.gca().legend(('analytical','numerical')) \n",
    "\n",
    "plt.subplot(122)\n",
    "plt.plot(x_analytic, derivative_ana, 'r', markersize=3)\n",
    "plt.plot(x_numeric[1:-1], derivative_num_forw[1:-1], 'b', linewidth=2)\n",
    "plt.plot(x_numeric[1:-1], derivative_num_back[1:-1], 'g', linewidth=2)\n",
    "plt.plot(x_numeric[1:-1], derivative_num_centr[1:-1], 'black', linewidth=2)\n",
    "plt.title('Derivatives', fontdict=font)\n",
    "plt.xlabel('x', fontdict=font)\n",
    "plt.ylabel('derivative', fontdict=font)\n",
    "plt.gca().legend(('analytical',\n",
    "                  'forward, $\\epsilon$ = {:3.2f}'.format(error_forw), \n",
    "                  'backward, $\\epsilon$ = {:3.2f}'.format(error_back), \n",
    "                  'central, $\\epsilon$ = {:3.2f}'.format(error_centr)))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img style=\"float: left;\" src=\"lecture_04/Slide21.PNG\" width=\"100%\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img style=\"float: left;\" src=\"lecture_04/Slide22.PNG\" width=\"100%\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A simple example in Python:\n",
    "**Use n = 5, 10, and 100 points to calculate the integral and compare the result.**\n",
    "\n",
    "$ \\int_{-1}^{1}{\\sqrt{1 + x^3}dx} \\approx 1.9527572 $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Note this is a Python code block. Write the code for the numerical evaluation of the integral here\n",
    "import numpy as np\n",
    "\n",
    "\n",
    "def my_func(x):\n",
    "    return np.sqrt(1 + x**3)\n",
    "\n",
    "num_p = 5\n",
    "int_low = -1\n",
    "int_high = 1\n",
    "int_len = int_high - int_low\n",
    "dx = int_len / num_p\n",
    "\n",
    "int_left = 0\n",
    "int_right = 0\n",
    "int_mid = 0\n",
    "int_trap = 0\n",
    "\n",
    "a = int_low\n",
    "\n",
    "for ii in range(num_p):\n",
    "    \n",
    "### write your code here ###\n",
    "    \n",
    "print(\"Integral with Left Riemann sum:\\t\\t\", int_left)\n",
    "print(\"Integral with Right Riemann sum:\\t\", int_right)\n",
    "print(\"Integral with Midpoint rule:\\t\\t\", int_mid)\n",
    "print(\"Integral with Trapezoidal:\\t\\t\", int_trap)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img style=\"float: left;\" src=\"lecture_05/Slide30.PNG\" width=\"100%\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Note:** In the above example the approximation error only conists of the truncation error, in reality there is also a round-off error due to numerical rounding off. The round-off error typically increases when decreasing the stepsize. For real-world Geo-Engineering problems the round-off error is several order of mangitude smaller than the truncation error though! "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Ordinary Differential Equations\n",
    "\n",
    "In this section we will go through an example of a ODE. An ODE is a special case of the PDE, particularly it's a Differential Equation with only one indepdent variable. In constrast, a PDE is a Differential Equation with more than one indepdent variables. \n",
    "\n",
    "The example we will go through is the following:\n",
    "* $ y' = -100(y - cos(t)) - sin(t) = f(t, y(t)) $, with initial conditions: $ y(0) = 0 $.\n",
    "\n",
    "This equation has the analytical solution: \n",
    "* $ y(t) = cos(t) -e^{-100t} $\n",
    "\n",
    "**!!!** Convince yourself that this is true by filling the differential equation in (take the derivative and substitute back in the equation, also check for the initial condition) **!!!**\n",
    "\n",
    "**In some cases the analytical solution is impossible to derrive and we have to use numerical methods to find the solution of our function over time** (or e.g. space-time in case of some PDEs).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Explicit and Implicit strategies for solving the ODE in time:\n",
    "First, let us realize that the analytical solution at an point in time, $ t $, can be found by integrating the above equation in time, which results in the following integral equation: \n",
    "\n",
    "* $ y(t) = y(t_0) + \\int_{t_0}^{t}{f(\\tau, y(\\tau))d\\tau} $\n",
    "\n",
    "which in return can be written as:\n",
    "\n",
    "* $ y(t_{n+1}) = y(t_n) + \\int_{t_n}^{t_{n+1}}{f(\\tau, y(\\tau))d\\tau} $\n",
    "\n",
    "Our main objective is to approximate the integral. We will explain later in the lecture how to do numerical integration, but let's say that we use the **left-hand rectangle** rule which will give us the following approximation:\n",
    "\n",
    "* $ y(t_{n+1}) = y(t_n) + \\Delta t f(t_{n}, y(t_{n})) $\n",
    "\n",
    "This approximation is also known as the Forward Euler (or Explicit Euler) method.\n",
    "\n",
    "We could also use the **right-hand rectangle** rule which would have given us the following approximation:\n",
    "\n",
    "* $ y(t_{n+1}) = y(t_n) + \\Delta t f(t_{n+1}, y(t_{n+1})) $\n",
    "\n",
    "This approximation is also known as the Backward Euler (or Implicit Euler) method.\n",
    "\n",
    "**NOTE:** The main difference between the two is that in the Explicit method, the right-hand side (rhs) of the equation, $ f(t,y(t)) $ , is evaluated at $ t_n $ . While in the Implicit method, the rhs of the equation is evaluated at $ t_n+t $ .\n",
    "\n",
    "Please also note that we could have also derrived the Forward Euler method by using the Taylor series expansion of the function $ y(t + \\Delta t) $ in the neighbourhood of $ t_n $:\n",
    "\n",
    "$ y(t_n + \\Delta t) = y(t_n) + \\Delta t \\frac{dy(t_n)}{dt_n} + \\mathcal{O}(\\Delta t^2) = y(t_n) + \\Delta t f(t_{n}, y(t_{n})) + \\mathcal{O}(\\Delta t^2) $\n",
    "\n",
    "Since we know from the original equation that $ \\frac{dy}{dt} = f(t, y(t)) $. \n",
    "\n",
    "Another way is to approximate the derivative and evaluating the rhs at either $ t_n \\$ or \\$ t_{n+1} $: \n",
    "\n",
    "$ \\frac{dy}{dt} \\approx \\frac{y(t_{n+1}) - y(t_{n})}{\\Delta t} = f(t_{n}, y(t_{n})) $ (Forward Euler)\n",
    "\n",
    "$ \\frac{dy}{dt} \\approx \\frac{y(t_{n+1}) - y(t_{n})}{\\Delta t} = f(t_{n+1}, y(t_{n+1})) $ (Backward Euler)\n",
    "\n",
    "### Solving the actual problem:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Problem statement:\n",
    "y' = -100(y - cos(t)) - sin(t)\n",
    "y(0) = 0\n",
    "\n",
    "Solve with Forward and Backward Euler schemes, stepsize of 0.2 for a total time of 4. Plot both solutions on the\n",
    "interval [0, 4], on top of the exact solution which is given by:\n",
    "y(t) = cos(t) - e^(-100t)\n",
    "\n",
    "Forward Euler: y_n+1 = y_n + dt*f(t_n, y_n)\n",
    "Backward Euler: y_n+1 = y_n + dt*f(t_n+1, y_n+1)\n",
    "\"\"\"\n",
    "import numpy as np\n",
    "from math import e, cos, sin, ceil\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "\n",
    "# Define parameters:\n",
    "dt = 0.02\n",
    "tot_time = 4  \n",
    "num_steps = ceil(tot_time / dt)\n",
    "exact_steps = min(50000, 100 * num_steps)\n",
    "time_vec = np.linspace(0, tot_time, num_steps+1)\n",
    "time_vec_exact = np.linspace(0, tot_time, exact_steps+1)\n",
    "\n",
    "# Allocate mememory:\n",
    "sol_forward = np.zeros((num_steps+1,))\n",
    "sol_backard = np.zeros((num_steps+1,))\n",
    "\n",
    "# Perform time-integration\n",
    "for ii in range(1, num_steps+1):\n",
    "    # Forward Euler:\n",
    "    sol_forward[ii] = ???\n",
    "\n",
    "    # Backward Euler:\n",
    "    sol_backard[ii] = ???\n",
    "\n",
    "# Calculate exact solution:\n",
    "sol_exact = np.zeros((exact_steps+1,))\n",
    "for ii in range(exact_steps+1):\n",
    "    sol_exact[ii] = ???\n",
    "\n",
    "# plot with various axes scales\n",
    "plt.figure(num=None, figsize=(15, 5), dpi=80, facecolor='w', edgecolor='k')\n",
    "\n",
    "font = {'family': 'serif',\n",
    "        'color':  'darkred',\n",
    "        'weight': 'normal',\n",
    "        'size': 16,\n",
    "        }\n",
    "\n",
    "# Forward Euler\n",
    "plt.subplot(121)\n",
    "plt.plot(time_vec, sol_forward, 'r', time_vec_exact, sol_exact, 'g--')\n",
    "plt.title('Forward Euler dt = ' + str(dt), fontdict=font)\n",
    "plt.xlabel('t', fontdict=font)\n",
    "plt.ylabel('y', fontdict=font)\n",
    "plt.gca().legend(('forward','analytic'))\n",
    "\n",
    "plt.subplot(122)\n",
    "plt.plot(time_vec, sol_backard, 'r', time_vec_exact, sol_exact, 'g--')\n",
    "plt.title('Backward Euler dt = ' + str(dt), fontdict=font)\n",
    "plt.xlabel('t', fontdict=font)\n",
    "plt.ylabel('y', fontdict=font)\n",
    "plt.gca().legend(('backward','analytic'))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}