Simulation of the motion of a particle in a potential field.

Make changes in the parameters as necessary.

A sample for the motion of a test charge in the field of a list of configuration charges is provided.  If using an ordinary potential, such as x^2 or - 1 / r, be sure to give the particle a small initial velocity in a direction perpendicular to the expected motion, and keep all positions and momenta of order unity.

In any case, changes in inc will be necessary depending on the magnitude of the resulting velocities, and changes in ns will affect both the time required to run and the accuracy: larger ns produces more accuracy and takes longer to run.  If inc is changed, ns must be changed proportionately to maintain accuracy.  Also, if inc is changed, f (the number of frames) will have to be changed to simulate the same total elapsed time.

If you are using the sample below, all charges and locations (x coordinates in xc and y coordinates in yc) are to be placed in lists (denoted by {}).  Units are assumed to be Coulombs and Angstroms.  The only other parameters you will have to change for the sample are x0 and y0.

Off[General::spell1];

This notebook uses the Runge-Kutta (RK) algorithm for the numerical integration of "Hamilton's Equations", a set of first order partial differential equations describing the motion of a massive particle in the presence of a potential field.  It then produces an animation of that motion.  The user has control over the potential, initial conditions, number of frames and delta time values (both between frames and between RK iterations).

The notebook requires the animation package.

Needs["Graphics`Animation`"];

The independent variables are defined in the list v.  They are the standard Hamiltonian phase space variables, consisting of the coordinates and momenta in each dimension.  h is the Hamiltonian (sum of kinetic and potential energies).  rk is a Runge-Kutta function which iteratively integrates Hamilton's equations of motion.  Note that we attempt to minimize the number of expressions which need to be evaluated at each step.  The spaces surrounding the replacement operator are significant.  Also, the last expression is not followed by a semicolon so that its value will be passed back to the caller (otherwise NULL is passed back).

r is defined to be the distance from the origin.  It is provided as a convenience for the user who needs to use it in the definition of the potential.

v:={x,y,px,py};

h:=(px^2+py^2)/(2*N[m]) + u;

rk[v0_]:=(a=dt*N[s /. Thread[v->v0]    ];

		b=dt*N[s /. Thread[v->v0+a/2]];
			    			
		c=dt*N[s /. Thread[v->v0+b/2]];
			    			
		d=dt*N[s /. Thread[v->v0+c]  ];
			    			
		v0+(a+2*b+2*c+d)/6
	);
									
r:=Sqrt[x^2+y^2];

qc = {2,-1} * 1.6 * 10^-19;

xc = {-2,0} * 10^-10;

yc = {1,0} * 10^-10;

eps0 = 8.854*10^-12;

V = Sum[qc[[b]]/Sqrt[(x-xc[[b]])^2+(y-yc[[b]])^2], {b,Length[qc]}] / (4 Pi * eps0);

q = -1 * 1.6 * 10^-19;

The following cell contains values which may be changed by the user:

         u is the potential function;
         m is the mass;
         t0 is the start time;
         x0 and y0 are the initial position of the particle;
         px0 and py0 are its initial momenta in the x and y directions, respectively (with explicit 
                             mass factors so that you can think in terms of velocities);
         f is the number of animation frames to produce;
         xmax and ymax specify the size of the graph;
         inc is the simulated time between frames, and
         ns is the number of steps between frames.

u := q * V;

m = 9.31 * 10^-31;

t0 = 0.0;

x0 = 0 * 10^-10;        	y0 = 0.5 * 10^-10;

px0 = 0 * 10^-10 * m;  	py0 = 0 * 10^-10 * m;

f = 20;

xmax = 4 * 10^-10;	          ymax = 4 * 10^-10;

inc = 1 * 10^-17;		ns = 80;

The following cell does all of the work.  s is the system of Hamilton's Equations, in the form of a list of the symbolic values of each partial derivative.  path is a list of graphics objects which will be used to create the animation.  It is appended to to produce a list of the positions of the particle after each set of RK iterations.  z is used to determine the plotsize.  Note the use of the "xa" (,etc.) variables: it is extremely important that the symbolic variables used in the system of equations not have their values changed during the iterations.

After changing the above values as appropriate, evaluate the entire notebook.

Red particle trajectories indicate acceleration, while blue indicate deceleration.

DO NOT CHANGE ANYTHING IN THE CELL BELOW !

s={D[h,px], D[h,py], -D[h,x], -D[h,y]};
path={};
dt=N[inc/ns];
Animate[({xa,ya,pxa,pya}=Nest[rk,N[{x0,y0,px0,py0}],ns];
	AppendTo[path,If[(pxa^2+pya^2)>(px0^2+py0^2),
				RGBColor[1,0,0],RGBColor[0,0,1]]];
	AppendTo[path,Line[{{x0,y0},{xa,ya}}]];
	{ax,ay}={(pxa-px0)/m*inc,(pya-py0)/m*inc};
	{x0,y0,px0,py0}={xa,ya,pxa,pya};
	Graphics[{path},
		AspectRatio->Automatic,
		 PlotRange->{{-xmax,xmax},{-ymax,ymax}},
		 Axes->True, Ticks->None,
		 AxesOrigin->{-xmax,-ymax},
		 Prolog->{Text[FontForm[StringForm["t= ``", t+inc],
		 		     {"Times",9}],Scaled[{0,.925}],{-1,0}],
		 	Text[FontForm[StringForm["x= ``, y= ``", xa, ya],
		 		     {"Times",9}],Scaled[{0,.85}],{-1,0}],
		 	Text[FontForm[StringForm["vx= ``, vy=``", pxa/m, pya/m],
		 		     {"Times",9}],Scaled[{0,.775}],{-1,0}],
		 	Text[FontForm[StringForm["ax= ``, ay= ``", ax, ay],
		 		     {"Times",9}],Scaled[{0,.7}],{-1,0}]}]),
	{t,t0,t0+(f-1)*inc,inc}];