### Matlab (Updated on March 1, 2011)

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 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. 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.

```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) 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 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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%```
` `

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.