diff --git a/lecture_07/lecture_07.ipynb b/lecture_07/lecture_07.ipynb
index 08d29c4..66fa9a8 100644
--- a/lecture_07/lecture_07.ipynb
+++ b/lecture_07/lecture_07.ipynb
@@ -1435,22 +1435,23 @@
}
],
"metadata": {
+ "anaconda-cloud": {},
"kernelspec": {
- "display_name": "Octave",
- "language": "octave",
- "name": "octave"
+ "display_name": "Python [conda root]",
+ "language": "python",
+ "name": "conda-root-py"
},
"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"
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.5.2"
}
},
"nbformat": 4,
diff --git a/lecture_07/mod_secant.m b/lecture_07/mod_secant.m
index a3caa10..df97983 100644
--- a/lecture_07/mod_secant.m
+++ b/lecture_07/mod_secant.m
@@ -14,13 +14,13 @@
% ea = approximate relative error (%)
% iter = number of iterations
if nargin<3,error('at least 3 input arguments required'),end
-if nargin<4 || isempty(es),es=0.0001;end
+if nargin<4 || isempty(es),es=0.0001;end
if nargin<5 || isempty(maxit),maxit=50;end
-iter = 0;
+iter = 0
while (1)
xrold = xr;
dfunc=(func(xr+dx)-func(xr))./dx;
- xr = xr - func(xr)/dfunc;
+ xr = xr - func(xr)/dfunc; %calculate as an approximation
iter = iter + 1;
if xr ~= 0
ea = abs((xr - xrold)/xr) * 100;
diff --git a/lecture_10/GaussNaive.m b/lecture_10/GaussNaive.m
index d8c60d3..9153c13 100644
--- a/lecture_10/GaussNaive.m
+++ b/lecture_10/GaussNaive.m
@@ -13,7 +13,9 @@
% forward elimination
for k = 1:n-1
for i = k+1:n
- factor = Aug(i,k)/Aug(k,k);
+ factor = Aug(i,k)/Aug(k,k);
+ %this is naive because if a coefficient in the matrix is zero
+ %then divide by zero errot.
Aug(i,k:nb) = Aug(i,k:nb)-factor*Aug(k,k:nb);
end
end
diff --git "a/lecture_15/\340\244\267\001" "b/lecture_15/\340\244\267\001"
deleted file mode 100644
index e69de29..0000000
diff --git "a/lecture_15/\360&\341\001" "b/lecture_15/\360&\341\001"
deleted file mode 100644
index e69de29..0000000
diff --git a/lecture_17/.ipynb_checkpoints/lecture_17-checkpoint.ipynb b/lecture_17/.ipynb_checkpoints/lecture_17-checkpoint.ipynb
index 2fd6442..4d8aeeb 100644
--- a/lecture_17/.ipynb_checkpoints/lecture_17-checkpoint.ipynb
+++ b/lecture_17/.ipynb_checkpoints/lecture_17-checkpoint.ipynb
@@ -1,6 +1,2711 @@
{
- "cells": [],
- "metadata": {},
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": [
+ "setdefaults"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": [
+ "%plot --format svg"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "### Project Ideas so far\n",
+ "\n",
+ "- Nothing yet...probably something heat transfer related\n",
+ "\n",
+ "- Modeling Propulsion or Propagation of Sound Waves\n",
+ "\n",
+ "- Low Thrust Orbital Transfer\n",
+ "\n",
+ "- Tracking motion of a satellite entering orbit until impact\n",
+ "\n",
+ "- What ever you think is best.\n",
+ "\n",
+ "- You had heat transfer project as an option; that sounded cool\n",
+ "\n",
+ "- Heat transfer through a pipe\n",
+ "\n",
+ "- I would prefer to do something with beam/plate mechanics or vibrations than a heat transfer or thermo problem\n",
+ "\n",
+ "- Modeling Couette flow with a pressure gradient using a discretized form of the Navier-Stokes equation and comparing it to the analytical solution\n",
+ "\n",
+ "- Software to instruct a robotic arm to orient itself based on input from a gyroscope"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Questions from you:\n",
+ "\n",
+ "How was your spring break\n",
+ "\n",
+ "Are grades for hw 3 and 4 going to be posted?\n",
+ "\n",
+ "For class occasionally switching to doc cam and working through problems by hand might help get ideas across.\n",
+ "\n",
+ "Are you assigning those who do not have groups to groups?\n",
+ "\n",
+ "How do you graph a standard distribution in Matlab?\n",
+ "\n",
+ "What is the longest code you have written?\n",
+ "\n",
+ "how to approach probability for hw 5?\n",
+ "??\n",
+ "\n",
+ "Will you be assigning groups to people who do not currently have one? \n",
+ "\n",
+ "what are some basic guidelines we should use to brainstorm project ideas?\n",
+ "\n",
+ "Are you a fan of bananas?\n",
+ "\n",
+ "Going through code isn't the most helpful, because you can easily lose interest. But I am not sure what else you can do.\n",
+ "\n",
+ "Has lecture 15 been posted yet? Looking in the repository I can't seem to find it."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Linear Least Squares Regression\n",
+ "\n",
+ "When approximating a set of data as a polynomial, there will always be error introduced (except in 2 cases). \n",
+ "\n",
+ "For a straight line, the actual data points, $y_{i}$ as a function of the independent variable, $x_{i}$, is:\n",
+ "\n",
+ "$y_{i}=a_{0}+a_{1}x_{i}+e_{i}$\n",
+ "\n",
+ "where $a_{0}$ and $a_{1}$ are the intercept and slope of the line and $e_{i}$ is the error between the approximate function and the recorded data point. \n",
+ "\n",
+ "We make the following assumptions in this analysis:\n",
+ "\n",
+ "1. Each x has a fixed value; it is not random and is known without error.\n",
+ "\n",
+ "2. The y values are independent random variables and all have the same variance.\n",
+ "\n",
+ "3. The y values for a given x must be normally distributed.\n",
+ "\n",
+ "The total error is \n",
+ "\n",
+ "$\\sum_{i=1}^{n}e_{i}=\\sum_{i=1}^{n}(y_{i}-a_{0}-a_{1}x_{i})$\n",
+ "\n",
+ "we don't care about the sign though. One approach has been demonstrated to provide a unique solution is minimizing the sum of squares error or\n",
+ "\n",
+ "$S_{r}=\\sum_{i=1}^{n}e_{i}^{2}=\\sum_{i=1}^{n}(y_{i}-a_{0}-a_{1}x_{i})^{2}$\n",
+ "\n",
+ "Where, $S_{r}$ is the sum of squares error (SSE). \n",
+ "\n",
+ "$\\frac{\\partial S_{r}}{\\partial a_{0}}=-2\\sum(y_{i}-a_{0}-a_{1}x_{i})$\n",
+ "\n",
+ "$\\frac{\\partial S_{r}}{\\partial a_{1}}=-2\\sum(y_{i}-a_{0}-a_{1}x_{i})x_{i}$\n",
+ "\n",
+ "The minimum $S_{r}$ occurrs when the partial derivatives are 0. \n",
+ "\n",
+ "$\\sum y_{i}= \\sum a_{0}+\\sum a_{1}x_{i}$\n",
+ "\n",
+ "$\\sum x_{i}y_{i}= \\sum a_{0}x_{i}+\\sum a_{1}x_{i}^{2}$\n",
+ "\n",
+ "$\\left[\\begin{array}{c}\n",
+ "\\sum y_{i}\\\\\n",
+ "\\sum x_{i}y_{i}\\end{array}\\right]=\n",
+ "\\left[\\begin{array}{cc}\n",
+ "n & \\sum x_{i}\\\\\n",
+ "\\sum x_{i} & \\sum x_{i}^{2}\\end{array}\\right]\n",
+ "\\left[\\begin{array}{c}\n",
+ "a_{0}\\\\\n",
+ "a_{1}\\end{array}\\right]$\n",
+ "\n",
+ "\n",
+ "$b=Ax$\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Example \n",
+ "\n",
+ "Find drag coefficient with best-fit line to experimental data\n",
+ "\n",
+ "|i | v (m/s) | F (N) |\n",
+ "|---|---|---|\n",
+ "|1 | 10 | 25 |\n",
+ "|2 | 20 | 70 |\n",
+ "|3 | 30 | 380|\n",
+ "|4 | 40 | 550|\n",
+ "|5 | 50 | 610|\n",
+ "|6 | 60 | 1220|\n",
+ "|7 | 70 | 830 |\n",
+ "|8 |80 | 1450|"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 68,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "a =\n",
+ "\n",
+ " -234.286\n",
+ " 19.470\n",
+ "\n"
+ ]
+ },
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ "
+
+
+
diff --git a/lecture_20/lecture_20.ipynb b/lecture_20/lecture_20.ipynb
index da8c4ac..639f684 100644
--- a/lecture_20/lecture_20.ipynb
+++ b/lecture_20/lecture_20.ipynb
@@ -2169,6 +2169,7 @@
}
],
"metadata": {
+ "anaconda-cloud": {},
"kernelspec": {
"display_name": "Octave",
"language": "octave",
diff --git a/lecture_20/lecture_20.pdf b/lecture_20/lecture_20.pdf
new file mode 100644
index 0000000..3b6d6ca
Binary files /dev/null and b/lecture_20/lecture_20.pdf differ
+
+
+
+
+
+
+
+In [1]:
+
+
+
+
+
+
+setdefaults
+
+
+
+
+
+In [2]:
+
+
+
+
+
+
+%plot --format svg
+
+
+
+
+
+
+
+
+Questions from last class¶
The interp2 function uses splines to interpolate between data points. What are three options for interpolation:
+-
+
cubic spline
+
+piecewise cubic spline
+
+linear spline
+
+quadratic spline
+
+fourth-order spline
+
+
Numerical integration is a general application of the Newton-Cotes formulas. What is the first order approximation of the Newton-Cotes formula? *
+-
+
trapezoidal rule
+
+Simpson's 1/3 rule
+
+Simpson's 3/8 rule
+
+linear approximation of integral
+
+constant approximation of integral (sum(f(x)*dx))
+
+
+
+
+
+
+
+
+
+Questions from you¶
-
+
Is spring here to stay?
+ +
+The time is now.
+
+Final Project Sheet?
+-
+
- coming this evening/tomorrow +
+can you provide some sort of hw answer key?
+-
+
- we can go through some of the HW +
+What's the most 1337 thing you've ever done?
+-
+
- sorry, I'm n00b to this +
+Can we do out more examples by hand (doc cam or drawing on computer notepad) instead of with pre-written code?
+-
+
- forthcoming +
+Favorite movie?
+-
+
- Big Lebowski +
+
+
+
+
+
+
+
+
+
+
+In [49]:
+
+
+
+
+
+
+fe_c=load('structural_steel_psi.jpg.dat');
+fe_cr =load('stainless_steel_psi.jpg.dat');
+
+fe_c_toughness=trapz(fe_c(:,1),fe_c(:,2));
+fe_cr_toughness=trapz(fe_cr(:,1),fe_cr(:,2));
+
+fprintf('toughness of structural steel is %1.1f psi\n',fe_c_toughness)
+fprintf('toughness of stainless steel is %1.1f psi',fe_cr_toughness)
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+Gauss Quadrature (for functions)¶
Evaluating an integral, we assumed a polynomial form for each Newton-Cotes approximation.
+If we can evaluate the function at any point, it makes more sense to choose points more wisely rather than just using endpoints
+
Let us set up two unknown constants, $c_{0}$ and $x_{0}$ and determine a wise place to evaluate f(x) such that
+$I=c_{0}f(x_{0})$
+and I is exact for polynomial of n=0, 1
+$\int_{a}^{b}1dx=b-a=c_{0}$
+$\int_{a}^{b}xdx=\frac{b^2-a^2}{2}=c_{0}x_{0}$
+so $c_{0}=b-a$ and $x_{0}=\frac{b+a}{2}$
+$I=\int_{a}^{b}f(x)dx \approx (b-a)f\left(\frac{b+a}{2}\right)$
+ +
+
+
+
+
+In [147]:
+
+
+
+
+
+
+f1=@(x) x+1
+f2=@(x) 1/2*x.^2+x+1
+f3=@(x) 1/6*x.^3+1/2*x.^2+x
+plot(linspace(-18,18),f3(linspace(-18,18)))
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+In [148]:
+
+
+
+
+
+
+fprintf('integral of f1 from 2 to 3 = %f',f2(3)-f2(2))
+fprintf('integral of f1 from 2 to 3 ~ %f',(3-2)*f1(3/2+2/2))
+
+fprintf('integral of f2 from 2 to 3 = %f',f3(3)-f3(2))
+fprintf('integral of f2 from 2 to 3 ~ %f',(3-2)*f2(3/2+2/2))
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+This process is called Gauss Quadrature. Usually, the bounds are fixed at -1 and 1 instead of a and b
+$I=c_{0}f(x_{0})$
+and I is exact for polynomial of n=0, 1
+$\int_{-1}^{1}1dx=b-a=c_{0}$
+$\int_{-1}^{1}xdx=\frac{1^2-(-1)^2}{2}=c_{0}x_{0}$
+so $c_{0}=2$ and $x_{0}=0$
+$I=\int_{-1}^{1}f(x)dx \approx 2f\left(0\right)$
+Now, integrals can be performed with a change of variable
+a=2
+b=3
+x= 2 to 3
+or $x_{d}=$ -1 to 1
+$x=a_{1}+a_{2}x_{d}$
+at $x_{d}=-1$, x=a
+at $x_{d}=1$, x=b
+so
+$x=\frac{(b+a) +(b-a)x_{d}}{2}$
+$dx=\frac{b-a}{2}dx_{d}$
+$\int_{2}^{3}x+1dx=\int_{-1}^{1}\left(\frac{(2+3) +(3-2)x_{d}}{2} ++1\right) +\frac{3-2}{2}dx_{d}$
+$\int_{2}^{3}x+1dx=\int_{-1}^{1}\left(\frac{5 +x_{d}}{2} ++1\right) +\frac{3-2}{2}dx_{d}$
+ +
+
+
+
+
+
+
+
+$\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$
+ +
+
+
+
+
+In [15]:
+
+
+
+
+
+
+function I=gauss_1pt(func,a,b)
+ % Gauss quadrature using single point
+ % exact for n<1 polynomials
+ c0=2;
+ xd=0;
+ dx=(b-a)/2;
+ x=(b+a)/2+(b-a)/2*xd;
+ I=func(x).*dx*c0;
+end
+
+
+
+
+
+In [13]:
+
+
+
+
+
+
+gauss_1pt(f1,2,3)
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+General Gauss weights and points¶
If you need to evaluate an integral, to increase accuracy, increase number of Gauss points¶
Adaptive Quadrature¶
Matlab/Octave built-in functions use two types of adaptive quadrature to increase accuracy of integrals of functions.
+-
+
+quad: Simpson quadrature good for nonsmooth functions
+
+quadl: Lobatto quadrature good for smooth functions
+
q = quad(fun, a, b, tol, trace, p1, p2, …)
+fun : function to be integrates
+a, b: integration bounds
+tol: desired absolute tolerance (default: 10-6)
+trace: flag to display details or not
+p1, p2, …: extra parameters for fun
+quadl has the same arguments
+
+
+
+
+
+In [149]:
+
+
+
+
+
+
+% integral of quadratic
+quad(f2,2,3)
+quadl(f2,2,3)
+f3(3)-f3(2)
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+Numerical Differentiation¶
Expanding the Taylor Series:
+$f(x_{i+1})=f(x_{i})+f'(x_{i})h+\frac{f''(x_{i})}{2!}h^2+\cdots$
+ +
+
+
+
+
+In [82]:
+
+
+
+
+
+
+x=linspace(-pi,pi);
+y_smooth=sin(x);
+y_noise =y_smooth+rand(size(x))*0.1-0.05;
+plot(x,y_smooth,x,y_noise)
+title('Low noise in sin wave')
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+In [98]:
+
+
+
+
+
+
+% Central Difference derivative
+
+dy_smooth=zeros(size(x));
+dy_smooth([1,end])=NaN;
+dy_smooth(2:end-1)=(y_smooth(3:end)-y_smooth(1:end-2))/2/(x(2)-x(1));
+
+dy_noise=zeros(size(x));
+dy_noise([1,end])=NaN;
+dy_noise(2:end-1)=(y_noise(3:end)-y_noise(1:end-2))/2/(x(2)-x(1));
+
+plot(x,dy_smooth,x,dy_noise)
+title('Noise Amplified with derivative')
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+Reduce noise¶
Options:
+-
+
Fit a function and take derivative
+a. splines won't help much
+b. best fit curve (better)
+
+Smooth data (does not matter if you smooth before/after derivative)
+
+
+
+
+
+
+In [99]:
+
+
+
+
+
+
+y_spline_i1=interp1(x,y_noise,x+0.1);
+y_spline_in1=interp1(x,y_noise,x-0.1);
+dy_spline=(y_spline_i1-y_spline_in1)/0.2;
+plot(x,dy_spline,x,dy_noise)
+legend('deriv of spline','no change')
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+In [100]:
+
+
+
+
+
+
+Z=[sin(x')];
+a=Z\y_noise'
+plot(x,a*sin(x),x,y_noise)
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+In [103]:
+
+
+
+
+
+
+plot(x,a*cos(x),x,dy_smooth,'o')
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+In [120]:
+
+
+
+
+
+
+a=[1,2,3,4,5]
+pa=[a(4/2:-1:1) a a(end:-1:end-4/2+1)]
+size(pa)
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+In [137]:
+
+
+
+
+
+
+% Smooth data
+N=10; % average data between 10 points (forward/backward)
+y_data=[y_noise(N/2:-1:1) y_noise y_noise(end:-1:end-N/2+1)];
+y_filter=y_data;
+for i=6:length(x)
+ y_filter(i)=mean(y_data(i-5:i+5));
+end
+y_filter=y_filter(N/2:end-N/2-1);
+size(y_filter)
+plot(x,y_filter,x,y_noise,'.')
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+In [138]:
+
+
+
+
+
+
+dy_filter=zeros(size(x));
+dy_filter([1,end])=NaN;
+dy_filter(2:end-1)=(y_filter(3:end)-y_filter(1:end-2))/2/(x(2)-x(1));
+
+plot(x,dy_smooth,x,dy_filter,x,dy_noise,'.')
+title('Noise Amplified with derivative')
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+In [ ]:
+
+
+
+
+
+
+
+