Taco Tuesdays

HW6/.~lock.Primary_Energy_monthly.csv#

@@ -0,0 +1 @@
+,ryan,fermi,31.03.2017 16:47,file:///home/ryan/.config/libreoffice/4;

HW6/README.md

@@ -0,0 +1,88 @@
+# 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})=1/m\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)x_{k}^{i}$
+
+   where $x_{k}^{i} is the k-th value of temperature raised to the i-th power (0, and 1)
+
+   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. 
+
+
+

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

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
+

HW6/problem_1_data.m

@@ -0,0 +1,15 @@
+
+% part a
+xa=[1 2 3 4 5]';
+ya=[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]';
+

final_project/README.md

@@ -0,0 +1,103 @@
+# ME 3255 - Final Project
+## Due May 1 by 11:59pm
+
+In this project you are going to solve for the shape of a beam under different loading
+conditions. The shape of the beam varies along the x-axis and as a function of time.  
+
+Notes: Label the plots with legends, x- and y-axis labels and make sure the plots are easy
+to read (you can use the `setdefaults.m` script we have used in class). All functions
+should have a help file and your README.md should describe each file in your repository
+and provide a description of each problem and each solution (use `#`-headings in your file
+to show the start of new problems)
+
+You will be graded both on documentation and implementation of the solutions. 
+
+![Diagram of beam and loading conditions](beam.png)
+
+We will use the Euler-Bernoulli beam equation to describe the shape of the beam, the
+differential equation that governs the solution is:
+
+$\frac{\partial^4 w}{\partial x^4}-\frac{P}{EI}\frac{\partial^2 w}{\partial
+x^2}+\frac{\rho A}{EI}\frac{\partial^2 w}{\partial t^2}=q(x)$     (1)
+
+Where w(x,t) is the displacement of the beam away from the neutral axis as a function of
+position along the beam, x, and time, t, P is the transverse loading of the beam, E is the
+Young's modulus, I is the second moment of Inertia of the beam, $\rho$ is the density, A
+is the cross-sectional area, and q(x) is the transverse distributed load (for a uniform
+pressure, it is the applied pressure times the width of the beam, in units of
+force/length).
+
+We can separate the function $w(x,t)=w(x)e^{i\omega t}$, now equation (1) becomes
+
+$\left(\frac{\partial^4 w}{\partial x^4}-\frac{P}{EI}\frac{\partial^2 w}{\partial
+x^2}-\frac{\rho A \omega^{2}}{EI}w\right)e^{i\omega t}=\frac{q(x)}{EI}$    (2)
+
+For the simply-supported beam shown in Figure 1, the boundary conditions are:
+
+$w(0)=w(L)=0$
+
+$w''(0)=w''(L)=0$
+
+The material is aluminum, E=70 GPa, $\rho$=2700 kg/m$^3$. The bar is 1-m-long with a base
+width, b=0.1 m, and height, h=0.01 m, and the second moment of inertia,
+I=$\frac{bh^3}{12}$.
+
+1. Analytically solve for the shape of the beam if q(x)=cst, P=0, and $\omega$=0 and
+create a function called `shape_simple_support.m` that returns the displacement w(x) given
+q and x
+
+```
+w=shape_simple_support(x,q);
+```
+
+a. Plot q vs the maximum deflection, $\delta x$, of the beam
+
+b. Use a Monte Carlo model to determine the mean and standard deviation  for the
+maximum deflection $\delta x$ if b and h are normally distributed random variables
+with 0.1 % standard deviations at q=50 N/m. 
+
+3. Now use the central difference approximation to set up a system of equations for the
+beam for q(x)=cst, P=0, and $\omega=0$. Use the boundary conditions with a numerical
+differentiation to determine the valuea of the end points
+
+    a. set up the system of equations for 6 segments as a function of q
+
+    b. set up the system of equations for 10 segments as a function of q 
+
+    c. set up the system of equations for 20 segments as a function of q
+
+    d. solve a-c for q=1,10,20,30,50 and plot the numerical results of q vs $\delta x$
+
+    e. Comment on the results from the analytical and numerical approaches (if you used
+    functions then provide help files, if you used scripts, then describe the steps used)
+
+4. Now set up the system of equations using a central difference method if P>0 and
+$\omega=0$
+
+    a. set up the system of equations for 6 segments as a function of q and P
+
+    b. set up the system of equations for 10 segments as a function of q and P
+
+    c. set up the system of equations for 20 segments as a function of q and P
+
+    d. solve a-c for q=1,10,20,30,50 and plot the numerical results of q vs $\delta x$ for
+    P=0, 100, 200, 300 (4 lines, labeled as P=0,P=100,...)
+
+5. Now set up an eigenvalue problem to solve for the natural frequencies of the simply
+supported beam if P=0 and q=0. 
+
+    a. set up the system of equations for 6 segments
+
+    b. set up the system of equations for 10 segments 
+
+    c. set up the system of equations for 20 segments
+
+    d. solve for the natural frequencies ($\omega_{1}$, $\omega_{2}$,...)
+
+    e. Plot the shape of the beam for the first 3 natural frequencies
+
+6. (Bonus 5pt) Create a function to return the system of equations for the eigenvalue
+problem as a function of P, if P>0. Then, plot the lowest natural frequency vs the applied
+load P. 
+
+