diff --git a/00_intro_numerical_methods-extra-examples.ipynb b/00_intro_numerical_methods-extra-examples.ipynb new file mode 100644 index 00000000..3258e754 --- /dev/null +++ b/00_intro_numerical_methods-extra-examples.ipynb @@ -0,0 +1,318 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": { + "hide_input": false, + "slideshow": { + "slide_type": "skip" + } + }, + "source": [ + "\n", + " \n", + "
\n", + " Text provided under a Creative Commons Attribution license, CC-BY. All code is made available under the FSF-approved MIT license. (c) Kyle T. Mandli
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "hide_input": false, + "slideshow": { + "slide_type": "skip" + } + }, + "outputs": [], + "source": [ + "from __future__ import print_function\n", + "\n", + "%matplotlib inline\n", + "import numpy\n", + "import matplotlib.pyplot as plt" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "skip" + } + }, + "source": [ + "Note to lecturers: This notebook is designed to work best as a classic Jupyter Notebook with nbextensions \n", + "* hide_input: to hide selected python cells particularly for just plotting\n", + "* RISE: Interactive js slide presentations" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "skip" + } + }, + "source": [ + "## Why should I care?" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "skip" + } + }, + "source": [ + "*Because I have to be here for a requirement...*" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "skip" + } + }, + "source": [ + "*or I might actually need to solve something important...like can I ever retire?*" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "skip" + } + }, + "source": [ + "### An example: The Retirement Problem\n", + "\n", + "$$A = \\frac{P}{r} \\left((1+r)^n - 1 \\right)$$\n", + "\n", + "$A$ is the total amount after n payments\n", + "\n", + "$P$ is the incremental payment\n", + "\n", + "$r$ is the interest rate per payment period\n", + "\n", + "$n$ is the number of payments\n", + "\n", + "\n", + "\n", + "Note that these can all be functions of $r$, $n$, and time" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "hide_input": true, + "slideshow": { + "slide_type": "subslide" + } + }, + "source": [ + "What is the minimum interest rate $r$ required to accrue \\$1M over 20 years saving \\\\$1000 per month?\n", + "\n", + "$$\n", + "A = \\frac{P}{r} \\left((1+r)^n - 1 \\right)\n", + "$$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "hide_input": false, + "jupyter": { + "source_hidden": true + }, + "slideshow": { + "slide_type": "skip" + } + }, + "outputs": [], + "source": [ + "def A(P, r, n):\n", + " return P / r * ((1 + r)**n - 1)\n", + "\n", + "P = 1000.\n", + "years = numpy.linspace(0,20,241)\n", + "n = 12*years\n", + "target = 1.e6" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "hide_input": false, + "jupyter": { + "source_hidden": true + }, + "slideshow": { + "slide_type": "skip" + } + }, + "outputs": [], + "source": [ + "# find the root using scipy's brentq method\n", + "from scipy.optimize import brentq\n", + "n_max = n[-1]\n", + "f = lambda r : target - A(P,r/12., n_max)\n", + "r_target = brentq(f,0.1,0.15)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "hide_input": true, + "jupyter": { + "source_hidden": true + }, + "slideshow": { + "slide_type": "fragment" + } + }, + "outputs": [], + "source": [ + "fig = plt.figure(figsize=(8,6))\n", + "fig.set_figwidth(fig.get_figwidth() * 2)\n", + "\n", + "axes = fig.add_subplot(1, 2, 1)\n", + "for r in [0.02, 0.05, 0.08, 0.1, 0.12, 0.15]:\n", + " axes.plot(years, A(P, r/12., n),label='r = {}'.format(r))\n", + "axes.plot(years, numpy.ones(years.shape) * target, 'k--',linewidth=2,label=\"Target\")\n", + "axes.legend(loc='best')\n", + "axes.set_xlabel(\"Years\",fontsize=18)\n", + "axes.set_ylabel(\"Annuity Value (Dollars)\",fontsize=18)\n", + "axes.set_title(\"The Forward Problem $A(P,r,n)$\",fontsize=20)\n", + "axes.grid()\n", + "\n", + "axes = fig.add_subplot(1, 2, 2)\n", + "r = numpy.linspace(1.e-6,0.15,100)\n", + "target = 1.e6\n", + "axes.plot(A(P,r/12.,n_max),r,linewidth=2)\n", + "axes.plot(numpy.ones(r.shape) * target, r,'k--',linewidth=2)\n", + "axes.scatter(A(P,r_target/12,n_max),r_target,s=100,c='r')\n", + "axes.set_xlabel('Annuity Value (Dollars)',fontsize=18)\n", + "axes.set_title('The Inverse Problem $r(A,P,n)$, $r\\geq$ {:.3f}'.format(r_target),fontsize=20)\n", + "axes.grid()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "source": [ + "This is a rootfinding problem for a function of a single variable" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "skip" + } + }, + "source": [ + "### Another Example: Boat race (numerical quadrature of arc length)\n", + "Given a river with coordinates $(x,f(x))$ find the total length actually rowed over a given interval $x\\in [0,X]$\n", + "\n", + "Example: a sinusoidal river $$f(x) = A \\sin x$$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "hide_input": true, + "slideshow": { + "slide_type": "skip" + } + }, + "outputs": [], + "source": [ + "A=.5\n", + "fig = plt.figure(figsize=(8,6))\n", + "x = numpy.linspace(0,10,100)\n", + "plt.plot(x,A*numpy.sin(x),'r',linewidth=2,)\n", + "plt.xlabel('distance (km)')\n", + "plt.grid()\n", + "plt.title('The Boat Race: $f(x)=A\\sin(x)')\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "skip" + } + }, + "source": [ + "We need to calculate the function $f(x)$'s arc-length from $[0, X]$ e.g.\n", + "\n", + "\\begin{align}\n", + " L(X) &= \\int_0^{X} \\sqrt{1 + |f'(x)|^2} dx\\\\\n", + " &= \\int_0^{X} \\sqrt{1 + A^2\\cos^2(x)} dx\\\\\n", + "\\end{align}\n", + "\n", + "In general the value of this integral cannot be solved analytically (the answer is actually an elliptic integral). We usually need to approximate the integral using a numerical quadrature." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "hide_input": true, + "slideshow": { + "slide_type": "skip" + } + }, + "source": [ + "## Why is this exciting?\n", + "\n", + "Computers have revolutionized our ability to explore...\n", + "\n", + "[Top 10 Algorithms of the 20th Century](http://www.siam.org/pdf/news/637.pdf?t=1&cn=ZmxleGlibGVfcmVjcw%3D%3D&refsrc=email&iid=658bdab6af614c83a8df1f8e02035eae&uid=755271476&nid=244+285282312)\n", + "\n", + "...and there is more to come!" + ] + } + ], + "metadata": { + "celltoolbar": "Slideshow", + "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.5" + }, + "latex_envs": { + "bibliofile": "biblio.bib", + "cite_by": "apalike", + "current_citInitial": 1, + "eqLabelWithNumbers": true, + "eqNumInitial": 0 + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/00_intro_numerical_methods.ipynb b/00_intro_numerical_methods.ipynb index 3dfaf964..4a479ccd 100644 --- a/00_intro_numerical_methods.ipynb +++ b/00_intro_numerical_methods.ipynb @@ -3,6 +3,7 @@ { "cell_type": "markdown", "metadata": { + "hide_input": false, "slideshow": { "slide_type": "skip" } @@ -18,7 +19,7 @@ "cell_type": "code", "execution_count": null, "metadata": { - "collapsed": true, + "hide_input": false, "slideshow": { "slide_type": "skip" } @@ -36,11 +37,13 @@ "cell_type": "markdown", "metadata": { "slideshow": { - "slide_type": "slide" + "slide_type": "skip" } }, "source": [ - "# Introduction and Motivation: Modeling and methods for scientific computing" + "Note to lecturers: This notebook is designed to work best as a classic Jupyter Notebook with nbextensions \n", + "* hide_input: to hide selected python cells particularly for just plotting\n", + "* RISE: Interactive js slide presentations" ] }, { @@ -51,80 +54,206 @@ } }, "source": [ - "## Why are we here?\n", + "# Introduction to Numerical Methods\n", "\n", - "It is only possible to solve relatively simple equations. Unfortunately, most nontrivial systems can not be solved algebraically. For example, determining the root of a fifth order polynomial like\n", - "$$x^5 + 3x^2+ 2x + 3 = 0$$\n", - "cannot be solved directly.\n", + "## Outline\n", "\n", - "In general, the root of most nonlinear systems,\n", - "$$\\vec{f}(x,y,z,t) = \\vec{0},$$\n", - "cannot be directly solved." + "* Introduction and Motivation\n", + "* Some examples of Numerical Problems\n", + "* Overview of Course materials\n", + "* Course Logistics\n", + "* **Prerequisites and \"Fluency\"**\n", + "* Overview of the course computational environment (Python, jupyter, git)\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { - "slide_type": "subslide" + "slide_type": "slide" } }, "source": [ - "Sometimes, a problem can be too big. It may not possible to calculate an exact answer in a convenient amount of time. It may be that a system is over constrained and obtaining an answer that matches all requirements is not possible.\n", + "## What are Numerical Methods?\n", "\n", - "![Linear Regression](./images/linear_regression.png)" + "* An extremely broad field -- much broader than this class" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { - "slide_type": "subslide" + "slide_type": "fragment" + } + }, + "source": [ + "* Rough definition: analysis and application of algorithms to allow computers to solve problems in math, science, engineering..." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "source": [ + "* Strictly speaking, numerical methods don't require computers (e.g. Newton's method is ~ 17 century)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "source": [ + "* But computers make things practical...Consider this class an introduction to computational math" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "hide_input": false, + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Why do we need numerical methods?" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "hide_input": false, + "slideshow": { + "slide_type": "fragment" } }, "source": [ - "Actually want an answer..." + "### Some problems have **no closed form solution** \n", + " \n", + " Example: find $f(x)=0$ for the 5th order polynomial\n", + " \n", + "$$\n", + " x^5 + 3x^2+ 2x + 3 = 0\n", + "$$" ] }, { "cell_type": "markdown", "metadata": { + "hide_input": false, "slideshow": { "slide_type": "subslide" } }, "source": [ - "We usually have experience and insights into a system and have a good idea what the solutions should be. At the same time the numerics can help us build our intuition and complement analytical methods. The numerics can also be used to test ideas and extend our understanding of a system. In doing so we can update mathematical descriptions to account for new aspects.\n", + "$$\n", + " p(x) = x^5 + 3x^2+ 2x + 3 = 0\n", + "$$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "hide_input": true, + "slideshow": { + "slide_type": "fragment" + } + }, + "outputs": [], + "source": [ + "p = lambda x: x**5 +3*x**2 +2*x +3\n", "\n", - "The process of describing physical systems using mathematical expressions is often referred to as mathematical modeling. It is a process in which intuition, simulation, and analysis are used to build expressions that mimic general behaviours of physical systems. The general practice is to start with a basic description and refine it over many iterations. \n", + "from scipy.optimize import brentq\n", "\n", - "Each step in the process introduces error. In this course we will focus on the errors associated with numerical approximation. It is important, though, to be able to determine the difference between errors associated with simplifications associated with the model and the errors associated with the approximation of mathematical expressions.\n", + "x = numpy.linspace(-2,2,100)\n", + "x0 = brentq(p,-2,2)\n", "\n", - "
\n", - "\n", - "
" + "plt.figure(figsize=(8,6))\n", + "plt.plot(x,p(x),linewidth=2)\n", + "plt.plot(x,numpy.zeros(x.shape),'k',linewidth=2)\n", + "plt.scatter(x0,p(x0),s=50,c='r')\n", + "plt.grid()\n", + "plt.xlabel('x',fontsize=16)\n", + "plt.ylabel('p(x)',fontsize=16)\n", + "plt.title('Root at x={}'.format(x0),fontsize=18)\n", + "plt.show()" ] }, { "cell_type": "markdown", "metadata": { + "hide_input": false, "slideshow": { - "slide_type": "slide" + "slide_type": "subslide" + } + }, + "source": [ + "The 5th order polynomial $p(x) = x^5 + 3x^2+ 2x + 3 = 0$ actually has 5 complex roots" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "hide_input": true, + "slideshow": { + "slide_type": "fragment" } }, + "outputs": [], "source": [ - "## Why should I care?" + "c = [1,0,0,3,2,3]\n", + "roots = numpy.roots(c)\n", + "fig = plt.figure(figsize=(8,6))\n", + "plt.scatter(numpy.real(roots),numpy.imag(roots),s=50,c='r')\n", + "x=numpy.linspace(-2,2)\n", + "plt.plot(x,numpy.zeros(x.shape),'k',linewidth=2)\n", + "plt.plot(numpy.zeros(x.shape),x,'k',linewidth=2)\n", + "\n", + "plt.xlabel('Re',fontsize=16)\n", + "plt.ylabel('Im',fontsize=16)\n", + "plt.grid()\n", + "plt.title('Complex roots of a 5th order polynomial')\n", + "plt.axis('square')\n", + "plt.xlim([-2,2])\n", + "plt.ylim([-2,2])\n", + "plt.show()" ] }, { "cell_type": "markdown", "metadata": { + "hide_input": false, + "slideshow": { + "slide_type": "fragment" + } + }, + "source": [ + "which require completely different numerical methods to find." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "hide_input": false, "slideshow": { "slide_type": "subslide" } }, "source": [ - "*Because I have to be here for a requirement...*" + "In general, the roots of most nonlinear systems,\n", + " \n", + "$$\n", + " \\mathbf{F}(x,y,z,t) = \\mathbf{0},\n", + "$$\n", + " \n", + "cannot be solved directly." ] }, { @@ -135,43 +264,66 @@ } }, "source": [ - "#### The Retirement Problem\n", "\n", - "$$A = \\frac{P}{r} \\left((1+r)^n - 1 \\right)$$\n", + "### Some problems are too big to be done by hand. \n", "\n", - "$P$ is the incremental payment\n", + "It may not possible to calculate an exact answer in a convenient amount of time. It may be that a system is over constrained and obtaining an answer that matches all requirements is not possible.\n", "\n", - "$r$ is the interest rate per payment period\n", - "\n", - "$n$ is the number of payments\n", + "\n", + " \n", + " \n", + " \n", + "
" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, + "source": [ + "### Some problems are simply not analytic\n", "\n", - "$A$ is the total amount after n payments\n", + "#### Data analysis of digital data\n", "\n", - "Note that these can all be functions of $r$, $n$, and time" + "Finding trends in real data represented without a closed form (analytical form).\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { - "collapsed": true, + "hide_input": true, "slideshow": { - "slide_type": "skip" + "slide_type": "fragment" } }, "outputs": [], "source": [ - "def A(P, r, n):\n", - " return P / r * ((1 + r)**n - 1)\n", - "\n", - "n = numpy.linspace(0, 20, 100)\n", - "target = 5000\n", - "for r in [0.02, 0.05, 0.08, 0.1, 0.12]:\n", - " plt.plot(n, A(100, r, n))\n", - "plt.plot(n, numpy.ones(n.shape) * target, 'k--')\n", - "plt.legend([\"r = 0.02\", \"r = 0.05\", \"r = 0.08\", \"r = 0.1\", \"r = 0.12\", \"Target\"], loc=2)\n", - "plt.xlabel(\"Years\")\n", - "plt.ylabel(\"Annuity Value (Dollars)\")\n", + "data = numpy.loadtxt(\"./data/sunspot.dat\")\n", + "\n", + "fig = plt.figure(figsize=(10, 5))\n", + "fig.set_figwidth(fig.get_figwidth() * 2)\n", + "\n", + "axes = fig.add_subplot(1, 2, 1)\n", + "axes.plot(data[:, 0], data[:, 1],linewidth=2)\n", + "axes.set_xlabel(\"Year\",fontsize=16)\n", + "axes.set_ylabel(\"Number\",fontsize=16)\n", + "axes.set_title(\"Number of Sunspots\",fontsize=18)\n", + "\n", + "axes = fig.add_subplot(1, 2, 2)\n", + "N = int(data.shape[0] / 2)\n", + "period = 1.0 / numpy.fft.fftfreq(data.shape[0])[1:N]\n", + "sunspot_freq = numpy.fft.fft(data[:, 1])[1:N]\n", + "freq = numpy.fft.fftfreq(data.shape[0])[1:N]\n", + "axes.plot(period, numpy.abs(sunspot_freq)**2,linewidth=2)\n", + "axes.set_xlabel(\"Years/Cycle\",fontsize=16)\n", + "axes.set_ylabel(\"Power Spectrum\",fontsize=16)\n", + "axes.set_title(\"Frequency of Sunspots\",fontsize=18)\n", + "axes.set_xlim((0, 50))\n", + "\n", "plt.show()" ] }, @@ -183,25 +335,36 @@ } }, "source": [ - "#### Boat race\n", - "Given a river (say a sinusoid) find the total length actually rowed over a given interval\n", + "### Sometimes you actually want an answer \n", "\n", - "$$f(x) = A \\sin x$$" + "(rather than show it exists or is unique)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { - "slide_type": "-" + "slide_type": "fragment" } }, "source": [ - "We need to calculate the function $f(x)$'s arc-length from $[0, 4 \\pi]$\n", - "\n", - "$$L = \\int_0^{4 \\pi} \\sqrt{1 + |f'(x)|^2} dx$$\n", + "##### Example: Solution of non-linear dynamical systems\n", "\n", - "In general the value of this integral cannot be solved analytically. We usually need to approximate the integral using a numerical quadrature." + "$$\n", + " \\frac{d\\mathbf{u}}{dt} = \\mathbf{F}(t,\\mathbf{u}),\\quad \\mathbf{u}(0)=\\mathbf{u}_0\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "source": [ + "* Can prove that if $\\mathbf{F}$ is sufficiently smooth (Lipshitz Continuous) over some domain that a unique solution exists for some local interval in time (Picard-Lindelhof Theorem)\n", + "* However, the theorem provides **no** way to efficiently find a solution" ] }, { @@ -212,19 +375,107 @@ } }, "source": [ - "#### Non-Linear population growth\n", + "### Example: Non-Linear population growth\n", "\n", - "Lotka-Volterra Predator-Prey model\n", + "Lotka-Volterra Predator-Prey model: Rabbits v. Foxes\n", "\n", "The variable $R$ represents the number of a prey animals in a population. The variable $F$ represents the number of a predators in a population. The interactions between the two can be approximated with the system of differential equations\n", + "\n", "$$\n", "\\begin{align}\n", - " \\frac{d R}{dt} &= R \\cdot (a - b \\cdot F), \\\\\n", - " \\frac{d F}{dt} &= F \\cdot (c \\cdot R + d).\n", + " \\frac{d R}{dt} &= R \\cdot (a - bF)\\\\\n", + " \\frac{d F}{dt} &= F \\cdot (-c + dR), \n", "\\end{align}\n", "$$\n", "\n", - "Questions an ecologist or a mathematician may have:\n", + "where $a,b,c,d$ are parameters of the problem" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "hide_input": false, + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [], + "source": [ + "# The Lotka Volterra predator prey problem using SciPy's \n", + "# ODE integrator solve_ivp\n", + "\n", + "from scipy.integrate import solve_ivp\n", + "\n", + "def lotkavolterra(t, u, a, b, c, d):\n", + " r,f = u\n", + " return [r*(a - b*f), f*(-c + d*r) ]\n", + "\n", + "a,b,c,d = (1.5, 1, 3, 1)\n", + "f = lambda t,u : lotkavolterra(t,u,a,b,c,d)\n", + "\n", + "time_span = [0, 15]\n", + "u_0 = [10, 5]\n", + "sol = solve_ivp(f, time_span , u_0, rtol=1.e-6, atol=1.e-9,dense_output = True)\n" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, + "source": [ + "### Numerical Solutions \n", + "\n", + "$$\n", + " R_0 = 10,\\, F_0 = 5\\quad a,b,c,d = (1.5, 1, 3, 1)\n", + "$$\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "hide_input": true, + "slideshow": { + "slide_type": "fragment" + } + }, + "outputs": [], + "source": [ + "t = numpy.linspace(0, 15, 300)\n", + "z = sol.sol(t)\n", + "\n", + "fig = plt.figure(figsize=(16,6))\n", + "axes = fig.add_subplot(1,2,1)\n", + "axes.plot(t,z[0],'r',label='rabbits')\n", + "axes.plot(t,z[1],'b',label='foxes')\n", + "axes.legend(loc='best',shadow=True)\n", + "axes.set_xlabel('Time',fontsize=16)\n", + "axes.set_ylabel('Population',fontsize=16)\n", + "axes.grid()\n", + "axes.set_title('Lotka-Volterra Predator Prey system',fontsize=18)\n", + "\n", + "axes = fig.add_subplot(1,2,2)\n", + "axes.plot(z[0],z[1],linewidth=1)\n", + "axes.grid()\n", + "axes.set_xlabel('Rabbits',fontsize=16)\n", + "axes.set_ylabel('Foxes',fontsize=16)\n", + "axes.set_title('Lotka-Volterra phase portrait',fontsize=18)\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, + "source": [ + "### Questions an ecologist or a mathematician might ask:\n", " - Where are the steady states?\n", " - Are the solutions to the system stable?\n", " - How do we solve the initial value problem?\n", @@ -241,54 +492,78 @@ } }, "source": [ - "#### Interpolation and Data Fitting\n", + "### Big Points\n", "\n", - "Finding trends in real data represented without a closed form (analytical form).\n", + "* Numerics **complement** analytic methods, but don't replace them\n", + " \n", + "* You need both to **understand** your problems\n", "\n", - "Sunspot counts" + "* Numerical methods are only a part of the bigger problem of **Mathematical Modeling**\n" ] }, { - "cell_type": "code", - "execution_count": null, + "cell_type": "markdown", "metadata": { - "collapsed": true, + "hide_input": true, "slideshow": { - "slide_type": "skip" + "slide_type": "subslide" } }, - "outputs": [], "source": [ - "data = numpy.loadtxt(\"./data/sunspot.dat\")\n", - "\n", - "fig = plt.figure()\n", - "fig.set_figwidth(fig.get_figwidth() * 2)\n", + "## Mathematical Modeling\n", "\n", - "axes = fig.add_subplot(1, 2, 1)\n", - "axes.plot(data[:, 0], data[:, 1])\n", - "axes.set_xlabel(\"Year\")\n", - "axes.set_ylabel(\"Number\")\n", - "axes.set_title(\"Number of Sunspots\")\n", + "\n", + " \n", + " \n", + " \n", + "
\n", + "
\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "source": [ + "* Murphy was a modeler..." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "source": [ + "* Strictly speaking this is not a course in Mathematical Modeling" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "skip" + } + }, + "source": [ + "We usually have experience and insights into a system and have a good idea what the solutions should be. At the same time the numerics can help us build our intuition and complement analytical methods. The numerics can also be used to test ideas and extend our understanding of a system. In doing so we can update mathematical descriptions to account for new aspects.\n", "\n", - "axes = fig.add_subplot(1, 2, 2)\n", - "N = int(data.shape[0] / 2)\n", - "period = 1.0 / numpy.fft.fftfreq(data.shape[0])[1:N]\n", - "sunspot_freq = numpy.fft.fft(data[:, 1])[1:N]\n", - "freq = numpy.fft.fftfreq(data.shape[0])[1:N]\n", - "axes.plot(period, numpy.abs(sunspot_freq)**2)\n", - "axes.set_xlabel(\"Years/Cycle\")\n", - "axes.set_ylabel(\"Power Spectrum\")\n", - "axes.set_title(\"Frequency of Sunspots\")\n", - "axes.set_xlim((0, 50))\n", + "The process of describing physical systems using mathematical expressions is often referred to as mathematical modeling. It is a process in which intuition, simulation, and analysis are used to build expressions that mimic general behaviours of physical systems. The general practice is to start with a basic description and refine it over many iterations. \n", "\n", - "plt.show()" + "Each step in the process introduces error. In this course we will focus on the errors associated with numerical approximation. It is important, though, to be able to determine the difference between errors associated with simplifications associated with the model and the errors associated with the approximation of mathematical expressions." ] }, { "cell_type": "markdown", "metadata": { + "hide_input": true, "slideshow": { - "slide_type": "subslide" + "slide_type": "skip" } }, "source": [ @@ -309,21 +584,267 @@ } }, "source": [ - "## Tools" + "## The issues:\n", + "### Accuracy and Efficiency\n", + "\n", + "- Numerical methods, invariably include an enormous range of approximations , each with attendant errors " + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "source": [ + "- Good numerical methods also return error estimates, and are stable in the presence of floating point error" ] }, { "cell_type": "markdown", "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "source": [ + "- The detailed analysis of algorithms and their errors is formally Numerical Analysis " + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## **This is not a class in numerical analysis** \n", + "\n", + "This is principally a Methods class where I will emphasize\n", + "\n", + "- Standard Methods and their errors\n", + "- Give insight into how they work (and when they don't)\n", + "- Give you practice implementing them to solve **your** problems\n", + "- Point you to high quality numerical software that you should be using in your own work\n", + "- Help you make good choices " + ] + }, + { + "cell_type": "markdown", + "metadata": { + "hide_input": true, "slideshow": { "slide_type": "subslide" } }, "source": [ - "### Computer Languages\n", + "\n", + "\n", + " \n", + "\n", + " \n", + "
\n", "\n" ] }, + { + "cell_type": "markdown", + "metadata": { + "hide_input": true, + "slideshow": { + "slide_type": "subslide" + } + }, + "source": [ + "\n", + "\n", + " \n", + "\n", + " \n", + "
\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Course Content\n", + "\n", + "Topics Covered: (see Syllabus on [Courseworks](https://courseworks2.columbia.edu/courses/95683/assignments/syllabus))\n", + "\n", + "1. Quick *Review* of computational Tools (Python,numpy,matplotlib,scipy,git)\n", + "- Sources of Error and Error Analysis\n", + "- Rootfinding/Optimization of non-linear functions $f(x)$ \n", + "- Interpolation\n", + "- Numerical integration (quadrature) and differentiation\n", + "- Solution of ODE Initial value problems\n", + "- Numerical Linear Algebra\n", + "- Solving systems of non-linear equations $\\mathbf{F}(\\mathbf{x})=\\mathbf{0}$\n", + "- ODE 2-point Boundary value problems (towards numerical PDE's)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, + "source": [ + "### Topics Not Covered: \n", + "\n", + "1. Optimization -- linear programming, constrained optimization\n", + "- Numerical Solution of PDE's (APMA E4301)\n", + "- Mathematical Modeling (maybe a little)\n", + "- Machine Learning/Data Science (APMA 4990 Intro Math of Data Science)\n", + "- Lot's of things that would be fun...\n" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Purpose of this course:\n", + "\n", + "Choose and understand critical methods that are essential for Numerical PDE's, modeling and scientific computation.\n", + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + "
HPC model of Windturbine turbulence by Andrew Kirby\n", + " FEniCS FEM model of a turbocharger\n", + "
" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "hide_input": true, + "slideshow": { + "slide_type": "subslide" + } + }, + "source": [ + "### A FEM model from Spiegelman\n", + "\n", + "Coupled fluid/solid dynamics model of magmatism in the Earth's mantle." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "hide_input": true, + "slideshow": { + "slide_type": "fragment" + } + }, + "outputs": [], + "source": [ + "from IPython.lib.display import YouTubeVideo\n", + "vid = YouTubeVideo('Pfp4snvqhus',rel=0,width=700,height=450, loop=1)\n", + "display(vid)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, + "source": [ + "### All of these problems require Interpolation, numerical quadrature,numerical linear algebra ,\n", + "\n", + "#### Plus PDE's, Vector Calculus, functional analysis,computational geometry, C++/Python programming" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "hide_input": true, + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Course Logistics\n", + "\n", + "- 2 Lectures per week\n", + "- ~1 Homework per week (All jupyter notebook based) (50%)\n", + "- 1 takehome Midterm (25%)\n", + "- 1 takehome Final (25%)\n", + "\n", + "\n", + "- No primary text but\n", + "- All notes available as Jupyter notebooks through github" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Course Pre-requisites\n", + "\n", + "### Numerical methods aren't particularly hard...\n", + "\n", + "* however...They require significant **fluency** in\n", + " * Calculus\n", + " * Vector Calculus\n", + " * ODE's\n", + " * Linear Algebra\n", + " * Programming: Python 3 (numpy, matplotlib, scipy)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, + "source": [ + "#### Computational Fluency \n", + "- the ability to:\n", + " - move smoothly between continous and discrete problems\n", + " - recognize and use a range of mathematical/computational techniques\n", + " - use best practices for software development \n", + " - keep track of errors and artifacts\n", + " - stay on top of complex, multi-part problems\n", + " - **Debug**\n", + " \n", + "- Anything that can go wrong...will go wrong" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Computational Tools/Environment" + ] + }, { "cell_type": "markdown", "metadata": { @@ -331,6 +852,21 @@ "slide_type": "subslide" } }, + "source": [ + "### Computer Languages\n", + "\n", + "* Compiled languages: C, C++, Fortran\n", + "* Commercial interactive modeling systems: Matlab, Mathematical, IDL\n", + "* Open-source interpreted languages: Python, R, Ruby, Julia" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "skip" + } + }, "source": [ "#### C, C++, Fortran\n", "\n", @@ -348,7 +884,7 @@ "cell_type": "markdown", "metadata": { "slideshow": { - "slide_type": "subslide" + "slide_type": "skip" } }, "source": [ @@ -371,7 +907,8 @@ } }, "source": [ - "#### Python\n", + "### Python (3.x)\n", + "\n", "##### Features and Project Goals:\n", " - Python is free (BSD-like license) and highly portable (Windows, Mac OS X, Linux, etc.)\n", " - Interactive interpreter\n", @@ -447,27 +984,25 @@ } }, "source": [ - "### Peer Review\n", + "### Jupyter Notebooks\n", "\n", - "**Why?**\n", - "In this class reviewing your peer's work can lead to\n", - " - Better understanding of a problem\n", - " - See alternative solutions to the same problem\n", - " - Learn how to read other people's code\n", - " - Hopefully learn some tips about your own coding style\n", - " - Another skill highly used in practice\n", - " - Demonstrate why readibility and discipline with respect to code matters." + "The notebook environment gives us a convenient means for working with code while annotating it. We will only cover the key concepts here and hope that you will explore on your own the environments.\n", + "\n", + "[Documentation](https://jupyter.readthedocs.io/en/latest/)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { - "slide_type": "slide" + "slide_type": "subslide" } }, "source": [ - "## Infrastructure" + "#### Toolbar\n", + "\n", + " - Notebooks are modal, they have an edit mode (editing cells) and command mode.\n", + " - Highly recommend taking the tour and checking out the help menu " ] }, { @@ -478,22 +1013,31 @@ } }, "source": [ - "### Python 3.x" + "#### Content types\n", + " - Markdown\n", + " - LaTeX -- $x^2 + y^2 = \\sin(x)$\n", + " - Python\n", + " - NumPy, SciPy, and other packages\n", + " - Other languages (R, Julia)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { - "slide_type": "subslide" + "slide_type": "slide" } }, "source": [ - "### Jupyter Notebooks\n", + "### Running the Notebooks\n", "\n", - "The notebook environment gives us a convenient means for working with code while annotating it. We will only cover the key concepts here and hope that you will explore on your own the environments.\n", + "To run and interact with the notebooks you need the Jupyter software and scipy software stack\n", "\n", - "[Documentation](https://jupyter.readthedocs.io/en/latest/)" + "#### Interactive Environments -- all web-based\n", + "\n", + "* Jupyter Lab\n", + "* Classical Jupyter Notebook\n", + "* Presentation modes (e.g. RISE)" ] }, { @@ -504,25 +1048,32 @@ } }, "source": [ - "#### Toolbar\n", - "\n", - " - Notebooks are modal, they have an edit mode (editing cells) and command mode.\n", - " - Highly recommend taking the tour and checking out the help menu " + "#### Standard Scipy software stack\n", + "* Python 3.x\n", + "* SciPy\n", + "* Numpy\n", + "* Matplotlib\n", + "* Sympy\n", + "* Pandas\n", + "* etc " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { - "slide_type": "subslide" + "slide_type": "slide" } }, "source": [ - "#### Content types\n", - " - Markdown\n", - " - LaTeX $x^2$\n", - " - Python\n", - " - NumPy, SciPy, and other packages" + "### Installations: A few options\n", + "\n", + "#### Install on your own machine \n", + " - [Notebook Quick Start Guide](https://jupyter-notebook-beginner-guide.readthedocs.io/en/latest/index.html)\n", + " - [Anaconda](http://anaconda.com/downloads)\n", + " - [Canopy](https://store.enthought.com/downloads/) \n", + "\n", + "Remember to grab the **Python 3.x** versions!" ] }, { @@ -533,7 +1084,13 @@ } }, "source": [ - "#### Jupyter Lab" + "#### The \"cloud\"\n", + "\n", + "Instead of running things locally on your machine there are a number of cloud services that you are welcome to use in order to get everything running easily:\n", + " - Columbia's CUIT Jupyter Hub server\n", + " - [Google's Colaboratory](https://colab.research.google.com)\n", + " - [Microsoft Azure Notebooks](https://notebooks.azure.com)\n", + " - [CoCalc](https://cocalc.com/)" ] }, { @@ -561,7 +1118,7 @@ "source": [ "**Clone** the repository\n", "\n", - "`$> git clone git://github.com/mandli/intro-numerical-methods`\n", + "`$> git clone git://github.com/mspieg/intro-numerical-methods`\n", "\n", "**Pull** in new changes\n", "\n", @@ -583,21 +1140,6 @@ "Also note that you can watch what changes were made and submit **issues** to the github project page if you find mistakes (PLEASE DO THIS!). It will be much easier to do so if you have your own account on github and work with your own clone. You can then submit a merge request to have your changes incorporated into other repositories." ] }, - { - "cell_type": "markdown", - "metadata": { - "slideshow": { - "slide_type": "slide" - } - }, - "source": [ - "### Installation\n", - "\n", - "A few options\n", - " 1. Install on your own machine\n", - " 2. Use a cloud service" - ] - }, { "cell_type": "markdown", "metadata": { @@ -606,42 +1148,29 @@ } }, "source": [ - "#### Your own machine\n", - "\n", - " - [Notebook Quick Start Guide](https://jupyter-notebook-beginner-guide.readthedocs.io/en/latest/index.html)\n", - " - [Anaconda](http://anaconda.com/downloads)\n", - " - [Canopy](https://store.enthought.com/downloads/) \n", - "\n", - "Remember to grab the 3.x versions!" - ] - }, - { - "cell_type": "markdown", - "metadata": { - "slideshow": { - "slide_type": "subslide" - } - }, - "source": [ - "#### The \"cloud\"\n", + "### BRING YOUR LAPTOP!\n", "\n", - "Instead of running things locally on your machine there are a number of cloud services that you are welcome to use in order to get everything running easily:\n", - " - [Google's Colaboratory](https://colab.research.google.com)\n", - " - [Microsoft Azure Notebooks](https://notebooks.azure.com)\n", - " - [CoCalc](https://cocalc.com/)" + "In class demos are better with your participation!" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { - "slide_type": "subslide" + "slide_type": "skip" } }, "source": [ - "### BRING YOUR LAPTOP!\n", + "### Peer Review\n", "\n", - "In class demos are better with your participation!" + "**Why?**\n", + "In this class reviewing your peer's work can lead to\n", + " - Better understanding of a problem\n", + " - See alternative solutions to the same problem\n", + " - Learn how to read other people's code\n", + " - Hopefully learn some tips about your own coding style\n", + " - Another skill highly used in practice\n", + " - Demonstrate why readibility and discipline with respect to code matters." ] } ], @@ -649,7 +1178,7 @@ "celltoolbar": "Slideshow", "kernelspec": { "display_name": "Python 3", - "language": "python3", + "language": "python", "name": "python3" }, "language_info": { @@ -662,7 +1191,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.6.5" + "version": "3.7.5" }, "latex_envs": { "bibliofile": "biblio.bib", @@ -673,5 +1202,5 @@ } }, "nbformat": 4, - "nbformat_minor": 1 + "nbformat_minor": 4 } diff --git a/images/Fenics_tc_vm.png b/images/Fenics_tc_vm.png new file mode 100644 index 00000000..00e60462 Binary files /dev/null and b/images/Fenics_tc_vm.png differ diff --git a/images/MakeGoodChoices.gif b/images/MakeGoodChoices.gif new file mode 100644 index 00000000..8420f8db Binary files /dev/null and b/images/MakeGoodChoices.gif differ diff --git a/images/RuiningMyLife.gif b/images/RuiningMyLife.gif new file mode 100644 index 00000000..2c8af72e Binary files /dev/null and b/images/RuiningMyLife.gif differ diff --git a/images/siemens-volrender0000_kirby_modeling.png b/images/siemens-volrender0000_kirby_modeling.png new file mode 100644 index 00000000..55ceae44 Binary files /dev/null and b/images/siemens-volrender0000_kirby_modeling.png differ