Competition 1

README.md

@@ -1,7 +1,6 @@
 # Computational Mechanics
 ## ME 3255 Spring 2017
 ### Github page: [https://github.uconn.edu/rcc02007/ME3255S2017.git]
-### Public (ipynb rendering)[https://github.com/cooperrc/ME3255S2017]
 
 ### Course Description
 This course introduces students to scientific programming utilizing Matlab/Octave.

final_project/README.md

@@ -56,7 +56,7 @@ 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
+2. 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
 
@@ -71,7 +71,7 @@ differentiation to determine the valuea of the end points
     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
+3. 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
@@ -83,7 +83,7 @@ $\omega=0$
     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
+4. 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
@@ -96,7 +96,7 @@ supported beam if P=0 and q=0.
 
     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
+5. (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. 
 

lecture_22/01_phugoid/01_01_Phugoid_Theory.ipynb

@@ -33,7 +33,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 1,
+   "execution_count": 2,
    "metadata": {
     "collapsed": false
    },
@@ -53,10 +53,10 @@
        "        "
       ],
       "text/plain": [
-       "<IPython.lib.display.YouTubeVideo at 0x7f514c38ffd0>"
+       "<IPython.lib.display.YouTubeVideo at 0x7fef9c611a50>"
       ]
      },
-     "execution_count": 1,
+     "execution_count": 2,
      "metadata": {},
      "output_type": "execute_result"
     }

lecture_23/.ipynb_checkpoints/lecture_23-checkpoint.ipynb

@@ -0,0 +1,6 @@
+{
+ "cells": [],
+ "metadata": {},
+ "nbformat": 4,
+ "nbformat_minor": 2
+}

lecture_23/boussinesq_lookup.m

@@ -0,0 +1,32 @@
+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;
+  if n < 1.3
+      f=fmn(:,2);
+    elseif n > 1.5
+        f=fmn(:,4);
+    else
+        f=fmn(:,3);
+  end
+  [~,i_fit]=sort(abs(m-fmn(:,1)));
+  x=fmn(i_fit(1:4),1);
+  y=f(i_fit(1:4));
+  f_out = Newtint(x,y,m);
+  sigma_z=q*f_out;
+end

lecture_23/coriolis.m

@@ -0,0 +1,63 @@
+function [t,r] = coriolis(L)
+  % In class we ran this function using L=41.8084 (the latitude of Storrs, CT). This
+  % function takes the latitude (L) in degrees and mass (m) and calculates the trajectory
+  % of a particle with a 100 N load directed North. The initial conditions are set as L
+  % (in radians) * radius of Earth, 0 m/s initial x-velocity, 0 m E-W position (add to
+  % -72.261319 degrees for longitude of Storrs), 0 m/s initial E-W velocity, 10 m initial
+  % altitude, 0 m/s initial z-velocity [L*pi/180*R 0 0 0 10 0], the force is given as 100
+  % N North, 0 West and 9.81*m z (neutrally buoyant) [100 0 9.81*m]
+  %
+  % the output of myode is ddr=[dx/dt d2x/dt2 dy/dt d2y/dt2 dz/dt d2z/dt2]' and the input
+  % to myode is r=[x dx/dt y dy/dt z dz/dt]'
+  % using ode23 solver solves for r as a function of time, here we solve from 0 to 200 s
+  % r(:,1) = x (the north-south position from 0 to 200 s)
+  % r(:,3) = x (the West-East position from 0 to 200 s)
+  % r(:,5) = x (the altitude from 0 to 200 s)
+
+  % define ordinary differential equation in terms of 6 first order ODE's
+  function ddr = myode(t,r,R,L)
+    g=9.81; % acceleration due to gravity m/s^2
+    l=10; % 10 m long cable
+    we=2*pi/23.934/3600; % rotation of Earth (each day is 23.934 hours long)
+    ddr=zeros(4,1); % initialize ddr
+
+    ddr(1) = r(2); % x North(+) South (-)
+    ddr(2) = 2*we*r(4).*sin(L)-g/l*r(1); % dx/dt 
+    ddr(3) = r(4); % y West (+) East (-)
+    ddr(4) = -2*we*(r(2).*sin(L))-g/l*r(3); % dy/dt
+  end
+
+  R=6378.1e3; % radius of Earth in m
+  L=L*pi/180;
+  [t,r]=ode45(@(t,r) myode(t,r,R,L),[0 30000], [1 0 0 0 ]);
+  figure()
+  z=-sqrt(10^2-r(:,1).^2-r(:,3).^2);
+  figure(1)
+  we=2*pi/23.934/3600; % rotation of Earth (each day is 23.934 hours long)
+
+  plot(t,tan(we*sin(L)*t),t,-r(:,3)./r(:,1),'.')
+  xlabel('time (s)','Fontsize',18)
+  ylabel('-y/x','Fontsize',18)
+  % Plot Coriolis effect path
+  figure(2)
+  title('Path at 0 hr, 4.1 hrs, 8.3 hrs','Fontsize',24) 
+  N=length(t);
+  i1=[1:100]; i2=[floor(N/2):floor(N/2)+100]; i3=[N-100:N];
+  plot3(r(i1,1),r(i1,3),z(i1))
+  hold on
+  plot3(r(i2,1),r(i2,3),z(i2),'k-')
+  plot3(r(i3,1),r(i3,3),z(i3),'g-')
+  xlabel('X (m)','Fontsize',18)
+  ylabel('Y (m)','Fontsize',18)
+  zlabel('Z (m)','Fontsize',18)
+  title('Coriolis acceleration Foucalt Pendulum')
+
+%  figure()
+%  % Plot Eotvos effect for deviation upwards
+%  plot(1e-3*(r(:,1)-r(1,1)),r(:,5))
+%  xlabel('North (+km)','Fontsize',18)
+%  ylabel('Altitude (+m)','Fontsize',18)
+%  title('Eotvos effect with north force')
+%
+end
+

lecture_23/eulode.m

@@ -0,0 +1,29 @@
+function [t,y] = eulode(dydt,tspan,y0,h,varargin)
+% eulode: Euler ODE solver
+%   [t,y] = eulode(dydt,tspan,y0,h,p1,p2,...):
+%           uses Euler's method to integrate an ODE
+% input:
+%   dydt = name of the M-file that evaluates the ODE
+%   tspan = [ti, tf] where ti and tf = initial and
+%   final values of independent variable
+%   y0 = initial value of dependent variable
+%   h = step size
+%   p1,p2,... = additional parameters used by dydt
+% output:
+%   t = vector of independent variable
+%   y = vector of solution for dependent variable
+if nargin<4,error('at least 4 input arguments required'),end
+ti = tspan(1);tf = tspan(2);
+if ~(tf>ti),error('upper limit must be greater than lower'),end
+t = (ti:h:tf)'; n = length(t);
+% if necessary, add an additional value of t
+% so that range goes from t = ti to tf
+if t(n)<tf
+  t(n+1) = tf;
+  n = n+1;
+end
+y=zeros(n,length(y0));
+y(1,:)=y0;
+for i = 1:n-1 %implement Euler's method
+  y(i+1,:) = y(i,:) + dydt(t(i),y(i,:),varargin{:})'*(t(i+1)-t(i));
+end