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