diff --git a/lecture_20/.ipynb_checkpoints/lecture_20-checkpoint.ipynb b/lecture_20/.ipynb_checkpoints/lecture_20-checkpoint.ipynb
new file mode 100644
index 0000000..2fd6442
--- /dev/null
+++ b/lecture_20/.ipynb_checkpoints/lecture_20-checkpoint.ipynb
@@ -0,0 +1,6 @@
+{
+ "cells": [],
+ "metadata": {},
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/lecture_20/gauss_weights.png b/lecture_20/gauss_weights.png
new file mode 100644
index 0000000..9e3f29d
Binary files /dev/null and b/lecture_20/gauss_weights.png differ
diff --git a/lecture_20/lecture_20.ipynb b/lecture_20/lecture_20.ipynb
new file mode 100644
index 0000000..50f858a
--- /dev/null
+++ b/lecture_20/lecture_20.ipynb
@@ -0,0 +1,2197 @@
+{
+ "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": [
+ "# Questions from last class\n",
+ "\n",
+ "The interp2 function uses splines to interpolate between data points. What are three options for interpolation:\n",
+ "\n",
+ "- cubic spline\n",
+ "\n",
+ "- piecewise cubic spline\n",
+ "\n",
+ "- linear spline\n",
+ "\n",
+ "- quadratic spline\n",
+ "\n",
+ "- fourth-order spline\n",
+ "\n",
+ "\n",
+ "\n",
+ "Numerical integration is a general application of the Newton-Cotes formulas. What is the first order approximation of the Newton-Cotes formula? *\n",
+ "\n",
+ "- trapezoidal rule\n",
+ "\n",
+ "- Simpson's 1/3 rule\n",
+ "\n",
+ "- Simpson's 3/8 rule\n",
+ "\n",
+ "- linear approximation of integral\n",
+ "\n",
+ "- constant approximation of integral (sum(f(x)*dx))\n",
+ "\n",
+ ""
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Questions from you\n",
+ "\n",
+ "- Is spring here to stay?\n",
+ " \n",
+ " - [Punsxatawney Phil](http://www.groundhog.org/)\n",
+ "\n",
+ "- The time is now.\n",
+ "\n",
+ "- Final Project Sheet?\n",
+ " \n",
+ " - coming this evening/tomorrow\n",
+ "\n",
+ "- can you provide some sort of hw answer key?\n",
+ " \n",
+ " - we can go through some of the HW \n",
+ "\n",
+ "- What's the most 1337 thing you've ever done?\n",
+ "\n",
+ " - sorry, I'm n00b to this\n",
+ "\n",
+ "- Can we do out more examples by hand (doc cam or drawing on computer notepad) instead of with pre-written code?\n",
+ "\n",
+ " - forthcoming\n",
+ "\n",
+ "- Favorite movie?\n",
+ " \n",
+ " - Big Lebowski\n",
+ " "
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Integrals in practice\n",
+ "\n",
+ "### Example: Compare toughness of two steels\n",
+ "\n",
+ "\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": 49,
+ "metadata": {
+ "collapsed": false
+ },
+ "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": {},
+ "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",
+ "\n",
+ "\n",
+ "Let us 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)$\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 147,
+ "metadata": {
+ "collapsed": false
+ },
+ "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": [
+ ""
+ ],
+ "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": 148,
+ "metadata": {
+ "collapsed": false
+ },
+ "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": {},
+ "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)$\n",
+ "\n",
+ "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\n",
+ "\n",
+ "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": {},
+ "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": 15,
+ "metadata": {
+ "collapsed": true
+ },
+ "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": 13,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "ans = 3.5000\r\n"
+ ]
+ }
+ ],
+ "source": [
+ "gauss_1pt(f1,2,3)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## General Gauss weights and points\n",
+ "\n",
+ "\n",
+ "\n",
+ "### 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",
+ "```\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 149,
+ "metadata": {
+ "collapsed": false
+ },
+ "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": "markdown",
+ "metadata": {},
+ "source": [
+ "# Numerical Differentiation\n",
+ "\n",
+ "Expanding the Taylor Series:\n",
+ "\n",
+ "$f(x_{i+1})=f(x_{i})+f'(x_{i})h+\\frac{f''(x_{i})}{2!}h^2+\\cdots$"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 82,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "x=linspace(-pi,pi);\n",
+ "y_smooth=sin(x);\n",
+ "y_noise =y_smooth+rand(size(x))*0.1-0.05;\n",
+ "plot(x,y_smooth,x,y_noise)\n",
+ "title('Low noise in sin wave')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 98,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "% Central Difference derivative\n",
+ "\n",
+ "dy_smooth=zeros(size(x));\n",
+ "dy_smooth([1,end])=NaN;\n",
+ "dy_smooth(2:end-1)=(y_smooth(3:end)-y_smooth(1:end-2))/2/(x(2)-x(1));\n",
+ "\n",
+ "dy_noise=zeros(size(x));\n",
+ "dy_noise([1,end])=NaN;\n",
+ "dy_noise(2:end-1)=(y_noise(3:end)-y_noise(1:end-2))/2/(x(2)-x(1));\n",
+ "\n",
+ "plot(x,dy_smooth,x,dy_noise)\n",
+ "title('Noise Amplified with derivative')"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Reduce noise\n",
+ "\n",
+ "Options:\n",
+ "\n",
+ "1. Fit a function and take derivative\n",
+ " \n",
+ " a. splines won't help much\n",
+ " \n",
+ " b. best fit curve (better)\n",
+ " \n",
+ "2. Smooth data (does not matter if you smooth before/after derivative)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 99,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "y_spline_i1=interp1(x,y_noise,x+0.1);\n",
+ "y_spline_in1=interp1(x,y_noise,x-0.1);\n",
+ "dy_spline=(y_spline_i1-y_spline_in1)/0.2;\n",
+ "plot(x,dy_spline,x,dy_noise)\n",
+ "legend('deriv of spline','no change')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 100,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "a = 0.99838\r\n"
+ ]
+ },
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "Z=[sin(x')];\n",
+ "a=Z\\y_noise'\n",
+ "plot(x,a*sin(x),x,y_noise)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 103,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "plot(x,a*cos(x),x,dy_smooth,'o')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 120,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "a =\n",
+ "\n",
+ " 1 2 3 4 5\n",
+ "\n",
+ "pa =\n",
+ "\n",
+ " 2 1 1 2 3 4 5 5 4\n",
+ "\n",
+ "ans =\n",
+ "\n",
+ " 1 9\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "a=[1,2,3,4,5]\n",
+ "pa=[a(4/2:-1:1) a a(end:-1:end-4/2+1)]\n",
+ "size(pa)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 137,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "ans =\n",
+ "\n",
+ " 1 100\n",
+ "\n"
+ ]
+ },
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "% Smooth data\n",
+ "N=10; % average data between 10 points (forward/backward)\n",
+ "y_data=[y_noise(N/2:-1:1) y_noise y_noise(end:-1:end-N/2+1)];\n",
+ "y_filter=y_data;\n",
+ "for i=6:length(x)\n",
+ " y_filter(i)=mean(y_data(i-5:i+5));\n",
+ "end\n",
+ "y_filter=y_filter(N/2:end-N/2-1);\n",
+ "size(y_filter)\n",
+ "plot(x,y_filter,x,y_noise,'.')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 138,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "dy_filter=zeros(size(x));\n",
+ "dy_filter([1,end])=NaN;\n",
+ "dy_filter(2:end-1)=(y_filter(3:end)-y_filter(1:end-2))/2/(x(2)-x(1));\n",
+ "\n",
+ "plot(x,dy_smooth,x,dy_filter,x,dy_noise,'.')\n",
+ "title('Noise Amplified with derivative')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": []
+ }
+ ],
+ "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/lecture_20/q1.png b/lecture_20/q1.png
new file mode 100644
index 0000000..749476d
Binary files /dev/null and b/lecture_20/q1.png differ
diff --git a/lecture_20/q2.png b/lecture_20/q2.png
new file mode 100644
index 0000000..d61cecd
Binary files /dev/null and b/lecture_20/q2.png differ
diff --git a/lecture_20/stainless_steel_psi.jpg.dat b/lecture_20/stainless_steel_psi.jpg.dat
new file mode 100644
index 0000000..dfa2a4a
--- /dev/null
+++ b/lecture_20/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/lecture_20/steel_psi.jpg b/lecture_20/steel_psi.jpg
new file mode 100644
index 0000000..5781642
Binary files /dev/null and b/lecture_20/steel_psi.jpg differ
diff --git a/lecture_20/structural_steel_psi.jpg.dat b/lecture_20/structural_steel_psi.jpg.dat
new file mode 100644
index 0000000..10c9ade
--- /dev/null
+++ b/lecture_20/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/lecture_20/trap_example.png b/lecture_20/trap_example.png
new file mode 100644
index 0000000..facc398
Binary files /dev/null and b/lecture_20/trap_example.png differ