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 stepv=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 conditionsfor 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; endplot(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 readablexdot = 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 zarrayfor 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 zarrayfor 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 zarrayfor 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.
8b) Simple pedulum using DAEs and ODE45 (lecture March 13, 2014).
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, 201113) Animation of sunrise and sunset (no dynamics, just animation) done by a student.
14) Double pendulum, from lecture on May 2, 2013 (zipped folder)