ex926.m

ch9例题04_龙格库塔电子运动轨迹 matlab 程序

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;