diff --git a/16_splines/.ipynb_checkpoints/16_splines-checkpoint.ipynb b/16_splines/.ipynb_checkpoints/16_splines-checkpoint.ipynb new file mode 100644 index 0000000..4f7a14f --- /dev/null +++ b/16_splines/.ipynb_checkpoints/16_splines-checkpoint.ipynb @@ -0,0 +1,1439 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 69, + "metadata": { + "collapsed": true, + "slideshow": { + "slide_type": "skip" + } + }, + "outputs": [], + "source": [ + "setdefaults" + ] + }, + { + "cell_type": "code", + "execution_count": 70, + "metadata": { + "collapsed": true, + "slideshow": { + "slide_type": "skip" + } + }, + "outputs": [], + "source": [ + "%plot --format svg" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "# Splines" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, + "source": [ + "\n", + "\n", + "Following interpolation discussion, instead of estimating 9 data points with an eighth-order polynomial, it makes more sense to fit sections of the curve to lower-order polynomials:\n", + "\n", + "0. zeroth-order (nearest neighbor)\n", + "\n", + "1. first-order (linear interpolation)\n", + "\n", + "2. third-order (cubic interpolation)\n", + "\n", + "Matlab and Octave have built-in functions for 1D and 2D interpolation:\n", + "\n", + "`interp1`\n", + "\n", + "`interp2`" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "'interp1' is a function from the file /usr/share/octave/4.0.0/m/general/interp1.m\n", + "\n", + " -- Function File: YI = interp1 (X, Y, XI)\n", + " -- Function File: YI = interp1 (Y, XI)\n", + " -- Function File: YI = interp1 (..., METHOD)\n", + " -- Function File: YI = interp1 (..., EXTRAP)\n", + " -- Function File: YI = interp1 (..., \"left\")\n", + " -- Function File: YI = interp1 (..., \"right\")\n", + " -- Function File: PP = interp1 (..., \"pp\")\n", + "\n", + " One-dimensional interpolation.\n", + "\n", + " Interpolate input data to determine the value of YI at the points\n", + " XI. If not specified, X is taken to be the indices of Y ('1:length\n", + " (Y)'). If Y is a matrix or an N-dimensional array, the\n", + " interpolation is performed on each column of Y.\n", + "\n", + " The interpolation METHOD is one of:\n", + "\n", + " \"nearest\"\n", + " Return the nearest neighbor.\n", + "\n", + " \"previous\"\n", + " Return the previous neighbor.\n", + "\n", + " \"next\"\n", + " Return the next neighbor.\n", + "\n", + " \"linear\" (default)\n", + " Linear interpolation from nearest neighbors.\n", + "\n", + " \"pchip\"\n", + " Piecewise cubic Hermite interpolating\n", + " polynomial--shape-preserving interpolation with smooth first\n", + " derivative.\n", + "\n", + " \"cubic\"\n", + " Cubic interpolation (same as \"pchip\").\n", + "\n", + " \"spline\"\n", + " Cubic spline interpolation--smooth first and second\n", + " derivatives throughout the curve.\n", + "\n", + " Adding '*' to the start of any method above forces 'interp1' to\n", + " assume that X is uniformly spaced, and only 'X(1)' and 'X(2)' are\n", + " referenced. This is usually faster, and is never slower. The\n", + " default method is \"linear\".\n", + "\n", + " If EXTRAP is the string \"extrap\", then extrapolate values beyond\n", + " the endpoints using the current METHOD. If EXTRAP is a number,\n", + " then replace values beyond the endpoints with that number. When\n", + " unspecified, EXTRAP defaults to 'NA'.\n", + "\n", + " If the string argument \"pp\" is specified, then XI should not be\n", + " supplied and 'interp1' returns a piecewise polynomial object. This\n", + " object can later be used with 'ppval' to evaluate the\n", + " interpolation. There is an equivalence, such that 'ppval (interp1\n", + " (X, Y, METHOD, \"pp\"), XI) == interp1 (X, Y, XI, METHOD, \"extrap\")'.\n", + "\n", + " Duplicate points in X specify a discontinuous interpolant. There\n", + " may be at most 2 consecutive points with the same value. If X is\n", + " increasing, the default discontinuous interpolant is\n", + " right-continuous. If X is decreasing, the default discontinuous\n", + " interpolant is left-continuous. The continuity condition of the\n", + " interpolant may be specified by using the options \"left\" or \"right\"\n", + " to select a left-continuous or right-continuous interpolant,\n", + " respectively. Discontinuous interpolation is only allowed for\n", + " \"nearest\" and \"linear\" methods; in all other cases, the X-values\n", + " must be unique.\n", + "\n", + " An example of the use of 'interp1' is\n", + "\n", + " xf = [0:0.05:10];\n", + " yf = sin (2*pi*xf/5);\n", + " xp = [0:10];\n", + " yp = sin (2*pi*xp/5);\n", + " lin = interp1 (xp, yp, xf);\n", + " near = interp1 (xp, yp, xf, \"nearest\");\n", + " pch = interp1 (xp, yp, xf, \"pchip\");\n", + " spl = interp1 (xp, yp, xf, \"spline\");\n", + " plot (xf,yf,\"r\", xf,near,\"g\", xf,lin,\"b\", xf,pch,\"c\", xf,spl,\"m\",\n", + " xp,yp,\"r*\");\n", + " legend (\"original\", \"nearest\", \"linear\", \"pchip\", \"spline\");\n", + "\n", + " See also: pchip, spline, interpft, interp2, interp3, interpn.\n", + "\n", + "Additional help for built-in functions and operators is\n", + "available in the online version of the manual. Use the command\n", + "'doc ' to search the manual index.\n", + "\n", + "Help and information about Octave is also available on the WWW\n", + "at http://www.octave.org and via the help@octave.org\n", + "mailing list.\n" + ] + } + ], + "source": [ + "help interp1" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Choosing " + ] + }, + { + "cell_type": "code", + "execution_count": 75, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "fragment" + } + }, + "outputs": [ + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "Gnuplot\n", + "Produced by GNUPLOT 5.0 patchlevel 3 \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\t\n", + "\t\t-1.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-1\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-0.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t1\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t1.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-3\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-2\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-1\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t1\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t2\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t3\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\tdata\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\tdata\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\n", + "\t\n", + "\tlinear\n", + "\n", + "\t\n", + "\t\tlinear\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tcubic spline\n", + "\n", + "\t\n", + "\t\tcubic spline\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tpiecewise cubic\n", + "\n", + "\t\n", + "\t\tpiecewise cubic\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "x=linspace(-pi,pi,9);\n", + "xi=linspace(-pi,pi,100);\n", + "y=sin(x);\n", + "yi_lin=interp1(x,y,xi,'linear');\n", + "yi_spline=interp1(x,y,xi,'spline'); \n", + "yi_cubic=interp1(x,y,xi,'cubic');\n", + "plot(x,y,'o',xi,yi_lin,xi,yi_spline,xi,yi_cubic)\n", + "axis([-pi,pi,-1.5,1.5])\n", + "legend('data','linear','cubic spline','piecewise cubic','Location','NorthWest')\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Example: Accelerate then hold velocity\n", + "\n", + "Test driving a car, the accelerator is pressed, then released, then pressed again for 20-second intervals, until speed is 120 mph. Here the time is given as vector t in seconds and the velocity is in mph. " + ] + }, + { + "cell_type": "code", + "execution_count": 82, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "Gnuplot\n", + "Produced by GNUPLOT 5.0 patchlevel 3 \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t20\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t40\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t60\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t80\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t100\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t120\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t140\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t20\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t40\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t60\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t80\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t100\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t120\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\n", + "\t\n", + "\t\tv (mph)\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\tt (s)\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\tdata\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\tdata\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\n", + "\t\n", + "\tremoved data point\n", + "\n", + "\t\n", + "\t\tremoved data point\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\n", + "\t\n", + "\tlinear\n", + "\n", + "\t\n", + "\t\tlinear\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tcubic spline\n", + "\n", + "\t\n", + "\t\tcubic spline\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tpiecewise cubic\n", + "\n", + "\t\n", + "\t\tpiecewise cubic\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "t=[0 20 40 68 80 84 96 104 110]';\n", + "v=[0 20 20 80 80 100 100 125 125]';\n", + "tt=linspace(0,110)';\n", + "v_lin=interp1(t,v,tt);\n", + "v_spl=interp1(t,v,tt,'spline');\n", + "v_cub=interp1(t,v,tt,'cubic');\n", + "\n", + "plot(t,v,'o',56,38,'s',tt,v_lin,tt,v_spl,tt,v_cub)\n", + "xlabel('t (s)')\n", + "ylabel('v (mph)')\n", + "legend('data','removed data point','linear','cubic spline','piecewise cubic','Location','NorthWest')" + ] + }, + { + "cell_type": "code", + "execution_count": 83, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "Gnuplot\n", + "Produced by GNUPLOT 5.0 patchlevel 3 \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\t\n", + "\t\t-4\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-2\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t2\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t4\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t6\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t8\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t20\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t40\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t60\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t80\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t100\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t120\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\n", + "\t\n", + "\t\tdv/dt (mph/s)\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\tt (s)\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\tlinear\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\tlinear\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tcubic spline\n", + "\n", + "\t\n", + "\t\tcubic spline\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tpiecewise cubic\n", + "\n", + "\t\n", + "\t\tpiecewise cubic\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "t=[0 20 40 56 68 80 84 96 104 110]';\n", + "v=[0 20 20 38 80 80 100 100 125 125]';\n", + "tt=linspace(0,110)';\n", + "v_lin=interp1(t,v,tt);\n", + "v_spl=interp1(t,v,tt,'spline');\n", + "v_cub=interp1(t,v,tt,'cubic');\n", + "\n", + "\n", + "plot(tt(2:end),diff(v_lin)./diff(tt),tt(2:end),diff(v_spl)./diff(tt),tt(2:end),diff(v_cub)./diff(tt))\n", + "xlabel('t (s)')\n", + "ylabel('dv/dt (mph/s)')\n", + "legend('linear','cubic spline','piecewise cubic','Location','NorthWest')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Choose spline wisely\n", + "\n", + "For example of sin(x), not very important\n", + "\n", + "For stop-and-hold examples, the $C^{2}$-continuity should not be preserved. You don't need smooth curves.\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Numerical Integration\n", + "\n", + "A definite integral is defined by \n", + "\n", + "$I=\\int_{a}^{b}f(x)dx$\n", + "\n", + "To determine the mass of an object with varying density, you can perform a summation\n", + "\n", + "mass=$\\sum_{i=1}^{n}\\rho_{i}\\Delta V_{i}$\n", + "\n", + "or taking the limit as $\\Delta V \\rightarrow dV=dxdydz$\n", + "\n", + "mass=$\\int_{0}^{h}\\int_{0}^{w}\\int_{0}^{l}\\rho(x,y,z)dxdydz$\n", + "\n", + "## Newton-Cotes Formulas\n", + "\n", + "$I=\\int_{a}^{b}f(x)dx=\\int_{a}^{b}f_{n}(x)dx$\n", + "\n", + "where $f_{n}$ is an n$^{th}$-order polynomial approximation of f(x)\n", + "\n", + "## First-Order: Trapezoidal Rule\n", + "\n", + "$I=\\int_{a}^{b}f(x)dx\\approx \\int_{a}^{b}\\left(f(a)+\\frac{f(b)-f(a)}{b-a}(x-a)\\right)dx$\n", + "\n", + "$I\\approx(b-a)\\frac{f(a)+f(b)}{2}$" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "I_trap = 0.78540\n", + "I_act = 1.00000\n" + ] + }, + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "Gnuplot\n", + "Produced by GNUPLOT 5.0 patchlevel 3 \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.2\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.4\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.6\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.8\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t1\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t1\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t1.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t2\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t2.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t3\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t3.5\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\n", + "\n", + "\tgnuplot_plot_1a\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tgnuplot_plot_2a\n", + "\n", + "\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "x=linspace(0,pi)';\n", + "plot(x,sin(x),[0,pi/2],sin([0,pi/2]))\n", + "I_trap=mean(sin(([0,pi/2]))*(diff([0,pi/2])))\n", + "I_act = -(cos(pi/2)-cos(0))\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Improve estimate with more points\n", + "\n", + "$I=\\int_{a}^{b}f(x)dx=\\int_{a}^{a+\\Delta x}f(x)dx+\\int_{a+\\Delta x}^{a+2\\Delta x}f(x)dx+ \\cdots \\int_{b-\\Delta x}^{b}f(x)dx$\n", + "\n", + "$I\\approx\\Delta x\\frac{f(a)+f(a+\\Delta x)}{2}+\\Delta x\\frac{f(a+\\Delta x)+f(a+2\\Delta x)}{2}\n", + "+\\cdots \\Delta x\\frac{f(b-\\Delta x)+f(b)}{2}$\n", + "\n", + "$I\\approx \\frac{\\Delta x}{2}\\left(f(a)+2\\sum_{i=1}^{n-1}f(a+i\\Delta x) +f(b)\\right)$" + ] + }, + { + "cell_type": "code", + "execution_count": 86, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "For 70 steps\n", + "trapezoid approximation of integral is 0.99 \n", + " actual integral is 1.00\n" + ] + } + ], + "source": [ + "N=70;\n", + "I_trap=trap(@(x) sin(x),0,pi/2,N);\n", + "fprintf('For %i steps\\ntrapezoid approximation of integral is %1.2f \\n actual integral is %1.2f',N,I_trap,I_act)" + ] + }, + { + "cell_type": "code", + "execution_count": 92, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "Gnuplot\n", + "Produced by GNUPLOT 5.0 patchlevel 3 \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.2\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.4\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.6\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.8\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t1\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t1\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t1.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t2\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t2.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t3\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t3.5\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\n", + "\n", + "\tgnuplot_plot_1a\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tgnuplot_plot_2a\n", + "\n", + "\t\t \n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "x=linspace(0,pi);\n", + "plot(x,sin(x),linspace(0,pi/2,N),sin(linspace(0,pi/2,N)),'-o')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Increase accuracy = Increase polynomial order\n", + "\n", + "### Simpson's Rules\n", + "\n", + "When integrating f(x) and using a second order polynomial, this is known as **Simpson's 1/3 Rule**\n", + "\n", + "$I=\\frac{h}{3}(f(x_{0})+4f(x_{1})+f(x_{2}))$\n", + "\n", + "where a=$x_{0}$, b=$x_{2}$, and $x_{1}=\\frac{a+b}{2}$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This can be used with n=3 or multiples of 2 intervals\n", + "\n", + "$I=\\int_{x_{0}}^{x_{2}}f(x)dx+\\int_{x_{2}}^{x_{4}}f(x)dx+\\cdots +\\int_{x_{n-2}}^{x_{n}}f(x)dx$\n", + "\n", + "$I=(b-a)\\frac{f(x_{0})+4\\sum_{i=1,3,5}^{n-1}f(x_{i})+2\\sum_{i=2,4,6}^{n-2}f(x_{i})+f(x_{n})}{3n}$" + ] + }, + { + "cell_type": "code", + "execution_count": 93, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "ans = 1.6235\n", + "Is_1_3 = 1.0023\n" + ] + } + ], + "source": [ + "f=@(x) 0.2+25*x-200*x.^2+675*x.^3-900*x.^4+400*x.^5;\n", + "simpson3(f,0,0.8,4)\n", + "Is_1_3=simpson3(@(x) sin(x),0,pi/2,2)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## General Newton-Cotes formulae\n", + "\n", + "![Newton-Cotes Table](newton_cotes.png)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] + } + ], + "metadata": { + "celltoolbar": "Slideshow", + "kernelspec": { + "display_name": "Octave", + "language": "octave", + "name": "octave" + }, + "language_info": { + "file_extension": ".m", + "help_links": [ + { + "text": "MetaKernel Magics", + "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md" + } + ], + "mimetype": "text/x-octave", + "name": "octave", + "version": "0.19.14" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/16_splines/.ipynb_checkpoints/lecture 19-checkpoint.ipynb b/16_splines/.ipynb_checkpoints/lecture 19-checkpoint.ipynb new file mode 100644 index 0000000..2911efe --- /dev/null +++ b/16_splines/.ipynb_checkpoints/lecture 19-checkpoint.ipynb @@ -0,0 +1,497 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "setdefaults" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "%plot --format svg" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Splines (Brief introduction before next section)\n", + "\n", + "Following interpolation discussion, instead of estimating 9 data points with an eighth-order polynomial, it makes more sense to fit sections of the curve to lower-order polynomials:\n", + "\n", + "0. zeroth-order (nearest neighbor)\n", + "\n", + "1. first-order (linear interpolation)\n", + "\n", + "2. third-order (cubic interpolation)\n", + "\n", + "Matlab and Octave have built-in functions for 1D and 2D interpolation:\n", + "\n", + "`interp1`\n", + "\n", + "`interp2`" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "'interp1' is a function from the file /usr/share/octave/4.0.0/m/general/interp1.m\n", + "\n", + " -- Function File: YI = interp1 (X, Y, XI)\n", + " -- Function File: YI = interp1 (Y, XI)\n", + " -- Function File: YI = interp1 (..., METHOD)\n", + " -- Function File: YI = interp1 (..., EXTRAP)\n", + " -- Function File: YI = interp1 (..., \"left\")\n", + " -- Function File: YI = interp1 (..., \"right\")\n", + " -- Function File: PP = interp1 (..., \"pp\")\n", + "\n", + " One-dimensional interpolation.\n", + "\n", + " Interpolate input data to determine the value of YI at the points\n", + " XI. If not specified, X is taken to be the indices of Y ('1:length\n", + " (Y)'). If Y is a matrix or an N-dimensional array, the\n", + " interpolation is performed on each column of Y.\n", + "\n", + " The interpolation METHOD is one of:\n", + "\n", + " \"nearest\"\n", + " Return the nearest neighbor.\n", + "\n", + " \"previous\"\n", + " Return the previous neighbor.\n", + "\n", + " \"next\"\n", + " Return the next neighbor.\n", + "\n", + " \"linear\" (default)\n", + " Linear interpolation from nearest neighbors.\n", + "\n", + " \"pchip\"\n", + " Piecewise cubic Hermite interpolating\n", + " polynomial--shape-preserving interpolation with smooth first\n", + " derivative.\n", + "\n", + " \"cubic\"\n", + " Cubic interpolation (same as \"pchip\").\n", + "\n", + " \"spline\"\n", + " Cubic spline interpolation--smooth first and second\n", + " derivatives throughout the curve.\n", + "\n", + " Adding '*' to the start of any method above forces 'interp1' to\n", + " assume that X is uniformly spaced, and only 'X(1)' and 'X(2)' are\n", + " referenced. This is usually faster, and is never slower. The\n", + " default method is \"linear\".\n", + "\n", + " If EXTRAP is the string \"extrap\", then extrapolate values beyond\n", + " the endpoints using the current METHOD. If EXTRAP is a number,\n", + " then replace values beyond the endpoints with that number. When\n", + " unspecified, EXTRAP defaults to 'NA'.\n", + "\n", + " If the string argument \"pp\" is specified, then XI should not be\n", + " supplied and 'interp1' returns a piecewise polynomial object. This\n", + " object can later be used with 'ppval' to evaluate the\n", + " interpolation. There is an equivalence, such that 'ppval (interp1\n", + " (X, Y, METHOD, \"pp\"), XI) == interp1 (X, Y, XI, METHOD, \"extrap\")'.\n", + "\n", + " Duplicate points in X specify a discontinuous interpolant. There\n", + " may be at most 2 consecutive points with the same value. If X is\n", + " increasing, the default discontinuous interpolant is\n", + " right-continuous. If X is decreasing, the default discontinuous\n", + " interpolant is left-continuous. The continuity condition of the\n", + " interpolant may be specified by using the options \"left\" or \"right\"\n", + " to select a left-continuous or right-continuous interpolant,\n", + " respectively. Discontinuous interpolation is only allowed for\n", + " \"nearest\" and \"linear\" methods; in all other cases, the X-values\n", + " must be unique.\n", + "\n", + " An example of the use of 'interp1' is\n", + "\n", + " xf = [0:0.05:10];\n", + " yf = sin (2*pi*xf/5);\n", + " xp = [0:10];\n", + " yp = sin (2*pi*xp/5);\n", + " lin = interp1 (xp, yp, xf);\n", + " near = interp1 (xp, yp, xf, \"nearest\");\n", + " pch = interp1 (xp, yp, xf, \"pchip\");\n", + " spl = interp1 (xp, yp, xf, \"spline\");\n", + " plot (xf,yf,\"r\", xf,near,\"g\", xf,lin,\"b\", xf,pch,\"c\", xf,spl,\"m\",\n", + " xp,yp,\"r*\");\n", + " legend (\"original\", \"nearest\", \"linear\", \"pchip\", \"spline\");\n", + "\n", + " See also: pchip, spline, interpft, interp2, interp3, interpn.\n", + "\n", + "Additional help for built-in functions and operators is\n", + "available in the online version of the manual. Use the command\n", + "'doc ' to search the manual index.\n", + "\n", + "Help and information about Octave is also available on the WWW\n", + "at http://www.octave.org and via the help@octave.org\n", + "mailing list.\n" + ] + } + ], + "source": [ + "help interp1" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "'interp2' is a function from the file /usr/share/octave/4.0.0/m/general/interp2.m\n", + "\n", + " -- Function File: ZI = interp2 (X, Y, Z, XI, YI)\n", + " -- Function File: ZI = interp2 (Z, XI, YI)\n", + " -- Function File: ZI = interp2 (Z, N)\n", + " -- Function File: ZI = interp2 (Z)\n", + " -- Function File: ZI = interp2 (..., METHOD)\n", + " -- Function File: ZI = interp2 (..., METHOD, EXTRAP)\n", + "\n", + " Two-dimensional interpolation.\n", + "\n", + " Interpolate reference data X, Y, Z to determine ZI at the\n", + " coordinates XI, YI. The reference data X, Y can be matrices, as\n", + " returned by 'meshgrid', in which case the sizes of X, Y, and Z must\n", + " be equal. If X, Y are vectors describing a grid then 'length (X)\n", + " == columns (Z)' and 'length (Y) == rows (Z)'. In either case the\n", + " input data must be strictly monotonic.\n", + "\n", + " If called without X, Y, and just a single reference data matrix Z,\n", + " the 2-D region 'X = 1:columns (Z), Y = 1:rows (Z)' is assumed.\n", + " This saves memory if the grid is regular and the distance between\n", + " points is not important.\n", + "\n", + " If called with a single reference data matrix Z and a refinement\n", + " value N, then perform interpolation over a grid where each original\n", + " interval has been recursively subdivided N times. This results in\n", + " '2^N-1' additional points for every interval in the original grid.\n", + " If N is omitted a value of 1 is used. As an example, the interval\n", + " [0,1] with 'N==2' results in a refined interval with points at [0,\n", + " 1/4, 1/2, 3/4, 1].\n", + "\n", + " The interpolation METHOD is one of:\n", + "\n", + " \"nearest\"\n", + " Return the nearest neighbor.\n", + "\n", + " \"linear\" (default)\n", + " Linear interpolation from nearest neighbors.\n", + "\n", + " \"pchip\"\n", + " Piecewise cubic Hermite interpolating\n", + " polynomial--shape-preserving interpolation with smooth first\n", + " derivative.\n", + "\n", + " \"cubic\"\n", + " Cubic interpolation (same as \"pchip\").\n", + "\n", + " \"spline\"\n", + " Cubic spline interpolation--smooth first and second\n", + " derivatives throughout the curve.\n", + "\n", + " EXTRAP is a scalar number. It replaces values beyond the endpoints\n", + " with EXTRAP. Note that if EXTRAPVAL is used, METHOD must be\n", + " specified as well. If EXTRAP is omitted and the METHOD is\n", + " \"spline\", then the extrapolated values of the \"spline\" are used.\n", + " Otherwise the default EXTRAP value for any other METHOD is \"NA\".\n", + "\n", + " See also: interp1, interp3, interpn, meshgrid.\n", + "\n", + "Additional help for built-in functions and operators is\n", + "available in the online version of the manual. Use the command\n", + "'doc ' to search the manual index.\n", + "\n", + "Help and information about Octave is also available on the WWW\n", + "at http://www.octave.org and via the help@octave.org\n", + "mailing list.\n" + ] + } + ], + "source": [ + "help interp2" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "Gnuplot\n", + "Produced by GNUPLOT 5.0 patchlevel 3 \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\t\n", + "\t\t-1.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-1\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-0.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t1\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t1.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-3\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-2\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-1\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t1\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t2\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t3\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\tdata\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\tdata\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\n", + "\t\n", + "\tlinear\n", + "\n", + "\t\n", + "\t\tlinear\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tcubic spline\n", + "\n", + "\t\n", + "\t\tcubic spline\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tpiecewise cubic\n", + "\n", + "\t\n", + "\t\tpiecewise cubic\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "x=linspace(-pi,pi,9);\n", + "xi=linspace(-pi,pi,100);\n", + "y=sin(x);\n", + "yi_lin=interp1(x,y,xi,'linear');\n", + "yi_spline=interp1(x,y,xi,'spline'); \n", + "yi_cubic=interp1(x,y,xi,'cubic');\n", + "plot(x,y,'o',xi,yi_lin,xi,yi_spline,xi,yi_cubic)\n", + "axis([-pi,pi,-1.5,1.5])\n", + "legend('data','linear','cubic spline','piecewise cubic','Location','NorthWest')\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Example: Accelerate then hold velocity\n", + "\n", + "Here the time is given as vector t in seconds and the " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "t=[0 2 40 56 68 80 84 96 104 110]';\n", + "v=[0 20 20 38 80 80 100 100 125 125]';\n" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Octave", + "language": "octave", + "name": "octave" + }, + "language_info": { + "file_extension": ".m", + "help_links": [ + { + "text": "MetaKernel Magics", + "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md" + } + ], + "mimetype": "text/x-octave", + "name": "octave", + "version": "0.19.14" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/16_splines/16_splines.ipynb b/16_splines/16_splines.ipynb new file mode 100644 index 0000000..1cbc4bb --- /dev/null +++ b/16_splines/16_splines.ipynb @@ -0,0 +1,928 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": true, + "slideshow": { + "slide_type": "skip" + } + }, + "outputs": [], + "source": [ + "setdefaults" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "collapsed": true, + "slideshow": { + "slide_type": "skip" + } + }, + "outputs": [], + "source": [ + "%plot --format svg" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "# Splines" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, + "source": [ + "\n", + "\n", + "Splines fit sections of the curve to lower-order polynomials:\n", + "\n", + "0. zeroth-order (nearest neighbor)\n", + "\n", + "1. first-order (linear interpolation)\n", + "\n", + "2. third-order (cubic interpolation)\n", + "\n", + "Matlab and Octave have built-in functions for 1D and 2D interpolation:\n", + "\n", + "`interp1`\n", + "\n", + "`interp2`" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "'interp1' is a function from the file /usr/share/octave/4.0.0/m/general/interp1.m\n", + "\n", + " -- Function File: YI = interp1 (X, Y, XI)\n", + " -- Function File: YI = interp1 (Y, XI)\n", + " -- Function File: YI = interp1 (..., METHOD)\n", + " -- Function File: YI = interp1 (..., EXTRAP)\n", + " -- Function File: YI = interp1 (..., \"left\")\n", + " -- Function File: YI = interp1 (..., \"right\")\n", + " -- Function File: PP = interp1 (..., \"pp\")\n", + "\n", + " One-dimensional interpolation.\n", + "\n", + " Interpolate input data to determine the value of YI at the points\n", + " XI. If not specified, X is taken to be the indices of Y ('1:length\n", + " (Y)'). If Y is a matrix or an N-dimensional array, the\n", + " interpolation is performed on each column of Y.\n", + "\n", + " The interpolation METHOD is one of:\n", + "\n", + " \"nearest\"\n", + " Return the nearest neighbor.\n", + "\n", + " \"previous\"\n", + " Return the previous neighbor.\n", + "\n", + " \"next\"\n", + " Return the next neighbor.\n", + "\n", + " \"linear\" (default)\n", + " Linear interpolation from nearest neighbors.\n", + "\n", + " \"pchip\"\n", + " Piecewise cubic Hermite interpolating\n", + " polynomial--shape-preserving interpolation with smooth first\n", + " derivative.\n", + "\n", + " \"cubic\"\n", + " Cubic interpolation (same as \"pchip\").\n", + "\n", + " \"spline\"\n", + " Cubic spline interpolation--smooth first and second\n", + " derivatives throughout the curve.\n", + "\n", + " Adding '*' to the start of any method above forces 'interp1' to\n", + " assume that X is uniformly spaced, and only 'X(1)' and 'X(2)' are\n", + " referenced. This is usually faster, and is never slower. The\n", + " default method is \"linear\".\n", + "\n", + " If EXTRAP is the string \"extrap\", then extrapolate values beyond\n", + " the endpoints using the current METHOD. If EXTRAP is a number,\n", + " then replace values beyond the endpoints with that number. When\n", + " unspecified, EXTRAP defaults to 'NA'.\n", + "\n", + " If the string argument \"pp\" is specified, then XI should not be\n", + " supplied and 'interp1' returns a piecewise polynomial object. This\n", + " object can later be used with 'ppval' to evaluate the\n", + " interpolation. There is an equivalence, such that 'ppval (interp1\n", + " (X, Y, METHOD, \"pp\"), XI) == interp1 (X, Y, XI, METHOD, \"extrap\")'.\n", + "\n", + " Duplicate points in X specify a discontinuous interpolant. There\n", + " may be at most 2 consecutive points with the same value. If X is\n", + " increasing, the default discontinuous interpolant is\n", + " right-continuous. If X is decreasing, the default discontinuous\n", + " interpolant is left-continuous. The continuity condition of the\n", + " interpolant may be specified by using the options \"left\" or \"right\"\n", + " to select a left-continuous or right-continuous interpolant,\n", + " respectively. Discontinuous interpolation is only allowed for\n", + " \"nearest\" and \"linear\" methods; in all other cases, the X-values\n", + " must be unique.\n", + "\n", + " An example of the use of 'interp1' is\n", + "\n", + " xf = [0:0.05:10];\n", + " yf = sin (2*pi*xf/5);\n", + " xp = [0:10];\n", + " yp = sin (2*pi*xp/5);\n", + " lin = interp1 (xp, yp, xf);\n", + " near = interp1 (xp, yp, xf, \"nearest\");\n", + " pch = interp1 (xp, yp, xf, \"pchip\");\n", + " spl = interp1 (xp, yp, xf, \"spline\");\n", + " plot (xf,yf,\"r\", xf,near,\"g\", xf,lin,\"b\", xf,pch,\"c\", xf,spl,\"m\",\n", + " xp,yp,\"r*\");\n", + " legend (\"original\", \"nearest\", \"linear\", \"pchip\", \"spline\");\n", + "\n", + " See also: pchip, spline, interpft, interp2, interp3, interpn.\n", + "\n", + "Additional help for built-in functions and operators is\n", + "available in the online version of the manual. Use the command\n", + "'doc ' to search the manual index.\n", + "\n", + "Help and information about Octave is also available on the WWW\n", + "at http://www.octave.org and via the help@octave.org\n", + "mailing list.\n" + ] + } + ], + "source": [ + "help interp1" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "### Choosing the \"right\" spline\n", + "\n", + "Represent sin(x) with a few points and a spline" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "fragment" + } + }, + "outputs": [ + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "Gnuplot\n", + "Produced by GNUPLOT 5.0 patchlevel 3 \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\t\n", + "\t\t-1.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-1\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-0.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t1\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t1.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-3\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-2\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-1\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t1\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t2\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t3\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\tdata\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\tdata\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\n", + "\t\n", + "\tlinear\n", + "\n", + "\t\n", + "\t\tlinear\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tcubic spline\n", + "\n", + "\t\n", + "\t\tcubic spline\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tpiecewise cubic\n", + "\n", + "\t\n", + "\t\tpiecewise cubic\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "x=linspace(-pi,pi,9);\n", + "xi=linspace(-pi,pi,100);\n", + "y=sin(x);\n", + "yi_lin=interp1(x,y,xi,'linear');\n", + "yi_spline=interp1(x,y,xi,'spline'); \n", + "yi_cubic=interp1(x,y,xi,'cubic');\n", + "plot(x,y,'o',xi,yi_lin,xi,yi_spline,xi,yi_cubic)\n", + "axis([-pi,pi,-1.5,1.5])\n", + "legend('data','linear','cubic spline','piecewise cubic','Location','NorthWest')\n" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "### Example: Accelerate then hold velocity\n", + "\n", + "Test driving a car, the accelerator is pressed, then released, then pressed again for 20-second intervals, until speed is 120 mph. Here the time is given as vector t in seconds and the velocity is in mph. " + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [ + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "Gnuplot\n", + "Produced by GNUPLOT 5.0 patchlevel 3 \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t20\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t40\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t60\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t80\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t100\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t120\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t140\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t20\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t40\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t60\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t80\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t100\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t120\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\n", + "\t\n", + "\t\tv (mph)\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\tt (s)\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\tdata\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\tdata\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\n", + "\t\n", + "\tremoved data point\n", + "\n", + "\t\n", + "\t\tremoved data point\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\n", + "\t\n", + "\tlinear\n", + "\n", + "\t\n", + "\t\tlinear\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tcubic spline\n", + "\n", + "\t\n", + "\t\tcubic spline\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tpiecewise cubic\n", + "\n", + "\t\n", + "\t\tpiecewise cubic\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "t=[0 20 40 68 80 84 96 104 110]';\n", + "v=[0 20 20 80 80 100 100 125 125]';\n", + "tt=linspace(0,110)';\n", + "v_lin=interp1(t,v,tt);\n", + "v_spl=interp1(t,v,tt,'spline');\n", + "v_cub=interp1(t,v,tt,'cubic');\n", + "\n", + "plot(t,v,'o',56,38,'s',tt,v_lin,tt,v_spl,tt,v_cub)\n", + "xlabel('t (s)')\n", + "ylabel('v (mph)')\n", + "legend('data','removed data point','linear','cubic spline','piecewise cubic','Location','NorthWest')" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [ + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "Gnuplot\n", + "Produced by GNUPLOT 5.0 patchlevel 3 \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\t\n", + "\t\t-4\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-2\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t2\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t4\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t6\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t8\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t20\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t40\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t60\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t80\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t100\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t120\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\n", + "\t\n", + "\t\tdv/dt (mph/s)\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\tt (s)\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\tlinear\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\tlinear\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tcubic spline\n", + "\n", + "\t\n", + "\t\tcubic spline\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tpiecewise cubic\n", + "\n", + "\t\n", + "\t\tpiecewise cubic\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "t=[0 20 40 56 68 80 84 96 104 110]';\n", + "v=[0 20 20 38 80 80 100 100 125 125]';\n", + "tt=linspace(0,110)';\n", + "v_lin=interp1(t,v,tt);\n", + "v_spl=interp1(t,v,tt,'spline');\n", + "v_cub=interp1(t,v,tt,'cubic');\n", + "\n", + "\n", + "plot(tt(2:end),diff(v_lin)./diff(tt),tt(2:end),diff(v_spl)./diff(tt),tt(2:end),diff(v_cub)./diff(tt))\n", + "xlabel('t (s)')\n", + "ylabel('dv/dt (mph/s)')\n", + "legend('linear','cubic spline','piecewise cubic','Location','NorthWest')" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "### Choose spline wisely\n", + "\n", + "For example of sin(x), not very important\n", + "\n", + "For stop-and-hold examples, the $C^{2}$-continuity should not be preserved. You don't need smooth curves.\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "collapsed": true, + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "# Thanks" + ] + } + ], + "metadata": { + "celltoolbar": "Slideshow", + "kernelspec": { + "display_name": "Octave", + "language": "octave", + "name": "octave" + }, + "language_info": { + "file_extension": ".m", + "help_links": [ + { + "text": "MetaKernel Magics", + "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md" + } + ], + "mimetype": "text/x-octave", + "name": "octave", + "version": "0.19.14" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/16_splines/Newtint.m b/16_splines/Newtint.m new file mode 100644 index 0000000..e4c6c83 --- /dev/null +++ b/16_splines/Newtint.m @@ -0,0 +1,34 @@ +function yint = Newtint(x,y,xx) +% Newtint: Newton interpolating polynomial +% yint = Newtint(x,y,xx): Uses an (n - 1)-order Newton +% interpolating polynomial based on n data points (x, y) +% to determine a value of the dependent variable (yint) +% at a given value of the independent variable, xx. +% input: +% x = independent variable +% y = dependent variable +% xx = value of independent variable at which +% interpolation is calculated +% output: +% yint = interpolated value of dependent variable + +% compute the finite divided differences in the form of a +% difference table +n = length(x); +if length(y)~=n, error('x and y must be same length'); end +b = zeros(n,n); +% assign dependent variables to the first column of b. +b(:,1) = y(:); % the (:) ensures that y is a column vector. +for j = 2:n + for i = 1:n-j+1 + b(i,j) = (b(i+1,j-1)-b(i,j-1))/(x(i+j-1)-x(i)); + end +end +%b +% use the finite divided differences to interpolate +xt = 1; +yint = b(1,1); +for j = 1:n-1 + xt = xt*(xx-x(j)); + yint = yint+b(1,j+1)*xt; +end diff --git a/16_splines/newton_cotes.png b/16_splines/newton_cotes.png new file mode 100644 index 0000000..b734d11 Binary files /dev/null and b/16_splines/newton_cotes.png differ diff --git a/16_splines/octave-workspace b/16_splines/octave-workspace new file mode 100644 index 0000000..8c437bb Binary files /dev/null and b/16_splines/octave-workspace differ diff --git a/16_splines/simpson3.m b/16_splines/simpson3.m new file mode 100644 index 0000000..2ae0c81 --- /dev/null +++ b/16_splines/simpson3.m @@ -0,0 +1,23 @@ +function I = simpson3(func,a,b,n,varargin) +% simpson3: composite simpson's 1/3 rule +% I = simpson3(func,a,b,n,pl,p2,...): +% composite trapezoidal rule +% input: +% func = name of function to be integrated +% a, b = integration limits +% n = number of segments (default = 100) +% pl,p2,... = additional parameters used by func +% output: +% I = integral estimate +if nargin<3,error('at least 3 input arguments required'),end +if ~(b>a),error('upper bound must be greater than lower'),end +if nargin<4|isempty(n),n=100;end +x = a; h = (b - a)/n; + +xvals=linspace(a,b,n+1); +fvals=func(xvals,varargin{:}); +s=fvals(1); +s = s + 4*sum(fvals(2:2:end-1)); +s = s + 2*sum(fvals(3:2:end-2)); +s = s + fvals(end); +I = (b - a) * s/(3*n); diff --git a/16_splines/trap.m b/16_splines/trap.m new file mode 100644 index 0000000..85b8685 --- /dev/null +++ b/16_splines/trap.m @@ -0,0 +1,22 @@ +function I = trap(func,a,b,n,varargin) +% trap: composite trapezoidal rule quadrature +% I = trap(func,a,b,n,pl,p2,...): +% composite trapezoidal rule +% input: +% func = name of function to be integrated +% a, b = integration limits +% n = number of segments (default = 100) +% pl,p2,... = additional parameters used by func +% output: +% I = integral estimate +if nargin<3,error('at least 3 input arguments required'),end +if ~(b>a),error('upper bound must be greater than lower'),end +if nargin<4|isempty(n),n=100;end + +x = a; h = (b - a)/n; +xvals=linspace(a,b,n); +fvals=func(xvals,varargin{:}); +s=func(a,varargin{:}); +s = s + 2*sum(fvals(2:n-1)); +s = s + func(b,varargin{:}); +I = (b - a) * s/(2*n); diff --git a/17_integrals_and_derivatives/.ipynb_checkpoints/17_integrals-checkpoint.ipynb b/17_integrals_and_derivatives/.ipynb_checkpoints/17_integrals-checkpoint.ipynb new file mode 100644 index 0000000..03cd66f --- /dev/null +++ b/17_integrals_and_derivatives/.ipynb_checkpoints/17_integrals-checkpoint.ipynb @@ -0,0 +1,608 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": true, + "slideshow": { + "slide_type": "skip" + } + }, + "outputs": [], + "source": [ + "setdefaults" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": true, + "slideshow": { + "slide_type": "skip" + } + }, + "outputs": [], + "source": [ + "%plot --format svg" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "# Numerical Integrals" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, + "source": [ + "# Integrals in practice\n", + "\n", + "### Example: Compare toughness of Stainless steel to Structural steel\n", + "\n", + "![Stress-strain plot of steel](steel_psi.jpg)\n", + "\n", + "### Step 1 - G3Data to get points \n", + "\n", + "Use the plot shown to determine the toughness of stainless steel and the toughness of structural steel.\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "slide" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "toughness of structural steel is 10.2 psi\n", + "toughness of stainless steel is 18.6 psi\n" + ] + } + ], + "source": [ + "fe_c=load('structural_steel_psi.jpg.dat');\n", + "fe_cr =load('stainless_steel_psi.jpg.dat');\n", + "\n", + "fe_c_toughness=trapz(fe_c(:,1),fe_c(:,2));\n", + "fe_cr_toughness=trapz(fe_cr(:,1),fe_cr(:,2));\n", + "\n", + "fprintf('toughness of structural steel is %1.1f psi\\n',fe_c_toughness)\n", + "fprintf('toughness of stainless steel is %1.1f psi',fe_cr_toughness)\n" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "# Gauss Quadrature (for functions)\n", + "\n", + "Evaluating an integral, we assumed a polynomial form for each Newton-Cotes approximation.\n", + "\n", + "If we can evaluate the function at any point, it makes more sense to choose points more wisely rather than just using endpoints\n", + "\n", + "![trapezoidal example](trap_example.png)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "Set up two unknown constants, $c_{0}$ and $x_{0}$ and determine a *wise* place to evaluate f(x) such that \n", + "\n", + "$I=c_{0}f(x_{0})$\n", + "\n", + "and I is exact for polynomial of n=0, 1\n", + "\n", + "$\\int_{a}^{b}1dx=b-a=c_{0}$\n", + "\n", + "$\\int_{a}^{b}xdx=\\frac{b^2-a^2}{2}=c_{0}x_{0}$\n", + "\n", + "so $c_{0}=b-a$ and $x_{0}=\\frac{b+a}{2}$\n", + "\n", + "$I=\\int_{a}^{b}f(x)dx \\approx (b-a)f\\left(\\frac{b+a}{2}\\right)$" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "f1 =\n", + "\n", + "@(x) x + 1\n", + "\n", + "f2 =\n", + "\n", + "@(x) 1 / 2 * x .^ 2 + x + 1\n", + "\n", + "f3 =\n", + "\n", + "@(x) 1 / 6 * x .^ 3 + 1 / 2 * x .^ 2 + x\n", + "\n" + ] + }, + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "Gnuplot\n", + "Produced by GNUPLOT 5.0 patchlevel 3 \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\t\n", + "\t\t-1000\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-500\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t500\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t1000\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t1500\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-20\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-15\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-10\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t15\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t20\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\n", + "\n", + "\tgnuplot_plot_1a\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "f1=@(x) x+1\n", + "f2=@(x) 1/2*x.^2+x+1\n", + "f3=@(x) 1/6*x.^3+1/2*x.^2+x\n", + "plot(linspace(-18,18),f3(linspace(-18,18)))" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "slide" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "integral of f1 from 2 to 3 = 3.500000\n", + "integral of f1 from 2 to 3 ~ 3.500000\n", + "integral of f2 from 2 to 3 = 6.666667\n", + "integral of f2 from 2 to 3 ~ 6.625000\n" + ] + } + ], + "source": [ + "fprintf('integral of f1 from 2 to 3 = %f',f2(3)-f2(2))\n", + "fprintf('integral of f1 from 2 to 3 ~ %f',(3-2)*f1(3/2+2/2))\n", + "\n", + "fprintf('integral of f2 from 2 to 3 = %f',f3(3)-f3(2))\n", + "fprintf('integral of f2 from 2 to 3 ~ %f',(3-2)*f2(3/2+2/2))\n" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "This process is called **Gauss Quadrature**. Usually, the bounds are fixed at -1 and 1 instead of a and b\n", + "\n", + "$I=c_{0}f(x_{0})$\n", + "\n", + "and I is exact for polynomial of n=0, 1\n", + "\n", + "$\\int_{-1}^{1}1dx=b-a=c_{0}$\n", + "\n", + "$\\int_{-1}^{1}xdx=\\frac{1^2-(-1)^2}{2}=c_{0}x_{0}$\n", + "\n", + "so $c_{0}=2$ and $x_{0}=0$\n", + "\n", + "$I=\\int_{-1}^{1}f(x)dx \\approx 2f\\left(0\\right)$" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, + "source": [ + "Now, integrals can be performed with a change of variable\n", + "\n", + "a=2\n", + "\n", + "b=3\n", + "\n", + "x= 2 to 3\n", + "\n", + "or $x_{d}=$ -1 to 1\n", + "\n", + "$x=a_{1}+a_{2}x_{d}$\n", + "\n", + "at $x_{d}=-1$, x=a\n", + "\n", + "at $x_{d}=1$, x=b" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, + "source": [ + "so \n", + "\n", + "$x=\\frac{(b+a) +(b-a)x_{d}}{2}$\n", + "\n", + "$dx=\\frac{b-a}{2}dx_{d}$\n", + "\n", + "$\\int_{2}^{3}x+1dx=\\int_{-1}^{1}\\left(\\frac{(2+3) +(3-2)x_{d}}{2}\n", + "+1\\right)\n", + "\\frac{3-2}{2}dx_{d}$\n", + "\n", + "$\\int_{2}^{3}x+1dx=\\int_{-1}^{1}\\left(\\frac{5 +x_{d}}{2}\n", + "+1\\right)\n", + "\\frac{3-2}{2}dx_{d}$" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "source": [ + "$\\int_{2}^{3}x+1dx=\\int_{-1}^{1}\\left(\\frac{7}{4}+\\frac{1}{4}x_{d}\\right)dx_{d}=2\\frac{7}{4}=3.5$" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": { + "collapsed": true, + "slideshow": { + "slide_type": "slide" + } + }, + "outputs": [], + "source": [ + "function I=gauss_1pt(func,a,b)\n", + " % Gauss quadrature using single point\n", + " % exact for n<1 polynomials\n", + " c0=2;\n", + " xd=0;\n", + " dx=(b-a)/2;\n", + " x=(b+a)/2+(b-a)/2*xd;\n", + " I=func(x).*dx*c0;\n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "fragment" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "ans = 3.5000\r\n" + ] + } + ], + "source": [ + "gauss_1pt(f1,2,3)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## General Gauss weights and points\n", + "\n", + "![Gauss quadrature table](gauss_weights.png)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "### If you need to evaluate an integral, to increase accuracy, increase number of Gauss points\n", + "\n", + "### Adaptive Quadrature\n", + "\n", + "Matlab/Octave built-in functions use two types of adaptive quadrature to increase accuracy of integrals of functions. \n", + "\n", + "1. `quad`: Simpson quadrature good for nonsmooth functions\n", + "\n", + "2. `quadl`: Lobatto quadrature good for smooth functions\n", + "\n", + "```matlab\n", + "q = quad(fun, a, b, tol, trace, p1, p2, …)\n", + "fun : function to be integrates\n", + "a, b: integration bounds\n", + "tol: desired absolute tolerance (default: 10-6)\n", + "trace: flag to display details or not\n", + "p1, p2, …: extra parameters for fun\n", + "quadl has the same arguments\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "ans = 6.6667\n", + "ans = 6.6667\n", + "ans = 6.6667\n" + ] + } + ], + "source": [ + "% integral of quadratic\n", + "quad(f2,2,3)\n", + "quadl(f2,2,3)\n", + "f3(3)-f3(2)" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "f_c =\n", + "\n", + "@(x) cosh (x / 30) + exp (-10 * x) .* erf (x)\n", + "\n", + "ans = 2.0048\n", + "ans = 2.0048\n" + ] + } + ], + "source": [ + "f_c=@(x) cosh(x/30)+exp(-10*x).*erf(x)\n", + "\n", + "quad(f_c,1,3)\n", + "quadl(f_c,1,3)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "# Thanks" + ] + } + ], + "metadata": { + "celltoolbar": "Slideshow", + "kernelspec": { + "display_name": "Octave", + "language": "octave", + "name": "octave" + }, + "language_info": { + "file_extension": ".m", + "help_links": [ + { + "text": "MetaKernel Magics", + "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md" + } + ], + "mimetype": "text/x-octave", + "name": "octave", + "version": "0.19.14" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/17_integrals_and_derivatives/.ipynb_checkpoints/lecture_20-checkpoint.ipynb b/17_integrals_and_derivatives/.ipynb_checkpoints/lecture_20-checkpoint.ipynb new file mode 100644 index 0000000..2fd6442 --- /dev/null +++ b/17_integrals_and_derivatives/.ipynb_checkpoints/lecture_20-checkpoint.ipynb @@ -0,0 +1,6 @@ +{ + "cells": [], + "metadata": {}, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/17_integrals_and_derivatives/17_integrals.ipynb b/17_integrals_and_derivatives/17_integrals.ipynb new file mode 100644 index 0000000..03cd66f --- /dev/null +++ b/17_integrals_and_derivatives/17_integrals.ipynb @@ -0,0 +1,608 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": true, + "slideshow": { + "slide_type": "skip" + } + }, + "outputs": [], + "source": [ + "setdefaults" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": true, + "slideshow": { + "slide_type": "skip" + } + }, + "outputs": [], + "source": [ + "%plot --format svg" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "# Numerical Integrals" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, + "source": [ + "# Integrals in practice\n", + "\n", + "### Example: Compare toughness of Stainless steel to Structural steel\n", + "\n", + "![Stress-strain plot of steel](steel_psi.jpg)\n", + "\n", + "### Step 1 - G3Data to get points \n", + "\n", + "Use the plot shown to determine the toughness of stainless steel and the toughness of structural steel.\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "slide" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "toughness of structural steel is 10.2 psi\n", + "toughness of stainless steel is 18.6 psi\n" + ] + } + ], + "source": [ + "fe_c=load('structural_steel_psi.jpg.dat');\n", + "fe_cr =load('stainless_steel_psi.jpg.dat');\n", + "\n", + "fe_c_toughness=trapz(fe_c(:,1),fe_c(:,2));\n", + "fe_cr_toughness=trapz(fe_cr(:,1),fe_cr(:,2));\n", + "\n", + "fprintf('toughness of structural steel is %1.1f psi\\n',fe_c_toughness)\n", + "fprintf('toughness of stainless steel is %1.1f psi',fe_cr_toughness)\n" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "# Gauss Quadrature (for functions)\n", + "\n", + "Evaluating an integral, we assumed a polynomial form for each Newton-Cotes approximation.\n", + "\n", + "If we can evaluate the function at any point, it makes more sense to choose points more wisely rather than just using endpoints\n", + "\n", + "![trapezoidal example](trap_example.png)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "Set up two unknown constants, $c_{0}$ and $x_{0}$ and determine a *wise* place to evaluate f(x) such that \n", + "\n", + "$I=c_{0}f(x_{0})$\n", + "\n", + "and I is exact for polynomial of n=0, 1\n", + "\n", + "$\\int_{a}^{b}1dx=b-a=c_{0}$\n", + "\n", + "$\\int_{a}^{b}xdx=\\frac{b^2-a^2}{2}=c_{0}x_{0}$\n", + "\n", + "so $c_{0}=b-a$ and $x_{0}=\\frac{b+a}{2}$\n", + "\n", + "$I=\\int_{a}^{b}f(x)dx \\approx (b-a)f\\left(\\frac{b+a}{2}\\right)$" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "f1 =\n", + "\n", + "@(x) x + 1\n", + "\n", + "f2 =\n", + "\n", + "@(x) 1 / 2 * x .^ 2 + x + 1\n", + "\n", + "f3 =\n", + "\n", + "@(x) 1 / 6 * x .^ 3 + 1 / 2 * x .^ 2 + x\n", + "\n" + ] + }, + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "Gnuplot\n", + "Produced by GNUPLOT 5.0 patchlevel 3 \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\t\n", + "\t\t-1000\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-500\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t500\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t1000\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t1500\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-20\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-15\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-10\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t15\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t20\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\n", + "\n", + "\tgnuplot_plot_1a\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "f1=@(x) x+1\n", + "f2=@(x) 1/2*x.^2+x+1\n", + "f3=@(x) 1/6*x.^3+1/2*x.^2+x\n", + "plot(linspace(-18,18),f3(linspace(-18,18)))" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "slide" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "integral of f1 from 2 to 3 = 3.500000\n", + "integral of f1 from 2 to 3 ~ 3.500000\n", + "integral of f2 from 2 to 3 = 6.666667\n", + "integral of f2 from 2 to 3 ~ 6.625000\n" + ] + } + ], + "source": [ + "fprintf('integral of f1 from 2 to 3 = %f',f2(3)-f2(2))\n", + "fprintf('integral of f1 from 2 to 3 ~ %f',(3-2)*f1(3/2+2/2))\n", + "\n", + "fprintf('integral of f2 from 2 to 3 = %f',f3(3)-f3(2))\n", + "fprintf('integral of f2 from 2 to 3 ~ %f',(3-2)*f2(3/2+2/2))\n" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "This process is called **Gauss Quadrature**. Usually, the bounds are fixed at -1 and 1 instead of a and b\n", + "\n", + "$I=c_{0}f(x_{0})$\n", + "\n", + "and I is exact for polynomial of n=0, 1\n", + "\n", + "$\\int_{-1}^{1}1dx=b-a=c_{0}$\n", + "\n", + "$\\int_{-1}^{1}xdx=\\frac{1^2-(-1)^2}{2}=c_{0}x_{0}$\n", + "\n", + "so $c_{0}=2$ and $x_{0}=0$\n", + "\n", + "$I=\\int_{-1}^{1}f(x)dx \\approx 2f\\left(0\\right)$" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, + "source": [ + "Now, integrals can be performed with a change of variable\n", + "\n", + "a=2\n", + "\n", + "b=3\n", + "\n", + "x= 2 to 3\n", + "\n", + "or $x_{d}=$ -1 to 1\n", + "\n", + "$x=a_{1}+a_{2}x_{d}$\n", + "\n", + "at $x_{d}=-1$, x=a\n", + "\n", + "at $x_{d}=1$, x=b" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, + "source": [ + "so \n", + "\n", + "$x=\\frac{(b+a) +(b-a)x_{d}}{2}$\n", + "\n", + "$dx=\\frac{b-a}{2}dx_{d}$\n", + "\n", + "$\\int_{2}^{3}x+1dx=\\int_{-1}^{1}\\left(\\frac{(2+3) +(3-2)x_{d}}{2}\n", + "+1\\right)\n", + "\\frac{3-2}{2}dx_{d}$\n", + "\n", + "$\\int_{2}^{3}x+1dx=\\int_{-1}^{1}\\left(\\frac{5 +x_{d}}{2}\n", + "+1\\right)\n", + "\\frac{3-2}{2}dx_{d}$" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "source": [ + "$\\int_{2}^{3}x+1dx=\\int_{-1}^{1}\\left(\\frac{7}{4}+\\frac{1}{4}x_{d}\\right)dx_{d}=2\\frac{7}{4}=3.5$" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": { + "collapsed": true, + "slideshow": { + "slide_type": "slide" + } + }, + "outputs": [], + "source": [ + "function I=gauss_1pt(func,a,b)\n", + " % Gauss quadrature using single point\n", + " % exact for n<1 polynomials\n", + " c0=2;\n", + " xd=0;\n", + " dx=(b-a)/2;\n", + " x=(b+a)/2+(b-a)/2*xd;\n", + " I=func(x).*dx*c0;\n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "fragment" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "ans = 3.5000\r\n" + ] + } + ], + "source": [ + "gauss_1pt(f1,2,3)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## General Gauss weights and points\n", + "\n", + "![Gauss quadrature table](gauss_weights.png)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "### If you need to evaluate an integral, to increase accuracy, increase number of Gauss points\n", + "\n", + "### Adaptive Quadrature\n", + "\n", + "Matlab/Octave built-in functions use two types of adaptive quadrature to increase accuracy of integrals of functions. \n", + "\n", + "1. `quad`: Simpson quadrature good for nonsmooth functions\n", + "\n", + "2. `quadl`: Lobatto quadrature good for smooth functions\n", + "\n", + "```matlab\n", + "q = quad(fun, a, b, tol, trace, p1, p2, …)\n", + "fun : function to be integrates\n", + "a, b: integration bounds\n", + "tol: desired absolute tolerance (default: 10-6)\n", + "trace: flag to display details or not\n", + "p1, p2, …: extra parameters for fun\n", + "quadl has the same arguments\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "ans = 6.6667\n", + "ans = 6.6667\n", + "ans = 6.6667\n" + ] + } + ], + "source": [ + "% integral of quadratic\n", + "quad(f2,2,3)\n", + "quadl(f2,2,3)\n", + "f3(3)-f3(2)" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": { + "collapsed": false, + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "f_c =\n", + "\n", + "@(x) cosh (x / 30) + exp (-10 * x) .* erf (x)\n", + "\n", + "ans = 2.0048\n", + "ans = 2.0048\n" + ] + } + ], + "source": [ + "f_c=@(x) cosh(x/30)+exp(-10*x).*erf(x)\n", + "\n", + "quad(f_c,1,3)\n", + "quadl(f_c,1,3)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "# Thanks" + ] + } + ], + "metadata": { + "celltoolbar": "Slideshow", + "kernelspec": { + "display_name": "Octave", + "language": "octave", + "name": "octave" + }, + "language_info": { + "file_extension": ".m", + "help_links": [ + { + "text": "MetaKernel Magics", + "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md" + } + ], + "mimetype": "text/x-octave", + "name": "octave", + "version": "0.19.14" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/17_integrals_and_derivatives/gauss_weights.png b/17_integrals_and_derivatives/gauss_weights.png new file mode 100644 index 0000000..9e3f29d Binary files /dev/null and b/17_integrals_and_derivatives/gauss_weights.png differ diff --git a/17_integrals_and_derivatives/octave-workspace b/17_integrals_and_derivatives/octave-workspace new file mode 100644 index 0000000..8c437bb Binary files /dev/null and b/17_integrals_and_derivatives/octave-workspace differ diff --git a/17_integrals_and_derivatives/stainless_steel_psi.jpg.dat b/17_integrals_and_derivatives/stainless_steel_psi.jpg.dat new file mode 100644 index 0000000..dfa2a4a --- /dev/null +++ b/17_integrals_and_derivatives/stainless_steel_psi.jpg.dat @@ -0,0 +1,12 @@ +1.0208494318e-05 1.6556901722 +0.00241601032192 33.1999376148 +0.00420249682757 53.9506164087 +0.00603492155765 82.1777412288 +0.00844582763241 114.552704897 +0.00959768607462 122.017666367 +0.0207793901842 141.840010208 +0.0369377352739 161.610673548 +0.0574942399989 177.181817537 +0.0774314294019 181.959392878 +0.100609815751 174.241771174 +0.117644389936 156.618719826 diff --git a/17_integrals_and_derivatives/steel_psi.jpg b/17_integrals_and_derivatives/steel_psi.jpg new file mode 100644 index 0000000..5781642 Binary files /dev/null and b/17_integrals_and_derivatives/steel_psi.jpg differ diff --git a/17_integrals_and_derivatives/structural_steel_psi.jpg.dat b/17_integrals_and_derivatives/structural_steel_psi.jpg.dat new file mode 100644 index 0000000..10c9ade --- /dev/null +++ b/17_integrals_and_derivatives/structural_steel_psi.jpg.dat @@ -0,0 +1,13 @@ +1.0208494318e-05 1.6556901722 +0.00180179924712 23.2370851913 +0.00242111456908 34.0306538399 +0.00298938741945 36.5170602372 +0.00410551613155 38.1670081313 +0.0113042060414 39.7537909669 +0.026807506079 42.9158720819 +0.0450807109082 46.8799580317 +0.063896667352 49.1768692533 +0.0937667217264 50.5282186886 +0.134122601181 48.4475999405 +0.194912483429 42.0009357786 +0.224198952211 38.3737301413 diff --git a/17_integrals_and_derivatives/trap_example.png b/17_integrals_and_derivatives/trap_example.png new file mode 100644 index 0000000..facc398 Binary files /dev/null and b/17_integrals_and_derivatives/trap_example.png differ