A Matlab notebook

We can mix whatever Word can do (text, graphics, etc.) with live Matlab code and results, including results from graphs.

We have taken the code from page32.m imported it into a Word-file (page32.doc), and cut it into segements.

Code was highlighted and turned into Matlab input cells using the Notebook menu (Define Input Cell submenu). The text containing the code turned green, every command surrounded by [ ] brackets. The execution blocks (groups) were highlighted and joined together using the Group Cells submenu from Notebook.

The Cell Group was executed manually (Evaluate Cell submenu). Matlab results are returned in blue. The graphs come with handles which can be used to re-size them.

% recursion for the Kepler problem with drag force (page32)

k=1

m=1

dt=0.02

t_max=50

F_x=inline('-k*x/(x^2+y^2)^1.5','x','y','k');

F_y=inline('-k*y/(x^2+y^2)^1.5','x','y','k');

% generalize the drag force to allow dependence on powers of v

% this requires two separate components

Fdr_x=inline('-0.01*vx','vx','vy');

Fdr_y=inline('-0.01*vy','vx','vy');

%Fdr_x=inline('-0.01*vx*sqrt(vx^2+vy^2)^3','vx','vy');

%Fdr_y=inline('-0.01*vy*sqrt(vx^2+vy^2)^3','vx','vy');

N=floor(t_max/dt)

xV(1)=1

xV(2)=1

yV(1)=0

yV(2)=1.15*dt

for j=2 : N-1

vx=(xV(j)-xV(j-1))/dt;

vy=(yV(j)-yV(j-1))/dt;

xV(j+1)=2*xV(j)-xV(j-1)+dt^2/m*( F_x(xV(j),yV(j),k) + Fdr_x(vx,vy) );

yV(j+1)=2*yV(j)-yV(j-1)+dt^2/m*( F_y(xV(j),yV(j),k) + Fdr_y(vx,vy) );

end

for j=1 : N

if floor(j/10)*10==j

xP(j/10)=xV(j);

yP(j/10)=yV(j);

end

end

plot(xP,yP,'b+'),title('Stroboscopic trajectory plot')

% the figure switch directs Matlab to keep this figure window

% Note that the plot windows will be stacked on top of each other.

k =

1

m =

1

dt =

0.0200

t_max =

50

N =

2500

xV =

1

xV =

1 1

yV =

0

yV =

0 0.0230

Now this is the interface I like. It allows me to describe my results!

The spiralling orbit shows that… (here goes your text with explanations)

% now define expressions for kinetic and potential energies:

for j=2 : N-1

t(j)=(j-1)*dt;

r=sqrt(xV(j)^2+yV(j)^2);

v_x=(xV(j+1)-xV(j-1))/(2*dt);

v_y=(yV(j+1)-yV(j-1))/(2*dt);

T=0.5*m*(v_x^2+v_y^2);

V=-k/r;

E_m(j)=T+V;

KE(j)=T;

L_z(j)=m*(xV(j)*v_y-yV(j)*v_x);

Ecc(j)=sqrt(1+2*E_m(j)*L_z(j)^2/m/k^2);

end

plot(t,[E_m; KE],'.'),

title('Kinetic Energy (top), Total Energy (bottom) versus Time')

plot(t,L_z,'r.'),title('Angular Momentum versus Time')

% the next plot shows whether the trajectory loses ellipticity:

% a circular orbit has zero eccentricity.

plot(t,Ecc,'r.'),title('Eccentricity versus Time')