www.pudn.com > ch9_ex.rar > ex926.m

function ex926 
clear all,clf,hold off 
global B E; 
e=1.6e-19;        %Charge of electron, coulomb 
m=9.1e-31;        %Electron mass, kg 
B=[0;0;0.1]*e/m;      %Maganetic field strength, tesla 
E=[0;2e4;0]*e/m;      %Electic field strength, vplt/meter 
h=0.5*1e-12;      %Time step, second 
v0=1e5*[-10;2;0.1];   %Initial velocity of electron, meter/sec 
t_final=2*1e-9;     %Time initialization, sec 
 
options=odeset; 
options.relTol=1e-6; 
%[t,v]=ode45(@fun926,[0:h:t_final],v0); 
[t,v]=ode45(@fun926,[0,t_final],v0,options); 
 
figure(1) 
plot(t,v(:,1),t,v(:,2),t,v(:,3)*100) 
xlabel('t'); ylabel('Velocities'); 
 
figure(2) 
plot3(v(:,1),v(:,2),v(:,3)) 
xlabel('Vx'); ylabel('Vy'); zlabel('Vz'); 
 
xyz(1,:)=[0,0,0]; 
for i=2:length(t)-1 
   xyz(i,:)=xyz(i-1,:)+0.5*(v(i-1,:)+v(i,:))*h; 
end 
 
figure(3) 
comet3(xyz(:,1),xyz(:,2),xyz(:,3)) 
 
figure(4) 
plot3(xyz(:,1),xyz(:,2),xyz(:,3)) 
xlabel('x'); ylabel('y'); zlabel('z'); 
view([0,0,1]) 
 
%---------------------------------------------------------- 
function dv=fun926(t,v) 
global B E; 
dv=[v(2)*B(3)-v(3)*B(2); 
      -v(1)*B(3)+v(3)*B(1); 
      v(1)*B(2)-v(2)*B(1)]+E;