{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "# App. A: Sequences and difference equations\n", "\n", " **Hans Petter Langtangen**, Simula Research Laboratory and University of Oslo, Dept. of Informatics\n", "\n", "Date: **Aug 21, 2016**\n", "\n", "## Sequences\n", "\n", "*Sequences* is a central topic in mathematics:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "x_0,\\ x_1,\\ x_2,\\ \\ldots,\\ x_n,\\ldots,\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Example: all odd numbers" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "1, 3, 5, 7, \\ldots, 2n+1,\\ldots\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For this sequence we have a formula for the $n$-th term:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "x_n = 2n+1\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "and we can write the sequence more compactly as" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "(x_n)_{n=0}^\\infty,\\quad x_n = 2n+1\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Other examples of sequences" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "1,\\ 4,\\ 9,\\ 16,\\ 25,\\ \\ldots\\quad (x_n)_{n=0}^\\infty,\\ x_n=n^2\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "1,\\ {1\\over 2},\\ {1\\over3},\\ {1\\over4},\\ \\ldots\\quad (x_n)_{n=0}^\\infty,\\ x_n={1\\over {n+1}}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "1,\\ 1,\\ 2,\\ 6,\\ 24,\\ \\ldots\\quad (x_n)_{n=0}^\\infty,\\ x_n=n!\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "1,\\ 1+x,\\ 1+x+{1\\over2}x^2,\\ 1+x+{1\\over2}x^2+{1\\over6}x^3,\\ \\ldots\\quad (x_n)_{n=0}^\\infty,\\ x_n=\\sum_{j=0}^n {x^j\\over j!}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Finite and infinite sequences\n", "\n", " * Infinite sequences have an infinite number of terms ($n\\rightarrow\\infty$)\n", "\n", " * In mathematics, infinite sequences are widely used\n", "\n", " * In real-life applications, sequences are usually finite: $(x_n)_{n=0}^N$\n", "\n", " * Example: number of approved exercises every week in INF1100\n", " $x_0, x_1, x_2, \\ldots, x_{15}$\n", "\n", " * Example: the annual value of a loan\n", " $x_0, x_1, \\ldots, x_{20}$\n", "\n", "\n", "\n", "## Difference equations\n", "\n", " * For sequences occuring in modeling of real-world phenomena, there is seldom a formula for the $n$-th term\n", "\n", " * However, we can often set up one or more equations governing the sequence\n", "\n", " * Such equations are called difference equations\n", "\n", " * With a computer it is then very easy to generate the sequence by solving the difference equations\n", "\n", " * Difference equations have lots of applications and are very easy to solve on a computer, but often complicated or impossible to solve for $x_n$ (as a formula) by pen and paper!\n", "\n", " * The programs require only loops and arrays\n", "\n", "\n", "\n", "## Modeling interest rates\n", "\n", "**Problem:**\n", "\n", "Put $x_0$ money in a bank at year 0. What is the value after $N$ years if the interest rate is $p$ percent per year?\n", "\n", "\n", "\n", "**Solution:**\n", "\n", "The fundamental information relates the value at year $n$, $x_n$, to the value of the previous year, $x_{n-1}$:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "x_{n} = x_{n-1} + {p\\over 100}x_{n-1}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "How to solve for $x_n$? Start with $x_0$, compute $x_1$, $x_2$, ...\n", "\n", "\n", "\n", "## Simulating the difference equation for interest rates\n", "\n", "**What does it mean to simulate?**\n", "\n", "Solve math equations by repeating a simple procedure (relation) many times (boring, but well suited for a computer!)\n", "\n", "\n", "\n", "**Program for $x_{n} = x_{n-1} + (p/100)x_{n-1}$:**" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline\n", "\n", "from scitools.std import *\n", "x0 = 100 # initial amount\n", "p = 5 # interest rate\n", "N = 4 # number of years\n", "index_set = range(N+1)\n", "x = zeros(len(index_set))\n", "\n", "# Solution:\n", "x[0] = x0\n", "for n in index_set[1:]:\n", " x[n] = x[n-1] + (p/100.0)*x[n-1]\n", "print x\n", "plot(index_set, x, 'ro', xlabel='years', ylabel='amount')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## We do not need to store the entire sequence, but it is convenient for programming and later plotting\n", "\n", " * Previous program stores all the $x_n$ values in a NumPy array\n", "\n", " * To compute $x_n$, we only need one previous value, $x_{n-1}$\n", "\n", "Thus, we could only store the two last values in memory:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [], "source": [ "x_old = x0\n", "for n in index_set[1:]:\n", " x_new = x_old + (p/100.)*x_old\n", " x_old = x_new # x_new becomes x_old at next step" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "However, programming with an array `x[n]` is simpler, safer, and enables plotting the sequence, so we will continue to use arrays in the examples\n", "\n", "\n", "\n", "## Daily interest rate\n", "\n", " * A more relevant model is to add the interest every day\n", "\n", " * The interest rate per day is $r=p/D$ if $p$ is the annual interest rate and $D$ is the number of days in a year\n", "\n", " * A common model in business applies $D=360$, but $n$ counts exact (all) days\n", "\n", "Just a minor change in the model:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "x_{n} = x_{n-1} + {r\\over 100}x_{n-1}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "How can we find the number of days between two dates?" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import datetime\n", "date1 = datetime.date(2007, 8, 3) # Aug 3, 2007\n", "date2 = datetime.date(2008, 8, 4) # Aug 4, 2008\n", "diff = date2 - date1\n", "print diff.days" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Program for daily interest rate" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from scitools.std import *\n", "x0 = 100 # initial amount\n", "p = 5 # annual interest rate\n", "r = p/360.0 # daily interest rate\n", "import datetime\n", "date1 = datetime.date(2007, 8, 3)\n", "date2 = datetime.date(2011, 8, 3)\n", "diff = date2 - date1\n", "N = diff.days\n", "index_set = range(N+1)\n", "x = zeros(len(index_set))\n", "\n", "# Solution:\n", "x[0] = x0\n", "for n in index_set[1:]:\n", " x[n] = x[n-1] + (r/100.0)*x[n-1]\n", "print x\n", "plot(index_set, x, 'ro', xlabel='days', ylabel='amount')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## But the annual interest rate may change quite often...\n", "\n", "**Varying $p$ means $p_n$:**\n", "\n", " * Could not be handled in school (cannot apply $x_n =x_0(1+\\frac{p}{100})^n$)\n", "\n", " * A varying $p$ causes no problems in the program:\n", " just fill an array `p` with correct interest rate for day `n`\n", "\n", "Modified program:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [], "source": [ "p = zeros(len(index_set))\n", "# fill p[n] for n in index_set (might be non-trivial...)\n", "\n", "r = p/360.0 # daily interest rate\n", "x = zeros(len(index_set))\n", "\n", "x[0] = x0\n", "for n in index_set[1:]:\n", " x[n] = x[n-1] + (r[n-1]/100.0)*x[n-1]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Payback of a loan\n", "\n", " * A loan $L$ is paid back with a fixed amount $L/N$ every month over $N$ months + the interest rate of the loan\n", "\n", " * $p$: annual interest rate, $p/12:$ monthly rate\n", "\n", " * Let $x_n$ be the value of the loan at the end of month $n$\n", "\n", "The fundamental relation from one month to the text:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "x_n = x_{n-1} + {p\\over 12\\cdot 100}x_{n-1} -\n", "({p\\over 12\\cdot 100}x_{n-1} + {L\\over N})\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "which simplifies to" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "x_n = x_{n-1} - {L\\over N}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "($L/N$ makes the equation *nonhomogeneous*)\n", "\n", "\n", "\n", "## How to make a living from a fortune with constant consumption\n", "\n", " * We have a fortune $F$ invested with an annual interest rate of $p$ percent\n", "\n", " * Every year we plan to consume an amount $c_n$ ($n$ counts years)\n", "\n", " * Let $x_n$ be our fortune at year $n$\n", "\n", "A fundamental relation from one year to the other is" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "x_n = x_{n-1} + {p\\over 100}x_{n-1} - c_n\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Simplest possibility: keep $c_n$ constant, but inflation\n", "demands $c_n$ to increase...\n", "\n", "\n", "\n", "## How to make a living from a fortune with inflation-adjusted consumption\n", "\n", " * Assume $I$ percent inflation per year\n", "\n", " * Start with $c_0$ as $q$ percent of the interest the first year\n", "\n", " * $c_n$ then develops as money with interest rate $I$\n", "\n", "$x_n$ develops with rate $p$ but with a loss $c_n$ every year:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\begin{align*}\n", "x_n &= x_{n-1} + {p\\over 100}x_{n-1} - c_{n-1}, \\quad x_0=F,\\ c_0 = {pq\\over 10^4}F\\\\\n", "c_n &= c_{n-1} + {I\\over100}c_{n-1}\n", "\\end{align*}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This is a coupled system of *two* difference equations, but the programming\n", "is still simple: we update two arrays, first `x[n]`, then `c[n]`, inside the loop\n", "(good exercise!)\n", "\n", "\n", "\n", "## The mathematics of Fibonacci numbers\n", "\n", "No programming or math course is complete without an example on Fibonacci numbers:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "x_n = x_{n-1} + x_{n-2},\\quad x_0=1,\\ x_1=1\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Mathematical classification.**\n", "\n", "This is a *homogeneous difference equation of second order* (second order means three levels: $n$, $n-1$, $n-2$). This classification is important for mathematical solution technique, but not for simulation in a program.\n", "\n", "\n", "\n", "Fibonacci derived the sequence by modeling rat populations, but the sequence of numbers has a range of peculiar mathematical properties and has therefore attracted much attention from mathematicians.\n", "\n", "\n", "\n", "## Program for generating Fibonacci numbers" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [], "source": [ "N = int(sys.argv[1])\n", "from numpy import zeros\n", "x = zeros(N+1, int)\n", "x[0] = 1\n", "x[1] = 1\n", "for n in range(2, N+1):\n", " x[n] = x[n-1] + x[n-2]\n", " print n, x[n]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Fibonacci numbers can cause overflow in NumPy arrays\n", "\n", "Run the program with $N=50$:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [], "source": [ "2 2\n", "3 3\n", "4 5\n", "5 8\n", "6 13\n", "...\n", "45 1836311903\n", "Warning: overflow encountered in long_scalars\n", "46 -1323752223" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note:\n", "\n", " * Changing `int` to `long` or `int64` for array elements allows $N\\leq 91$\n", "\n", " * Can use `float96` (though $x_n$ is integer): $N \\leq 23600$\n", "\n", "\n", "\n", "## No overflow when using Python int types\n", "\n", " * Best: use Python scalars of type `int` - these automatically changes to `long` when overflow in `int`\n", "\n", " * The `long` type in Python has arbitrarily many digits (as many as required in a computation!)\n", "\n", " * Note: `long` for arrays is 64-bit integer (`int64`), while scalar `long` in Python is an integer with as \"infinitely\" many digits\n", "\n", "\n", "\n", "## Program with Python's int type for integers\n", "\n", "The program now avoids arrays and makes use of three `int` objects (which automatically changes to `long` when needed):" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import sys\n", "N = int(sys.argv[1])\n", "xnm1 = 1 # \"x_n minus 1\"\n", "xnm2 = 1 # \"x_n minus 2\"\n", "n = 2\n", "while n <= N:\n", " xn = xnm1 + xnm2\n", " print 'x_%d = %d' % (n, xn)\n", " xnm2 = xnm1\n", " xnm1 = xn\n", " n += 1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Run with $N=200$:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": false }, "outputs": [], "source": [ "x_2 = 2\n", "x_3 = 3\n", "...\n", "x_198 = 173402521172797813159685037284371942044301\n", "x_199 = 280571172992510140037611932413038677189525\n", "x_200 = 453973694165307953197296969697410619233826" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Limition: your computer's memory\n", "\n", "\n", "\n", "## New problem setting: exponential growth with limited environmental resources\n", "\n", "The model for growth of money in a bank has a solution of the type" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "x_n = x_0C^n \\quad (= x_0e^{n\\ln C})\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note:\n", "\n", " * This is exponential growth in time ($n$)\n", "\n", " * Populations of humans, animals, and cells also exhibit the same type of growth as long as there are unlimited resources (space and food)\n", "\n", " * Most environments can only support a maximum number $M$ of individuals\n", "\n", " * How can we model this limitation?\n", "\n", "\n", "\n", "## Modeling growth in an environment with limited resources\n", "\n", "Initially, when there are enough resources, the growth is exponential:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "x_n = x_{n-1} + {r\\over 100}x_{n-1}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The growth rate $r$ must decay to zero as $x_n$ approaches $M$. The simplest\n", "variation of $r(n)$ is a linear:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "r(n) = \\varrho \\left(1 - {x_n\\over M}\\right)\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Observe: $r(n)\\approx \\varrho$ for small $n$ when $x_n\\ll M$, and\n", "$r(n) \\rightarrow 0$ as $x_n\\rightarrow M$ and $n$ is big\n", "\n", "\n", "\n", "**Logistic growth model:**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "x_n = x_{n-1} + {\\varrho\\over 100} x_{n-1}\\left(1 - {x_{n-1}\\over M}\\right)\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "(This is a *nonlinear* difference equation)\n", "\n", "\n", "\n", "\n", "## The evolution of logistic growth\n", "\n", "\n", "In a program it is easy to introduce logistic instead of exponential growth, just replace" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": false }, "outputs": [], "source": [ "x[n] = x[n-1] + p/100.0)*x[n-1]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "by" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": false }, "outputs": [], "source": [ "x[n] = x[n-1] + (rho/100.0)*x[n-1]*(1 - x[n-1]/float(M))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "\n", "
\n", "\n", "\n", "\n", "\n", "\n", "## The factorial as a difference equation\n", "\n", "The factorial $n!$ is defined as" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "n(n-1)(n-2)\\cdots 1,\\quad 0!=1\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The following difference equation has $x_n=n!$ as solution and can be used to compute the factorial:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "x_n = nx_{n-1},\\quad x_0 = 1\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Difference equations must have an initial condition\n", "\n", " * In mathematics, it is much stressed that a difference equation for\n", " $x_n$ must have an initial condition $x_0$\n", "\n", " * The initial condition is obvious when programming: otherwise we cannot\n", " start the program ($x_0$ is needed to compute $x_n$)\n", "\n", " * However: if you forget `x[0] = x0` in the program, you get $x_0=0$\n", " (because `x = zeroes(N+1)`),\n", " which (usually) gives unintended results!\n", "\n", "\n", "\n", "## Have you ever though about how $\\sin x$ is really calculated?\n", "\n", " * How can you *calculate* $\\sin x$, $\\ln x$, $e^x$ without a calculator\n", " or program?\n", "\n", " * These functions were originally defined to have some desired mathematical\n", " properties, but without an algorithm for how to evaluate function values\n", "\n", " * Idea: approximate $\\sin x$, etc. by polynomials, since they are easy\n", " to calculate (sum, multiplication), but how??\n", "\n", "\n", "\n", "## Would you expect these fantastic mathematical results?\n", "\n", "**Amazing result by Gregory, 1667:**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\sin x = \\sum_{k=0}^\\infty (-1)^k\\frac{x^{2k+1}}{(2k+1)!}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Even more amazing result by Taylor, 1715:**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "f(x) = \\sum_{k=0}^\\infty \\frac{1}{k!}(\\frac{d^k}{dx^k} f(0))x^k\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For \"any\" $f(x)$, if we can differentiate, add, and multiply $x^k$,\n", "we can evaluate $f$ at any $x$ (!!!)\n", "\n", "\n", "\n", "## Taylor polynomials\n", "\n", "**Practical applications works with a truncated sum:**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "f(x) \\approx \\sum_{k=0}^{{\\color{red} N}}\n", "\\frac{1}{k!}(\\frac{d^k}{dx^k} f(0))x^k\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$N=1$ is *very* popular and has been essential in developing\n", "physics and technology\n", "\n", "\n", "\n", "**Example:**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\begin{align*}\n", "e^x &= \\sum_{k=0}^\\infty \\frac{x^k}{k!}\\\\\n", "&\\approx 1 + x + \\frac{1}{2}x^2 + \\frac{1}{6}x^3\\\\\n", "&\\approx 1 + x\n", "\\end{align*}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Taylor polynomials around an arbitrary point\n", "\n", "The previous Taylor polynomials are most accurate around $x=0$. Can\n", "make the polynomials accurate around any point $x=a$:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "f(x) \\approx \\sum_{k=0}^N \\frac{1}{k!}(\\frac{d^k}{dx^k} f(a))(x-a)^k\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Taylor polynomial as one difference equation\n", "\n", "\n", "The Taylor series for $e^x$ around $x=0$ reads" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "e^x= \\sum_{n=0}^\\infty {x^n\\over n!}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Define" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "e_{n}=\\sum_{k=0}^{n-1} \\frac{x^k}{k!} =\n", "\\sum_{k=0}^{n-2} \\frac{x^k}{k!} + \\frac{x^{n-1}}{(n-1)!}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can formulate the sum in $e_n$ as the following difference equation:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "e_n = e_{n-1} + \\frac{x^{n-1}}{(n-1)!},\\quad e_0=0\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## More efficient computation: the Taylor polynomial as two difference equations\n", "\n", "Observe:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\frac{x^n}{n!} = \\frac{x^{n-1}}{(n-1)!}\\cdot \\frac{x}{n}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let $a_n = x^n/n!$. Then we can efficiently compute $a_n$ via" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "a_n = a_{n-1}\\frac{x}{n},\\quad a_0 = 1\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we can update each term via the $a_n$ equation and sum the terms via\n", "the $e_n$ equation:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\begin{align*}\n", "e_n &= e_{n-1} + a_{n-1},\\quad e_0 = 0,\\ a_0 = 1\\\\\n", "a_n &= \\frac{x}{n} a_{n-1}\n", "\\end{align*}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "See the book for more details\n", "\n", "\n", "\n", "## Nonlinear algebraic equations\n", "\n", "**Generic form of any (algebraic) equation in $x$:**\n", "\n", "\\[ f(x)=0 \\]\n", "\n", "\n", "\n", "**Examples that can be solved by hand:**\n", "\n", " \\[ ax + b = 0 \\]\n", " \\[ ax^2 + bx + c = 0 \\]\n", " \\[ \\sin x + \\cos x = 1 \\]\n", "\n", "\n", "\n", " * Simple numerical algorithms can solve \"any\" equation $f(x)=0$\n", "\n", " * Safest: Bisection\n", "\n", " * Fastest: Newton's method\n", "\n", " * Don't like $f'(x)$ in Newton's method? Use the Secant method\n", "\n", " * Secant and Newton are difference equations!\n", "\n", "\n", "\n", "## Newton's method for finding zeros; illustration\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "## Newton's method for finding zeros; mathematics\n", "\n", "**Newton's method.**\n", "\n", "Simpson (1740) came up with the following general\n", "method for solving $f(x)=0$ (based on ideas by Newton):" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "x_n = x_{n-1} - {f(x_{n-1})\\over f'(x_{n-1})},\\quad x_0 \\hbox{ given}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note:\n", "\n", " * This is a (nonlinear!) difference equation\n", "\n", " * As $n\\rightarrow\\infty$, we hope that $x_n\\rightarrow x_s$, where $x_s$ solves $f(x_s)=0$\n", "\n", " * How to choose $N$ when what we want is $x_N$ close to $x_s$?\n", "\n", " * Need a slightly different program: simulate until $f(x)\\leq \\epsilon$, where $\\epsilon$ is a small tolerance\n", "\n", " * Caution: Newton's method may (easily) diverge, so $f(x)\\leq\\epsilon$ may\n", " never occur!\n", "\n", "\n", "\n", "## A program for Newton's method\n", "\n", "**Quick implementation:**" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def Newton(f, x, dfdx, epsilon=1.0E-7, max_n=100):\n", " n = 0\n", " while abs(f(x)) > epsilon and n <= max_n:\n", " x = x - f(x)/dfdx(x)\n", " n += 1\n", " return x, n, f(x)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note:\n", "\n", " * $f(x)$ is evaluated twice in each pass of the loop - only one evaluation is strictly necessary (can store the value in a variable and reuse it)\n", "\n", " * `f(x)/dfdx(x)` can give integer division\n", "\n", " * It could be handy to store the `x` and `f(x)` values in each iteration (for plotting or printing a convergence table)\n", "\n", "\n", "\n", "## An improved function for Newton's method\n", "\n", "Only one $f(x)$ call in each iteration, optional storage of $(x,f(x))$ values during the iterations, and ensured float division:" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def Newton(f, x, dfdx, epsilon=1.0E-7, max_n=100,\n", " store=False):\n", " f_value = f(x)\n", " n = 0\n", " if store: info = [(x, f_value)]\n", " while abs(f_value) > epsilon and n <= max_n:\n", " x = x - float(f_value)/dfdx(x)\n", " n += 1\n", " f_value = f(x)\n", " if store: info.append((x, f_value))\n", " if store:\n", " return x, info\n", " else:\n", " return x, n, f_value" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Application of Newton's method" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "e^{-0.1x^2}\\sin ({\\pi\\over 2}x) =0\n", "Solutions: $x=0, \\pm 2, \\pm 4, \\pm 6, \\ldots$\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Main program:**" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from math import sin, cos, exp, pi\n", "import sys\n", "\n", "def g(x):\n", " return exp(-0.1*x**2)*sin(pi/2*x)\n", "\n", "def dg(x):\n", " return -2*0.1*x*exp(-0.1*x**2)*sin(pi/2*x) + \\\n", " pi/2*exp(-0.1*x**2)*cos(pi/2*x)\n", "\n", "x0 = float(sys.argv[1])\n", "x, info = Newton(g, x0, dg, store=True)\n", "print 'Computed zero:', x\n", "\n", "# Print the evolution of the difference equation\n", "# (i.e., the search for the root)\n", "for i in range(len(info)):\n", " print 'Iteration %3d: f(%g)=%g' % (i, info[i][0], info[i][1])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Results from this test problem\n", "\n", "$x_0=1.7$ gives quick convergence towards the closest root $x=0$:" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": false }, "outputs": [], "source": [ "zero: 1.999999999768449\n", "Iteration 0: f(1.7)=0.340044\n", "Iteration 1: f(1.99215)=0.00828786\n", "Iteration 2: f(1.99998)=2.53347e-05\n", "Iteration 3: f(2)=2.43808e-10" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Start value $x_0=3$ (closest root $x=2$ or $x=4$):" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "collapsed": false }, "outputs": [], "source": [ "zero: 42.49723316011362\n", "Iteration 0: f(3)=-0.40657\n", "Iteration 1: f(4.66667)=0.0981146\n", "Iteration 2: f(42.4972)=-2.59037e-79" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## What happened here??\n", "\n", "Try the demo program `src/diffeq/Newton_movie.py` with $x_0=3$,\n", "$x\\in [-2,50]$ for plotting and numerical approximation of $f'(x)$:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " Terminal> python Newton_movie.py \"exp(-0.1*x**2)*sin(pi/2*x)\" \\\n", " numeric 3 -2 50\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Lesson learned:**\n", "\n", "Newton's method may work fine or give wrong results! You need to understand the method to interpret the results!\n", "\n", "\n", "\n", "## First step: we're moving to the right ($x=4$?)\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "## Second step: oops, too much to the right...\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "## Third step: disaster since we're \"done\" ($f(x)\\approx 0$)\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "## Programming with sound\n", "\n", "**Tones are sine waves:**\n", "\n", "A tone **A** (440 Hz) is a sine wave with frequency 440 Hz:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "s(t) = A\\sin\\left( 2\\pi f t\\right),\\quad f = 440\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "On a computer we represent $s(t)$ by a discrete set of points on the function curve (exactly as we do when we plot $s(t)$).\n", "CD quality needs 44100 samples per second.\n", "\n", "\n", "\n", "## Making a sound file with single tone (part 1)\n", "\n", " * $r$: sampling rate (samples per second, default 44100)\n", "\n", " * $f$: frequency of the tone\n", "\n", " * $m$: duration of the tone (seconds)\n", "\n", "Sampled sine function for this tone:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "s_n =\n", "A\\sin\\left( 2\\pi f {n\\over r}\\right),\n", "\\quad n=0,1,\\ldots, m\\cdot r\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Code (we use descriptive names: `frequency` $f$, `length` $m$, `amplitude` $A$, `sample_rate` $r$):" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import numpy\n", "def note(frequency, length, amplitude=1,\n", " sample_rate=44100):\n", " time_points = numpy.linspace(0, length,\n", " length*sample_rate)\n", " data = numpy.sin(2*numpy.pi*frequency*time_points)\n", " data = amplitude*data\n", " return data" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Making a sound file with single tone (part 2)\n", "\n", " * We have `data` as an array with `float` and unit amplitude\n", "\n", " * Sound data in a file should have 2-byte integers (`int16`) as data elements and amplitudes up to $2^{15}-1$ (max value for `int16` data)" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "collapsed": false }, "outputs": [], "source": [ "data = note(440, 2)\n", "data = data.astype(numpy.int16)\n", "max_amplitude = 2**15 - 1\n", "data = max_amplitude*data\n", "import scitools.sound\n", "scitools.sound.write(data, 'Atone.wav')\n", "scitools.sound.play('Atone.wav')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Reading sound from file\n", "\n", " * Let us read a sound file and add echo\n", "\n", " * Sound = array `s[n]`\n", "\n", " * Echo means to add a delay of the sound" ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# echo: e[n] = beta*s[n] + (1-beta)*s[n-b]\n", "\n", "def add_echo(data, beta=0.8, delay=0.002,\n", " sample_rate=44100):\n", " newdata = data.copy()\n", " shift = int(delay*sample_rate) # b (math symbol)\n", " for i in xrange(shift, len(data)):\n", " newdata[i] = beta*data[i] + (1-beta)*data[i-shift]\n", " return newdata" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Load data, add echo and play:" ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "collapsed": false }, "outputs": [], "source": [ "data = scitools.sound.read(filename)\n", "data = data.astype(float)\n", "data = add_echo(data, beta=0.6)\n", "data = data.astype(int16)\n", "scitools.sound.play(data)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Playing many notes\n", "\n", " * Each note is an array of samples from a sine with a frequency corresponding to the note\n", "\n", " * Assume we have several note arrays `data1`, `data2`, ...:" ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# put data1, data2, ... after each other in a new array:\n", "data = numpy.concatenate((data1, data2, data3, ...))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The start of \"Nothing Else Matters\" (Metallica):" ] }, { "cell_type": "code", "execution_count": 22, "metadata": { "collapsed": false }, "outputs": [], "source": [ "E1 = note(164.81, .5)\n", "G = note(392, .5)\n", "B = note(493.88, .5)\n", "E2 = note(659.26, .5)\n", "intro = numpy.concatenate((E1, G, B, E2, B, G))\n", "...\n", "song = numpy.concatenate((intro, intro, ...))\n", "scitools.sound.play(song)\n", "scitools.sound.write(song, 'tmp.wav')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Summary of difference equations\n", "\n", " * Sequence: $x_0$, $x_1$, $x_2$, $\\ldots$, $x_n$, $\\ldots$, $x_N$\n", "\n", " * Difference equation: relation between $x_n$, $x_{n-1}$ and maybe $x_{n-2}$ (or more terms in the \"past\") + known start value $x_0$ (and more values $x_1$, ... if more levels enter the equation)\n", "\n", "Solution of difference equations by simulation:" ] }, { "cell_type": "code", "execution_count": 23, "metadata": { "collapsed": false }, "outputs": [], "source": [ "index_set =