ex925.m

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

function ex925
clear all,clf,hold off

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.5e-12;      %Time step, second
v(:,1)=1e5*[-10;2;0.1];     %Initial velocity of electron, meter/sec
t(1)=0;     %Time initialization, sec
xyz(:,1)=[0;0;0];

for i=2:4000
    t(i)=h*(i-1);
    k1=vxv_(v(:,i-1),B)+E;
    k2=vxv_(v(:,i-1)+0.5*h*k1,B)+E;
    k3=vxv_(v(:,i-1)+0.5*h*k2,B)+E;
    k4=vxv_(v(:,i-1)+h*k3,B)+E;
    v(:,i)=v(:,i-1)+h*(k1+2*k2+2*k3+k4)/6;
    xyz(:,i)=xyz(:,i-1)+0.5*(v(:,i-1)+v(:,i))*h;
end
v'

figure(1)
comet3(xyz(1,:),xyz(2,:),xyz(3,:))

figure(2)
plot3(xyz(1,:),xyz(2,:),xyz(3,:))
xlabel('x'); ylabel('y'); zlabel('z');
view([0,0,1])

figure(3)
plot3(v(1,:),v(2,:),v(3,:))
xlabel('Vx'); ylabel('Vy'); zlabel('Vz');

figure(4)
plot(t,v(1,:),t,v(2,:),t,v(3,:)*100)
xlabel('t'); ylabel('Velocities');
text(t(90),v(1,90),'Vx')
text(t(1200),v(2,1200),'Vy')
text(t(500),v(3,500)*100,'Vz*100')

function c=vxv_(a,b)
  c=[a(2)*b(3)-a(3)*b(2);
      -a(1)*b(3)+a(3)*b(1);
      a(1)*b(2)-a(2)*b(1)];