Matlab (Updated on Feb 16, 2014)

You want to get good at plotting, basic vector math (example), solving ODEs (below) and eventually animation.
The samples here are not for you to copy, but to learn from. You should be able to write these
all from scratch on your own.

Much of the material on this page, starting at a lower level and ending with optimization (e.g., fmincon), is covered in this excellent looking 49 page pdf by Allan Bower.

ODES.

These three examples (below) are covered in the lecture and lecture video
from Tuesday Feb 8, 2011. Here's how it was done on Jan 29, 2013.
Version 2 with better memory allocation (1/31/13). Version 3 with damping, sinusoidal
forcing and resonance (2/5/13).

1) Here is a simple Euler Method solution for the harmonic oscillator.
It is not efficient. It is just simple. Master this before going on.

%ODE demo with spring-mass
 m=1; k = 1; tf = 2*pi*sqrt(m/k); 
 n = 20000; h = tf/n; %time step
v=0; x = 1; z=[ ]; % initial conditions, z is empty.
for i=1:n;
   a = -(k/m)*x; % F = ma
   x = x + h * v; 
   v = v + h * a;
   z = [z x];
 end

 plot(z)

2) Here's the same idea but putting aside memory ahead of time. (Filling in arrays with zeros.) Master this before going on.

%ODE demo with spring-mass
       m=1; k = 1; tf = 2*pi*sqrt(m/k); 
       n = 10000; h = tf/n; %time step
%Initialize
       x = zeros(n+1,1); v=zeros(n+1,1); t = zeros(n+1,1);
v(1)=0; x(1) = 1; % initial conditions
for i=1:n;
       a = -(k/m)*x(i); % F = ma
       t(i+1) = t(i) + h;
       x(i+1) = x(i) + h * v(i); 
       v(i+1) = v(i) + h * a;
end
plot(t,x); xlabel('t'); ylabel('x');
x(end)

3) Here we try to separate out the ODE from the Euler method. The key is to think of a differential equation that predicts the future.
The rate of change zdot of the state z depends on the present value of the state z, the time t, and any parameters p. The rate of change is
calculated in the "right hand side" called rhs. Master this before going on.

function springmass()
       disp(['Starting at: ' datestr(now)])
 m =1; k = 1; tf = 10*pi*sqrt(m/k);
       p.m = m; p.k = k;
       n = 1000; h = tf/n; %time step
 %Initialize
       tarray = zeros(n+1,1);
       zarray = zeros(n+1,2);
 x0 = 1; v0=0; t(1) = 0; % initial conditions
       z0 = [x0 v0]';
 zarray(1,:) = z0';
       %%%%%Euler's Method%%%%%%%%%%
       for i=1:n;
       z = zarray(i,:)'; t = tarray(i);
       znew = z + h*rhs(t,z,p); %THE KEY LINE!!
       zarray(i+1,:) = znew';
       tarray(i+1) = tarray(i) + h;
       end
       %%%%%%%%%%%%%%%%%%%%%%%%%%%%
 x = zarray(:,1); v=zarray(:,2);
       plot(x,v); xlabel('x'); ylabel('v');
x(end)
end
     
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %RIGHT HAND SIDE
   function zdot = rhs(t,z,p)
       m = p.m; k= p.k; %unpack p to make it more readable
       x = z(1); v = z(2); %unpack z to make it more readable
       xdot = v;
       vdot = -(k/m)*x;
       zdot = [xdot vdot]';
   end
%%%%%%%%%%%%%%%%%%%%%%%%%%%


4) Here we seperate out Euler's method into a seperate function, as
per lecture of Thursday Feb 10, 2011. Know this well (be able to write it on your own) before going on.

 

function springmass4()
 disp(['Starting at: ' datestr(now)])
 p.m =1; p.k = 1; tspan=[0 20];
   n=1000;
   
   x0 = 1; v0=0; z0 = [x0 v0]';
   
   %Call the ODE solver
 [tarray zarray] = eulermethod(tspan,z0,n,p);
 x = zarray(:,1); v=zarray(:,2);
   plot(tarray,x);
  x(end)
 end
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
   %Euler integrator
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%
   function [tarray zarray] = eulermethod(tspan,z0,n,p)
   tarray = linspace(tspan(1),tspan(2),n+1); % initialize tarray
   h = tarray(2)-tarray(1); %size of time step
   zarray = zeros(n+1,2); zarray(1,:) = z0'; % initialize zarray
   for i=1:n;
   %z = zarray(i,:)'; t = tarray(i);
   %znew = z + h*rhs(t,z,p); %THE KEY LINE!!
   %zarray(i+1,:) = znew';
   %The following line is like the 3 above
   zarray(i+1,:) = zarray(i,:) + h*rhs(tarray(i),zarray(i,:)',p)';
   end
   end
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %RIGHT HAND SIDE
   function zdot = rhs(t,z,p)
   m = p.m; k = p.k; % unpack p 
   x = z(1); v = z(2); % unpack z 
   xdot = v; % First ODE
   vdot = -(k/m)*x; % Second ODE
   zdot =[xdot vdot]'; % Pack up z
   end
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%

5) Same as above, but with 2 masses and 3 springs. This was shown in lecture on Thursday Feb 10, 2011. Here's the version from Feb 14, 2013 using the midpoint method. Here is how to derive the equations using symbolic algebra. Here is the version from Feb 13, 2014 using the midpoint method.

 

