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.
ODES.
These three examples are covered in the lecture and lecture video
from Tuesday Feb 8, 2011.
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. 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.
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) Sattelite example from lecture on Thurs Feb 17.
Almost the same as the normal modes example above.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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%6) Sattelite example from lecture on Thurs Feb 24. 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..
7) Two-planets example from lecture on Tues March 1. 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.
8) Rotation of drawing and Animation of same. From lecture on Tuesday March 15, 201.
9) Animation of a pendulum. Also, using ODE23. From lecture Thursday March 17, 2011. zipfile
10) Solution to ballistics problem from Prelim 2, question 3.