Skip to content

Commit

Permalink
Merge pull request #12 from rcc02007/master
Browse files Browse the repository at this point in the history
update 4/5/17
mrw13002 committed Apr 5, 2017
2 parents a2926ba + 598519d commit 569d51f
Showing 35 changed files with 8,331 additions and 7 deletions.
1 change: 1 addition & 0 deletions HW6/.~lock.Primary_Energy_monthly.csv#
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
,ryan,fermi,31.03.2017 16:47,file:///home/ryan/.config/libreoffice/4;
86 changes: 86 additions & 0 deletions HW6/README.md
Original file line number Diff line number Diff line change
@@ -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.



Binary file added HW6/README.pdf
Binary file not shown.
22 changes: 22 additions & 0 deletions HW6/boussinesq_lookup.m
Original file line number Diff line number Diff line change
@@ -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
27 changes: 27 additions & 0 deletions HW6/cost_logistic.m
Original file line number Diff line number Diff line change
@@ -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

15 changes: 15 additions & 0 deletions HW6/problem_1_data.m
Original file line number Diff line number Diff line change
@@ -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];

29 changes: 22 additions & 7 deletions lecture_16/lecture_16.ipynb
Original file line number Diff line number Diff line change
@@ -20644,7 +20644,7 @@
},
{
"cell_type": "code",
"execution_count": 123,
"execution_count": 200,
"metadata": {
"collapsed": false
},
@@ -20793,24 +20793,38 @@
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
"</g>\n",
"\t<g id=\"gnuplot_plot_1a\"><title>gnuplot_plot_1a</title>\n",
"<g color=\"white\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"3.00\">\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"1.00\">\n",
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"1.00\">\n",
"\t<path d=\"M80.2,121.7 L80.2,25.7 L246.0,25.7 L246.0,121.7 L80.2,121.7 Z \" stroke=\"black\"/></g>\n",
"\t<g id=\"gnuplot_plot_1a\"><title>data</title>\n",
"<g color=\"white\" fill=\"none\" stroke=\"black\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"3.00\">\n",
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"3.00\">\n",
"\t<g onmousemove=\"gnuplot_svg.showHypertext(evt,'')\" onmouseout=\"gnuplot_svg.hideHypertext()\"><title> </title>\n",
"\t<use color=\"rgb( 0, 0, 255)\" transform=\"translate(71.9,289.8) scale(12.00)\" xlink:href=\"#gpPt5\"/></g>\n",
"\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"end\" transform=\"translate(246.0,55.7)\">\n",
"\t\t<text><tspan font-family=\"{}\">data</tspan></text>\n",
"\t</g>\n",
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"3.00\">\n",
"\t<use color=\"rgb( 0, 0, 255)\" transform=\"translate(71.9,289.8) scale(12.00)\" xlink:href=\"#gpPt5\"/>\n",
"\t<use color=\"rgb( 0, 0, 255)\" transform=\"translate(138.1,283.6) scale(12.00)\" xlink:href=\"#gpPt5\"/>\n",
"\t<use color=\"rgb( 0, 0, 255)\" transform=\"translate(204.2,240.7) scale(12.00)\" xlink:href=\"#gpPt5\"/>\n",
"\t<use color=\"rgb( 0, 0, 255)\" transform=\"translate(270.4,217.2) scale(12.00)\" xlink:href=\"#gpPt5\"/>\n",
"\t<use color=\"rgb( 0, 0, 255)\" transform=\"translate(336.5,208.9) scale(12.00)\" xlink:href=\"#gpPt5\"/>\n",
"\t<use color=\"rgb( 0, 0, 255)\" transform=\"translate(402.7,124.6) scale(12.00)\" xlink:href=\"#gpPt5\"/>\n",
"\t<use color=\"rgb( 0, 0, 255)\" transform=\"translate(468.8,178.5) scale(12.00)\" xlink:href=\"#gpPt5\"/>\n",
"\t<use color=\"rgb( 0, 0, 255)\" transform=\"translate(535.0,92.8) scale(12.00)\" xlink:href=\"#gpPt5\"/>\n",
"\t<use color=\"rgb( 0, 0, 255)\" transform=\"translate(118.3,49.7) scale(12.00)\" xlink:href=\"#gpPt5\"/>\n",
"</g>\n",
"\t</g>\n",
"\t<g id=\"gnuplot_plot_2a\"><title>gnuplot_plot_2a</title>\n",
"\t<g id=\"gnuplot_plot_2a\"><title>best-fit</title>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"3.00\">\n",
"\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"end\" transform=\"translate(246.0,103.7)\">\n",
"\t\t<text><tspan font-family=\"{}\">best-fit</tspan></text>\n",
"\t</g>\n",
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"3.00\">\n",
"\t<path d=\"M71.9,298.7 L138.1,271.8 L204.2,244.9 L270.4,218.0 L336.5,191.0 L402.7,164.1 L468.8,137.2 L535.0,110.3 \" stroke=\"rgb( 0, 128, 0)\"/></g>\n",
"\t<path d=\"M91.4,97.7 L145.2,97.7 M71.9,298.7 L138.1,271.8 L204.2,244.9 L270.4,218.0 L336.5,191.0 L402.7,164.1 L468.8,137.2 L535.0,110.3 \" stroke=\"rgb( 0, 128, 0)\"/></g>\n",
"\t</g>\n",
"<g color=\"white\" fill=\"none\" stroke=\"rgb( 0, 128, 0)\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"2.00\">\n",
"</g>\n",
@@ -20851,6 +20865,7 @@
"a=A\\b\n",
"\n",
"plot(x,y,'o',x,a(1)+a(2)*x)\n",
"legend('data','best-fit','Location','NorthWest')\n",
"xlabel('Force (N)')\n",
"ylabel('velocity (m/s)')"
]
Loading

0 comments on commit 569d51f

Please sign in to comment.