{VERSION 4 0 "IBM INTEL NT" "4.0" } {USTYLETAB {CSTYLE "Maple Input" -1 0 "Courier" 0 1 255 0 0 1 0 1 0 0 1 0 0 0 0 1 }{CSTYLE "2D Math" -1 2 "Times" 0 1 0 0 0 0 0 0 2 0 0 0 0 0 0 1 }{CSTYLE "2D Input" 2 19 "" 0 1 255 0 0 1 0 0 0 0 0 0 0 0 0 1 } {CSTYLE "" -1 256 "" 1 16 108 191 108 0 0 1 0 0 0 0 0 0 0 1 }{CSTYLE " " -1 257 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 258 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 259 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 260 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 } {CSTYLE "" -1 261 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{PSTYLE "Normal " -1 0 1 {CSTYLE "" -1 -1 "Times" 1 12 0 0 0 1 2 2 2 2 2 2 1 1 1 1 }1 1 0 0 0 0 1 0 1 0 2 2 0 1 }} {SECT 0 {EXCHG {PARA 0 "" 0 "" {TEXT 256 13 "Lorentz force" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 836 "We simulate the \+ motion of a charged particle exposed to an electric and magnetic field . The situation corresponds to that encountered in a vacuum tube with \+ a long filament surrounded by a cylindrical anode with a solenoid plac ed on top. This is more complicated than the usual discussion of a cha rged particle with constant velocity (acquired, e.g., upon acceleratio n in an electric field) entering a homogeneous magnetic field (circula r trajectory) or an electric field aligned with the magnetic field (he lical motion as a result of independent influence of the fields on the position variables). In this example the motion is confined to a plan e, the particles are accelerated by an electric field acting in the ra dial direction. Once they acquire a non-zero velocity the Lorentz forc e deflects the particle away from the radial path." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 481 "The Verlet algorithm use d for the Coulomb problem requires a modification to account for the v elocity dependence of the Lorentz force. The equations of motion are w ritten in Cartesian coordinates (E and B are the magnitudes of the ele ctric field acting in the radial direction and the magnetic field acti ng in the z direction respectively). The radial dependence of the elec tric field strength was determined from a solution to the Laplace equa tion in cylindrical coordinates (see " }{TEXT 19 11 "Laplace.mws" } {TEXT -1 23 "), or from Gauss' law (" }{TEXT 19 9 "Gauss.mws" }{TEXT -1 2 ")." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT 258 1 "m" }{TEXT -1 18 " diff(x(t),t$2) = " }{TEXT 260 1 "q" }{TEXT -1 35 " [E*x(t)/(x(t)^2+y(t)^2) + vy(t)*B]" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT 259 1 "m" }{TEXT -1 18 " diff(y(t),t$2) = " } {TEXT 261 1 "q" }{TEXT -1 35 " [E*y(t)/(x(t)^2+y(t)^2) - vx(t)*B]" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 13 "with(linalg):" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 73 "The Cartesian components for the Lorentz \+ force are calculated as follows:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 22 "v:=vector([vx,vy,vz]);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 23 "Bf:=vector([Bx,By,Bz]);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 16 "crossprod(v,Bf);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 77 "Now we restrict the homogeneous magnetic field to be aligned with \+ the z axis:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 18 "subs(Bx=0,By =0,%);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 43 "We choose the mass (uni ts) and charge sign:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "m:= 1; q:=-1;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 131 "We worry about the \+ step-size for small distances r as the force is inversely proportional to r. We should try various small values." }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 20 "dt:=0.05; dt2:=dt^2;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 15 "vx0:=0; vy0:=0;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 181 "We have to start the trajectories at a finite distance (filament \+ thickness or larger). We can also choose finite initial velocities (th ey depend on the temperature of the filament)." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 60 "Let's measure distances i n multiples of the filament radius." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 30 "rx[0]:=1; rx[1]:=rx[0]+dt*vx0;" }}}{EXCHG {PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 30 "ry[0]:=1; ry[1]:=ry[0]+dt*vy0;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 15 "E:='E'; B:='B';" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 178 "The electric field varies with distance accordin g to the following functions. To obtain the electric force one needs t o multiply with the charge of the probe particle (electron)." }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "Ex:=-x/(x^2+y^2);" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "Ey:=-y/(x^2+y^2);" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 72 "Verlet:=proc(n1,E0,B0) local Fx,Fy,vx,vy; global rx,ry,dt,dt2,m,q,Ex,Ey;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 51 "vx:=(rx[n1]-rx[n1-1])/dt; vy:=(ry[n1]-ry[n1-1])/dt;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 51 "Fx:=evalf(q*subs(x=rx[n1],y=ry[n1 ],E0*Ex)+q*B0*vy);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 51 "Fy:=evalf(q*s ubs(x=rx[n1],y=ry[n1],E0*Ey)-q*B0*vx);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 37 "rx[n1+1]:=2*rx[n1]-rx[n1-1]+dt2/m*Fx;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 37 "ry[n1+1]:=2*ry[n1]-ry[n1-1]+dt2/m*Fy;" }}{PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 4 "end:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 7 "N:=200;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 21 "for i from 1 \+ to N do:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 18 "Verlet(i,1,2); od:" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 23 "nops(convert(rx,list));" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 14 "rx[40],ry[40];" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 66 "plot([seq([rx[i],ry[i]],i=0..N)],st yle=point,scaling=constrained);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 74 "This doesn't seem right. Reducing the time step shows a repetitive motion." }}{PARA 0 "" 0 "" {TEXT -1 91 "How can we fix the equations \+ to use the velocities implicitly? Let us write the recursions:" }} {PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 111 "rx:='rx': ry:='ry': dt:='dt': q:='q': i:='i': m:=1: B:='B': f:= 'f': g:='g': rx_ip1:='rx_ip1': ry_ip1:='ry_ip1':" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 65 "eq1:=m*(rx_ip1+rx[i-1]-2*rx[i])/dt^2=f+B*(ry_ ip1-ry[i-1])/(2*dt);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 65 "eq2 :=m*(ry_ip1+ry[i-1]-2*ry[i])/dt^2=g-B*(rx_ip1-rx[i-1])/(2*dt);" } {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 38 "sol:=solve( \{eq1,eq2\},\{rx_ip1,ry_ip1\});" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "assign(sol);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 7 "rx_i p1;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 20 "Now redefine Verlet:" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 80 "Verlet:=proc(n1,E0,B0) local Fx,Fy; global rx,ry,dt,dt2,m,q,Ex,Ey,rx_ip1,ry_ip1;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 43 "Fx:=evalf(q*subs(x=rx[n1],y=ry[n1],E0*Ex));" }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 43 "Fy:=evalf(q*subs(x=rx[n1],y=ry[n1], E0*Ey));" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 38 "rx[n1+1]:=subs(f=Fx,g=F y,B=B0,rx_ip1);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 38 "ry[n1+1]:=subs(f =Fy,g=Fy,B=B0,ry_ip1);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 4 "end:" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 6 "q:=-1;" }}}{EXCHG {PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 20 "dt:=0.05; dt2:=dt^2;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 15 "vx0:=0; vy0:=0;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 30 "rx[0]:=1; rx[1]:=rx[0]+dt*vx0;" }}}{EXCHG {PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 30 "ry[0]:=1; ry[1]:=ry[0]+dt*vy0;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 41 "Now we carry out a large number of steps. " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "N:=1000;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 21 "for i from 1 to N do:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 18 "Verlet(i,1,2); od:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 23 "nops(convert(rx,list));" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "rx[4],ry[4];" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 66 "plot([seq([rx[i],ry[i]],i=0..N)],style=point,scaling= constrained);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 190 "We see that the motion is that of a regularly shaped rosette. To achieve better accur acy (we note that the solution is slipping after a complete round) we \+ would need to reduce the step-size." }}{PARA 0 "" 0 "" {TEXT -1 72 "We can try cases where the solution gets further away from the filament: " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 9 "dt:=0.02;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "N:=5000;" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 21 "for i from 1 to N do:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 21 "Verlet(i,10,1/2); od:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 75 "plot([seq([rx[i*10],ry[i*10]],i=1..N/10)],style=point ,scaling=constrained);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 69 "Is this solution accurate? Let us re-compute with half the step-size:" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 9 "dt:=0.01;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 9 "N:=10000;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 21 "for i from 1 to N do:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 21 "Verlet(i,10,1/2); od:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 75 "plot([seq([rx[i*20],ry[i*20]],i=1..N/20)],style=point,scaling=co nstrained);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 27 "Now we should be c oncerned:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 10 "dt:=0.005;" }} }{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 9 "N:=20000;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 21 "for i from 1 to N do:" }}{PARA 0 "> " 0 " " {MPLTEXT 1 0 21 "Verlet(i,10,1/2); od:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 75 "plot([seq([rx[i*40],ry[i*40]],i=1..N/40)],style=point ,scaling=constrained);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 121 "It app ears as if the motion should be quite regular, i.e., the 'circles' sho uld be progressing steadily in the x-y plane." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 10 "dt:=0.002;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 9 "N:=50000;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 21 "for i from 1 to N do:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 21 "Verlet (i,10,1/2); od:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 78 "plot([se q([rx[i*100],ry[i*100]],i=1..N/100)],style=point,scaling=constrained); " }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 53 "We seem to be getting a diffe rent answer each time..." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 10 "dt:=0.001;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 10 "N:=100000;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 21 "for i from 1 to N do:" }}{PARA 0 "> " 0 " " {MPLTEXT 1 0 21 "Verlet(i,10,1/2); od:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 78 "plot([seq([rx[i*200],ry[i*200]],i=1..N/200)],style=po int,scaling=constrained);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 294 "Now that we have a satisfying answer (which takes a long time to compute) we turn to Maple's dsolve[numeric] engine. It can handle the solution much more quickly, as it does not use a fixed step-size. It slows dow n in the time-step only when the charged particle returns to the radia l distance " }{TEXT 257 1 "r" }{TEXT -1 77 "=1, where we expect intuit ively the solution to be most difficult to achieve." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "restart;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 26 "q:=-1; m:= 1: E0:=1; B0:=2;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 67 "de1:=m* diff(x(t),t$2)=-q*E0*x(t)/(x(t)^2+y(t)^2)+diff(y(t),t)*q*B0;" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 67 "de2:=m*diff(y(t),t$2)=-q*E0* y(t)/(x(t)^2+y(t)^2)-diff(x(t),t)*q*B0;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 65 "#dsolve(\{de1,de2,x(0)=1,D(x)(0)=0,y(0)=1,D(y)(0)=0\} ,\{x(t),y(t)\});" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 104 "Maple tries \+ for a long time to obtain an analytic solution - to no avail. The pair of ODEs is nonlinear." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 98 "s ol:=dsolve(\{de1,de2,x(0)=1,D(x)(0)=0,y(0)=1,D(y)(0)=0\},\{x(t),y(t)\} ,numeric,output=listprocedure):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 37 "X:=subs(sol,x(t)): Y:=subs(sol,y(t)):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 62 "plot(['X(t)','Y(t)',t=0..40],scaling=constrained ,thickness=2);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 91 "It pays off to \+ use a sophisticated differential equation solver. Let us try the other case:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "restart;" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 29 "q:=-1; m:=1: E0:=10; B0:=1/2 ;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 67 "de1:=m*diff(x(t),t$2)= -q*E0*x(t)/(x(t)^2+y(t)^2)+diff(y(t),t)*q*B0;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 67 "de2:=m*diff(y(t),t$2)=-q*E0*y(t)/(x(t)^2+y(t)^2) -diff(x(t),t)*q*B0;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 98 "sol: =dsolve(\{de1,de2,x(0)=1,D(x)(0)=0,y(0)=1,D(y)(0)=0\},\{x(t),y(t)\},nu meric,output=listprocedure):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 37 "X:=subs(sol,x(t)): Y:=subs(sol,y(t)):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 63 "plot(['X(t)','Y(t)',t=0..200],scaling=constrained, thickness=2);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 165 "For very long i ntegration times Maple's built-in numerical method may break down as w ell. One can verify the accuracy by trying different numerical methods , such as " }{TEXT 19 5 "lsode" }{TEXT -1 4 " or " }{TEXT 19 7 "dverk7 8" }{TEXT -1 85 " instead of the usual adaptive Runge-Kutta-Fehlberg 4 /5 method (cf the help pages on " }{TEXT 19 15 "dsolve[numeric]" } {TEXT -1 2 ")." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 113 "sol:=dso lve(\{de1,de2,x(0)=1,D(x)(0)=0,y(0)=1,D(y)(0)=0\},\{x(t),y(t)\},numeri c,method=dverk78,output=listprocedure):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 37 "X:=subs(sol,x(t)): Y:=subs(sol,y(t)):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 63 "plot(['X(t)','Y(t)',t=0..200],scali ng=constrained,thickness=2);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}{MARK "0 0 0" 13 }{VIEWOPTS 1 1 0 1 1 1803 1 1 1 1 } {PAGENUMBERS 0 1 2 33 1 1 }