function normalmodes1()
       disp(['Starting at: ' datestr(now)])
       p.m1 =1; p.m2 =1; 
       p.k1 = 1; p.k2 = 1; p.k3 = 1; 
       tspan=[0 20]; n=10000;
   
       x0 = [1 0]; 
       v0 = [0 0 ]; 
       z0 = [x0 v0]';
   
       %Call the ODE solver
       [tarray zarray] = eulermethod(tspan,z0,n,p);
       x1 = zarray(:,1); x2=zarray(:,2);
       plot(tarray,x1,tarray, x2); 
       disp(['Ending at: ' datestr(now)])
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     
%Euler integrator
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
       function [tarray zarray] = eulermethod(tspan,z0,n,p)
       tarray = linspace(tspan(1),tspan(2),n+1); % initialize tarray
       h = tarray(2)-tarray(1); %size of time step
       zarray = zeros(n+1,length(z0)); zarray(1,:) = z0'; % initialize zarray
   for i=1:n;
       z = zarray(i,:)'; t = tarray(i);
       znew = z + h*rhs(t,z,p); %THE KEY LINE!!
       zarray(i+1,:) = znew';
   end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
       %RIGHT HAND SIDE
       %%%%%%%%%%%%%%%%%%%%%%%%%%%%
       function zdot = rhs(t,z,p)
       m1 = p.m1; m2 = p.m2; k1 = p.k1; k2 = p.k2; k3 = p.k3; 
       x1 = z(1); x2 = z(2); v1 = z(3); v2 = z(4); % unpack z 
       x1dot = v1; x2dot = v2; 
       v1dot = (1/m1)*(-k1*x1 + k2*(x2-x1)); 
       v2dot = (1/m2)*(-k3*x2 + k2*(x1-x2));
       zdot = [x1dot x2dot v1dot v2dot]';
       end
       %%%%%%%%%%%%%%%%%%%%%%%%%%%%

6) Collisions of two masses in 1D. Solve momentum and restitution equations simultaneously using the backslash command.

7) Sattelite example from lecture on Thurs Feb 17, 2011. Almost the same as the normal modes example above.
Here's the same idea, using the mid-point method, from Thurs Feb 21, 2013: The arbitrary particle dynamics problem.

       function sattelite1()
       disp(['Starting at: ' datestr(now)])
       clear all
       p.m1 =1; p.m2 =1; p.G = 1; 
       tspan=[0 7]; n=10000;
   
       r0 = [1 0]; % x and y components
       v0 = [0 1 ]; % x and y components 
       z0 = [r0 v0]';
   
       %Call the ODE solver
       [tarray zarray] = eulermethod(tspan,z0,n,p);
       x = zarray(:,1); y = zarray(:,2);
       vx = zarray(:,3); vy = zarray(:,4);
   
       figure(1); plot(x,y,0,0,'*'); axis('equal')
       figure(2); plot(tarray,x, tarray,y)
       
       shg; % SHow Graphs (bring graphics windows on top)
       disp(['Ending at: ' datestr(now)])
       end
       %%%%%%%%%%%%%%%%%%%%%%%%%%%%
     
       %Euler integrator
       %%%%%%%%%%%%%%%%%%%%%%%%%%%%
       function [tarray zarray] = eulermethod(tspan,z0,n,p)
       
       tarray = linspace(tspan(1),tspan(2),n+1); % initialize tarray
       h = tarray(2)-tarray(1); %size of time step
       zarray = zeros(n+1,length(z0)); zarray(1,:) = z0'; % initialize zarray
         for i=1:n;
           z = zarray(i,:)'; t = tarray(i);
           znew = z + h*rhs(t,z,p); %THE KEY LINE!!
           zarray(i+1,:) = znew';
         end
       
       end
       %%%%%%%%%%%%%%%%%%%%%%%%%%%%

       %RIGHT HAND SIDE
       %%%%%%%%%%%%%%%%%%%%%%%%%%%%
       function zdot = rhs(t,z,p)
       %Assume mass 1 is held in place by a sky hook
       m1 = p.m1; m2 = p.m2; G=p.G;
       r = z(1:2); v = z(3:4); % unpack z, r is position of mass 2
       magr = sqrt(dot(r,r)); % magnitude of r
       ur = r/magr; % unit vector in r direction

       a    = -(G*m1/(magr^2))*ur; % that's vdot (m2 cancels)  %F = ma
   
       rdot = v; % velocity vector of m2, the first two ODEs
       vdot = a; % the second two ODEs
       
       zdot = [rdot;vdot]; % column vector of [xdot ydot vxdot vydot]'
       end
       %%%%%%%%%%%%%%%%%%%%%%%%%%%%
 

8) Sattelite example from lecture on Thurs Feb 24, 2011. Exact same problem as above with the same solution but broken into 3 files. The big new Matlab thing is FEVAL (function evaluation). Here are the 3 files.

Here is the same idea from lecture Thurs Feb 28, 2013, but with two planets. In this one folder are two right-hand-sides (rhs1 and rhs2) corresponding to gravity and spring-dashpot connections. Also, two solvers (Euler and Midpoint). This uses the command FEVAL. You should be able to write and debug all 5 of these functions on your own. Study and play with these first.

9) Two-planets example from lecture on Tues March 1, 2011. Similar to above but with two planets, hence 8 ODEs. The big new Matlab thing is RK2 (midpoint method, much more accurate for given number of steps). Here are the 3 files.

10) Rotation of drawing and Animation of same. From lecture on Tuesday March 15, 2011.
      Similar but from Thursday March 14, 2013 (also showing ODE45 and 'options' using ODESET)

11) Animation of a pendulum. Also, using ODE23. From lecture Thursday March 17, 2011. zipfile

12) Solution to ballistics problem from Prelim 2, question 3, 2011

13) Animation of sunrise and sunset (no dynamics, just animation) done by a student.

14) Double pendulum, from lecture on May 2, 2013 (zipped folder)