diff --git a/HW6/.~lock.Primary_Energy_monthly.csv# b/HW6/.~lock.Primary_Energy_monthly.csv#
new file mode 100644
index 0000000..eef216f
--- /dev/null
+++ b/HW6/.~lock.Primary_Energy_monthly.csv#
@@ -0,0 +1 @@
+,ryan,fermi,31.03.2017 16:47,file:///home/ryan/.config/libreoffice/4;
\ No newline at end of file
diff --git a/HW6/README.md b/HW6/README.md
new file mode 100644
index 0000000..1f3c86c
--- /dev/null
+++ b/HW6/README.md
@@ -0,0 +1,86 @@
+# Homework #6
+## due 4/14 by 11:59pm
+
+
+0. Create a new github repository called 'curve_fitting'.
+
+ a. Add rcc02007 and pez16103 as collaborators.
+
+ b. Clone the repository to your computer.
+
+
+1. Create a least-squares function called `least_squares.m` that accepts a Z-matrix and
+dependent variable y as input and returns the vector of best-fit constants, a, the
+best-fit function evaluated at each point $f(x_{i})$, and the coefficient of
+determination, r2.
+
+```matlab
+[a,fx,r2]=least_squares(Z,y);
+```
+
+ Test your function on the sets of data in script `problem_1_data.m` and show that the
+ following functions are the best fit lines:
+
+ a. y=0.3745+0.98644x+0.84564/x
+
+ b. y=-11.4887+7.143817x-1.04121 x^2+0.046676 x^3
+
+ c. y=4.0046e^(-1.5x)+2.9213e^(-0.3x)+1.5647e^(-0.05x)
+
+
+2. Use the Temperature and failure data from the Challenger O-rings in lecture_18
+(challenger_oring.csv). Your independent variable is temerature and your dependent
+variable is failure (1=fail, 0=pass). Create a function called `cost_logistic.m` that
+takes the vector `a`, and independent variable `x` and dependent variable `y`. Use the
+function, $\sigma(t)=\frac{1}{1+e^{-t}}$ where $t=a_{0}+a_{1}x$. Use the cost function,
+
+ $J(a_{0},a_{1})=\sum_{i=1}^{n}\left[-y_{i}\log(\sigma(t_{i}))-(1-y_{i})\log((1-\sigma(t_{i})))\right]$
+
+ and gradient
+
+ $\frac{\partial J}{\partial a_{i}}=
+ 1/m\sum_{k=1}^{N}\left(\sigma(t_{k})-y_{k}\right)t_{k}$
+
+ a. edit `cost_logistic.m` so that the output is `[J,grad]` or [cost, gradient]
+
+ b. use the following code to solve for a0 and a1
+
+```matlab
+% Set options for fminunc
+options = optimset('GradObj', 'on', 'MaxIter', 400);
+% Run fminunc to obtain the optimal theta
+% This function will return theta and the cost
+[theta, cost] = ...
+fminunc(@(a)(costFunction(a, x, y)), initial_a, options);
+```
+
+ c. plot the data and the best-fit logistic regression model
+
+```matlab
+plot(x,y, x, sigma(a(1)+a(2)*x))
+```
+
+3. The vertical stress under a corner of a rectangular area subjected to a uniform load of
+intensity $q$ is given by the solution of the Boussinesq's equation:
+
+ $\sigma_{z} =
+ \frac{q}{4\pi}\left(\frac{2mn\sqrt{m^{2}+n^{2}+1}}{m^{2}+n^{2}+1+m^{2}n^{2}}\frac{m^{2}+n^{2}+2}{m^{2}+n^{2}+1}+sin^{-1}\left(\frac{2mn\sqrt{m^{2}+n^{2}+1}}{m^{2}+n^{2}+1+m^{2}n^{2}}\right)\right)$
+
+ Typically, this equation is solved as a table of values where:
+
+ $\sigma_{z}=q f(m,n)$
+
+ where $f(m,n)$ is the influence value, q is the uniform load, m=a/z, n=b/z, a and b are
+ width and length of the rectangular area and z is the depth below the area.
+
+ a. Finish the function `boussinesq_lookup.m` so that when you enter a force, q,
+ dimensions of rectangular area a, b, and depth, z, it uses a third-order polynomial
+ interpolation of the four closest values of m to determine the stress in the vertical
+ direction, sigma_z=$\sigma_{z}$. Use a $0^{th}$-order, polynomial interpolation for
+ the value of n (i.e. round to the closest value of n).
+
+ b. Copy the `boussinesq_lookup.m` to a file called `boussinesq_spline.m` and use a
+ cubic spline to interpolate in two dimensions, both m and n, that returns sigma_z.
+
+
+
diff --git a/HW6/README.pdf b/HW6/README.pdf
new file mode 100644
index 0000000..3cc8991
Binary files /dev/null and b/HW6/README.pdf differ
diff --git a/HW6/boussinesq_lookup.m b/HW6/boussinesq_lookup.m
new file mode 100644
index 0000000..004c91c
--- /dev/null
+++ b/HW6/boussinesq_lookup.m
@@ -0,0 +1,22 @@
+function sigma_z=boussinesq_lookup(q,a,b,z)
+ % function that determines stress under corner of an a by b rectangular platform
+ % z-meters below the platform. The calculated solutions are in the fmn data
+ % m=fmn(:,1)
+ % in column 2, fmn(:,2), n=1.2
+ % in column 3, fmn(:,2), n=1.4
+ % in column 4, fmn(:,2), n=1.6
+
+ fmn= [0.1,0.02926,0.03007,0.03058
+ 0.2,0.05733,0.05894,0.05994
+ 0.3,0.08323,0.08561,0.08709
+ 0.4,0.10631,0.10941,0.11135
+ 0.5,0.12626,0.13003,0.13241
+ 0.6,0.14309,0.14749,0.15027
+ 0.7,0.15703,0.16199,0.16515
+ 0.8,0.16843,0.17389,0.17739];
+
+ m=a/z;
+ n=b/z;
+
+ %...
+end
diff --git a/HW6/cost_logistic.m b/HW6/cost_logistic.m
new file mode 100644
index 0000000..169718f
--- /dev/null
+++ b/HW6/cost_logistic.m
@@ -0,0 +1,27 @@
+function [J, grad] = cost_logistic(a, x, y)
+% cost_logistic Compute cost and gradient for logistic regression
+% J = cost_logistic(theta, X, y) computes the cost of using theta as the
+% parameter for logistic regression and the gradient of the cost
+% w.r.t. to the parameters.
+
+% Initialize some useful values
+N = length(y); % number of training examples
+
+% You need to return the following variables correctly
+J = 0;
+grad = zeros(size(theta));
+
+% ====================== YOUR CODE HERE ======================
+% Instructions: Compute the cost of a particular choice of a.
+% Compute the partial derivatives and set grad to the partial
+% derivatives of the cost w.r.t. each parameter in theta
+%
+% Note: grad should have the same dimensions as theta
+%
+
+
+
+% =============================================================
+
+end
+
diff --git a/HW6/problem_1_data.m b/HW6/problem_1_data.m
new file mode 100644
index 0000000..6af2956
--- /dev/null
+++ b/HW6/problem_1_data.m
@@ -0,0 +1,15 @@
+
+% part a
+xa=[1 2 3 4 5]';
+yb=[2.2 2.8 3.6 4.5 5.5]';
+
+% part b
+
+xb=[3 4 5 7 8 9 11 12]';
+yb=[1.6 3.6 4.4 3.4 2.2 2.8 3.8 4.6]';
+
+% part c
+
+xc=[0.5 1 2 3 4 5 6 7 9];
+yc=[6 4.4 3.2 2.7 2.2 1.9 1.7 1.4 1.1];
+
diff --git a/lecture_18/lecture_18.ipynb b/lecture_18/lecture_18.ipynb
index bf7731c..f1fd1b3 100644
--- a/lecture_18/lecture_18.ipynb
+++ b/lecture_18/lecture_18.ipynb
@@ -2,7 +2,7 @@
"cells": [
{
"cell_type": "code",
- "execution_count": 146,
+ "execution_count": 2,
"metadata": {
"collapsed": true
},
@@ -13,7 +13,7 @@
},
{
"cell_type": "code",
- "execution_count": 107,
+ "execution_count": 3,
"metadata": {
"collapsed": true
},
@@ -1959,7 +1959,7 @@
},
{
"cell_type": "code",
- "execution_count": 162,
+ "execution_count": 8,
"metadata": {
"collapsed": false
},
@@ -1969,7 +1969,7 @@
"output_type": "stream",
"text": [
"ln(2)=0.693147\n",
- "ln(2)~0.366349\n"
+ "ln(2)~0.693147\n"
]
},
{
@@ -2014,113 +2014,117 @@
"\n",
"\n",
"\t\n",
- "\t\t\n",
+ "\t\t\n",
"\t\n",
"\n",
"\n",
"\n",
"\n",
- "\t\t\n",
+ "\t\t\n",
"\t\t-2\n",
"\t\n",
"\n",
"\n",
- "\t\t\n",
- "\t\t0\n",
+ "\t\t\n",
+ "\t\t-1.5\n",
"\t\n",
"\n",
"\n",
- "\t\t\n",
- "\t\t2\n",
+ "\t\t\n",
+ "\t\t-1\n",
"\t\n",
"\n",
"\n",
- "\t\t\n",
- "\t\t4\n",
+ "\t\t\n",
+ "\t\t-0.5\n",
"\t\n",
"\n",
"\n",
- "\t\t\n",
- "\t\t6\n",
+ "\t\t\n",
+ "\t\t0\n",
"\t\n",
"\n",
"\n",
- "\t\t\n",
- "\t\t8\n",
+ "\t\t\n",
+ "\t\t0.5\n",
"\t\n",
"\n",
"\n",
- "\t\t\n",
- "\t\t10\n",
+ "\t\t\n",
+ "\t\t1\n",
"\t\n",
"\n",
"\n",
- "\t\t\n",
- "\t\t0\n",
+ "\t\t\n",
+ "\t\t1.5\n",
"\t\n",
"\n",
"\n",
- "\t\t\n",
- "\t\t10\n",
+ "\t\t\n",
+ "\t\t2\n",
"\t\n",
"\n",
"\n",
- "\t\t\n",
- "\t\t20\n",
+ "\t\t\n",
+ "\t\t0\n",
"\t\n",
"\n",
"\n",
- "\t\t\n",
- "\t\t30\n",
+ "\t\t\n",
+ "\t\t1\n",
"\t\n",
"\n",
"\n",
- "\t\t\n",
- "\t\t40\n",
+ "\t\t\n",
+ "\t\t2\n",
"\t\n",
"\n",
"\n",
- "\t\t\n",
- "\t\t50\n",
+ "\t\t\n",
+ "\t\t3\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\t\n",
+ "\t\t4\n",
"\t\n",
"\n",
"\n",
"\t\t\n",
- "\t\t60\n",
+ "\t\t5\n",
"\t\n",
"\n",
"\n",
"\n",
"\n",
- "\t\n",
+ "\t\n",
"\n",
"\n",
"\tgnuplot_plot_1a\n",
"\n",
"\n",
"\n",
- "\t\n",
+ "\t\n",
"\t\n",
"\tgnuplot_plot_2a\n",
"\n",
"\t \n",
- "\t\n",
+ "\t\n",
"\n",
"\t\n",
"\tgnuplot_plot_3a\n",
"\n",
"\t \n",
- "\t\n",
- "\t\n",
+ "\t\n",
"\t\n",
"\tgnuplot_plot_4a\n",
"\n",
- "\t\n",
+ "\t\n",
"\t\n",
"\n",
"\n",
@@ -2144,9 +2148,9 @@
"source": [
"\n",
"\n",
- "x=[0.2,3,10,20,50,60]; % define independent var's\n",
+ "x=[0.2,1,2,3,4]; % define independent var's\n",
"y=log(x); % define dependent var's\n",
- "xx=linspace(min(x),max(x));\n",
+ "xx=linspace(min(x),max(x)+1);\n",
"yy=zeros(size(xx));\n",
"for i=1:length(xx)\n",
" yy(i)=Newtint(x,y,xx(i));\n",
diff --git a/lecture_18/octave-workspace b/lecture_18/octave-workspace
index c651696..d40f9ac 100644
Binary files a/lecture_18/octave-workspace and b/lecture_18/octave-workspace differ
diff --git a/lecture_19/lecture 19.ipynb b/lecture_19/lecture 19.ipynb
index e6f9f76..fde5c6e 100644
--- a/lecture_19/lecture 19.ipynb
+++ b/lecture_19/lecture 19.ipynb
@@ -2,7 +2,7 @@
"cells": [
{
"cell_type": "code",
- "execution_count": 1,
+ "execution_count": 69,
"metadata": {
"collapsed": true
},
@@ -13,7 +13,7 @@
},
{
"cell_type": "code",
- "execution_count": 2,
+ "execution_count": 70,
"metadata": {
"collapsed": true
},
@@ -38,6 +38,283 @@
""
]
},
+ {
+ "cell_type": "code",
+ "execution_count": 74,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "a =\n",
+ "\n",
+ " 1.04210\n",
+ " 8.19609\n",
+ " 0.50283\n",
+ "\n"
+ ]
+ },
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "% for linear regression\n",
+ "x=linspace(0,20)';\n",
+ "y=x.^2 +10*exp(-2*x)+0.5*sinh(x/2)+rand(size(x))*200-100;\n",
+ "Z=[x.^2 exp(-2*x) sinh(x/2)];\n",
+ "a=Z\\y\n",
+ "\n",
+ "plot(x,y,'o',x,a(1)*x.^2+a(2)*exp(-2*x)+a(3)*sinh(x/2))"
+ ]
+ },
{
"cell_type": "markdown",
"metadata": {},
@@ -51,6 +328,8 @@
"source": [
"\n",
"\n",
+ "# Answer: It depends\n",
+ "\n",
"#### Other:\n",
"\n",
"Twice the amount of points needed\n",
@@ -101,7 +380,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
- "# Splines (Brief introduction before next section)\n",
+ "# Splines (Brief introduction before integrals)\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",
@@ -315,7 +594,7 @@
},
{
"cell_type": "code",
- "execution_count": 9,
+ "execution_count": 75,
"metadata": {
"collapsed": false
},
@@ -532,12 +811,12 @@
"source": [
"### Example: Accelerate then hold velocity\n",
"\n",
- "Here the time is given as vector t in seconds and the velocity is in mph. "
+ "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": 10,
+ "execution_count": 82,
"metadata": {
"collapsed": false
},
@@ -683,12 +962,12 @@
"\n",
"\n",
"\n",
- "\t\n",
+ "\t\n",
"\tdata\n",
"\n",
"\n",
"\n",
- "\t\n",
+ "\t\n",
"\t\tdata\n",
"\t\n",
"\n",
@@ -696,7 +975,6 @@
"\t\n",
"\t\n",
- "\tlinear\n",
+ "\tremoved data point\n",
"\n",
- "\t\n",
+ "\t\n",
+ "\t\tremoved data point\n",
+ "\t\n",
+ "\n",
+ "\n",
+ "\t\n",
+ "\t\n",
+ "\tlinear\n",
+ "\n",
+ "\t\n",
"\t\tlinear\n",
"\t\n",
"\n",
"\n",
- "\t\n",
+ "\t\n",
"\t\n",
- "\tcubic spline\n",
+ "\tcubic spline\n",
"\n",
- "\t\n",
+ "\t\n",
"\t\tcubic spline\n",
"\t\n",
"\n",
"\n",
- "\t\n",
+ "\t\n",
"\t\n",
- "\tpiecewise cubic\n",
+ "\tpiecewise cubic\n",
"\n",
- "\t\n",
+ "\t\n",
"\t\tpiecewise cubic\n",
"\t\n",
"\n",
"\n",
- "\t\n",
+ "\t\n",
"\t\n",
- "\n",
+ "\n",
"\n",
"\n",
"\n",
@@ -753,22 +1042,22 @@
}
],
"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",
+ "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',tt,v_lin,tt,v_spl,tt,v_cub)\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','linear','cubic spline','piecewise cubic','Location','NorthWest')"
+ "legend('data','removed data point','linear','cubic spline','piecewise cubic','Location','NorthWest')"
]
},
{
"cell_type": "code",
- "execution_count": 14,
+ "execution_count": 83,
"metadata": {
"collapsed": false
},
@@ -1206,7 +1495,7 @@
},
{
"cell_type": "code",
- "execution_count": 48,
+ "execution_count": 86,
"metadata": {
"collapsed": false
},
@@ -1215,21 +1504,21 @@
"name": "stdout",
"output_type": "stream",
"text": [
- "For 5 steps\n",
- "trapezoid approximation of integral is 0.79 \n",
+ "For 70 steps\n",
+ "trapezoid approximation of integral is 0.99 \n",
" actual integral is 1.00\n"
]
}
],
"source": [
- "N=5;\n",
+ "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": 43,
+ "execution_count": 92,
"metadata": {
"collapsed": false
},
@@ -1365,11 +1654,76 @@
"\t\n",
"\tgnuplot_plot_2a\n",
"\n",
- "\t \n",
+ "\t\t \n",
"\t\n",
- "\t\n",
"\t\n",
@@ -1393,8 +1747,8 @@
}
],
"source": [
- "\n",
- "plot(x,sin(x),linspace(0,pi/2,N),sin(linspace(0,pi/2,N)),'o')"
+ "x=linspace(0,pi);\n",
+ "plot(x,sin(x),linspace(0,pi/2,N),sin(linspace(0,pi/2,N)),'-o')"
]
},
{
@@ -1425,7 +1779,7 @@
},
{
"cell_type": "code",
- "execution_count": 68,
+ "execution_count": 93,
"metadata": {
"collapsed": false
},
diff --git a/lecture_19/octave-workspace b/lecture_19/octave-workspace
new file mode 100644
index 0000000..4bd89dd
Binary files /dev/null and b/lecture_19/octave-workspace differ
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