update 4/18/17

lecture_23/.ipynb_checkpoints/lecture_23-checkpoint.ipynb

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

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

lecture_23/north_coriolis.m

@@ -0,0 +1,53 @@
+function [t,r] = coriolis(L,m)
+  % 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,F,m,R)
+    Fx=F(1); Fy=F(2); Fz=F(3); % set each Force component
+    g=9.81; % acceleration due to gravity m/s^2
+    we=2*pi/23.934/3600; % rotation of Earth (each day is 23.934 hours long)
+    ddr=zeros(6,1); % initialize ddr
+
+    ddr(1) = r(2); % x North(+) South (-)
+    ddr(2) = 2*we*r(2).*sin(r(1)/R)-r(6).*cos(r(1)/R)+Fx/m; % dx/dt 
+    ddr(3) = r(4); % y West (+) East (-)
+    ddr(4) = -2*we*(r(2).*sin(r(1)/R)-r(6).*cos(r(1)/R))+Fy/m; % dy/dt
+    ddr(5) = r(6); % z altitude
+    ddr(6) = -2*we*r(4).*cos(r(1)/R)+Fz/m-g; % dz/dt
+  end
+
+  R=6378.1e3; % radius of Earth in m
+  % Applied force is 100 N in Northern direction and z-direction equals force of gravity
+  % (neutrally buoyant travel)
+  F= [ 100 0 9.81*m]; %Applied force in N
+  [t,r]=ode45(@(t,r) myode(t,r,[100 0 9.81*m],m, R),[0 200], [L*pi/180*R 0 0 0 10 0]);
+  figure()
+  % Plot Coriolis effect for deviation Westward
+  plot(r(:,3),1e-3*(r(:,1)-r(1,1)))
+  xlabel('West (+m)','Fontsize',18)
+  ylabel('North (+km)','Fontsize',18)
+  text(0,0, 'Storrs, CT')
+  title('Coriolis acceleration with north force')
+
+  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/predprey.m

@@ -0,0 +1,6 @@
+function yp=predprey(t,y,a,b,c,d)
+	% predator-prey model (Lotka-Volterra equations)
+	yp=zeros(2,1);
+	yp(1)=a*y(1)-b*y(1)*y(2);
+	yp(2)=-c*y(2)+d*y(1)*y(2);
+end

lecture_23/rk4sys.m

@@ -0,0 +1,53 @@
+function [tp,yp] = rk4sys(dydt,tspan,y0,h,varargin)
+% rk4sys: fourth-order Runge-Kutta for a system of ODEs
+%   [t,y] = rk4sys(dydt,tspan,y0,h,p1,p2,...): integrates
+%           a system of ODEs with fourth-order RK method
+% input:
+%   dydt = name of the M-file that evaluates the ODEs
+%   tspan = [ti, tf]; initial and final times with output
+%           generated at interval of h, or
+%   = [t0 t1 ... tf]; specific times where solution output
+%   y0 = initial values of dependent variables
+%   h = step size
+%   p1,p2,... = additional parameters used by dydt
+% output:
+%   tp = vector of independent variable
+%   yp = vector of solution for dependent variables
+if nargin<4,error('at least 4 input arguments required'), end
+if any(diff(tspan)<=0),error('tspan not ascending order'), end
+n = length(tspan);
+ti = tspan(1);tf = tspan(n);
+if n == 2
+  t = (ti:h:tf)'; n = length(t);
+  if t(n)<tf
+    t(n+1) = tf;
+    n = n+1;
+  end
+else
+  t = tspan;
+end
+tt = ti; y(1,:) = y0;
+np = 1; tp(np) = tt; yp(np,:) = y(1,:);
+i=1;
+while(1)
+  tend = t(np+1);
+  hh = t(np+1) - t(np);
+  if hh>h,hh = h;end
+  while(1)
+    if tt+hh>tend,hh = tend-tt;end
+    k1 = dydt(tt,y(i,:),varargin{:})';
+    ymid = y(i,:) + k1.*hh./2;
+    k2 = dydt(tt+hh/2,ymid,varargin{:})';
+    ymid = y(i,:) + k2*hh/2;
+    k3 = dydt(tt+hh/2,ymid,varargin{:})';
+    yend = y(i,:) + k3*hh;
+    k4 = dydt(tt+hh,yend,varargin{:})';
+    phi = (k1+2*(k2+k3)+k4)/6;
+    y(i+1,:) = y(i,:) + phi*hh;
+    tt = tt+hh;
+    i=i+1;
+    if tt>=tend,break,end
+  end
+  np = np+1; tp(np) = tt; yp(np,:) = y(i,:);
+  if tt>=tf,break,end
+end