In technical work, often we must deal with equations that do not have an algebraic solution. In this circumstance, we must turn to numerical methods to find approximate values at a set of points. Powerful programs using sophisticated methods are available, but they are often an overkill. In this article I’ll describe some garage mathematics to get quick answers. Usually, it takes only a few minutes if you are familiar with a script language like Perl or Python.
I will use a recent project as an illustration. The figure shows the Pierce geometry for a planar electron injector with gap width D. (The system extends an infinite distance out of the page.) The theory to determine the electrode shapes is discussed in Sect. 7.1 of my book Charged Particle Beams. The equation for the anode surface outside the beam volume is given in terms of the local coordinates (z, x) by
(r/D)^1.3333 cos(4θ/3) = 1,
where r^2 = z^2 + x^2, θ = atan(x/z). I wanted to find values of z at a set of x values. The coordinates could then be used in the Mesh program to define the anode. The surface at large values of x has little effect on fields inside the beam. I picked a practical range 0.0 < x < 2.0.
As in all numerical work, the first task is to reduce the equations to a non-dimensional form so that the calculations have the greatest generality. Here, the task is simple. If we use the variables R = r/D, X = x/D and Z = z/D, then the equations to solve are:
R^1.3333 cos(4θ/3) = 1, 
R^2 = X^2 + Y^2, 
θ = atan(X/Z) = asin(X/R). 
The point at the beam edge has coordinates X = 0.0, Z = 1.0 (R = 1.0, θ = 0.0).
Because we can’t solve the problem directly, we must use an iterative method that walks toward the solution in small steps. There are usually several choices to create an iterative loop, some of which may be unstable. I choose to express the equation in this form
R‘ = 1.0/[cos(4θ/3)]^0.75, 
Z‘ = √(R’^2 – X’^2). 
The iteration proceeds as follows:
- Start from the known point, X0 = 0.0, R0 = 1.0.
- Using the next value X1, approximate the angle as θ1 = asin(X1/R0).
- Solve Eqs. 4 and 5 to find R1.
- Recalculate θ1 = asin(X1/R1) and determine a refined value of R1.
- Repeat the procedure until the change in R1 between cycles is small. Find Z1 from Eq. 2.
Then move to the next point X2 using the value R1 for the initial estimate. Because the entire calculation takes less than a second, it is not necessary to worry about convergence rate or efficiency. Here is a complete Python program:
import math # Scaling d = 1.0 x0 = 0.0 NCycle = 30 XMin = 0.0 XMax = 2.0 NStep = 20 Dx = (XMax-XMin)/NStep RStart = 1.0 # Loop over X values for n in range(NStep+1): x = XMin + n*Dx R = RStart # Iteration loop for m in range(NCycle): theta = math.asin(x/R) Rnew = 1.0/math.pow(math.cos(1.33333*theta),0.75) # Averaging needed for stability R = 0.25*Rnew + 0.75*R z = math.sqrt(R*R-x*x) xdisplay = x0 + d*x zdisplay = d*z print '%6i %9.6f %9.6f' % (n,xdisplay,zdisplay) # Start search at the location of the previous search RStart = R
Most of the time there is some iterative form of the initial equation that gives convergence. This case was more subtle. Initial values were correct, but the calculation became unstable as X approached 1.0. I got valid results over the full range by increasing NStep and adding an averaging statement to slow convergence:
R = 0.25*Rnew + 0.75*R
It was not necessary to open a file for numerical output. I simply copied and pasted the output of the print statement from the Python interpreter window. The output is shown below. With a few simple changes, the program can be modified to write the coordinates in a form that can be integrated directly into a Mesh script.
Results for d = 1.0, x0 = 0.0 0 0.000000 1.000000 1 0.100000 1.001662 2 0.200000 1.006604 3 0.300000 1.014700 4 0.400000 1.025755 5 0.500000 1.039534 6 0.600000 1.055782 ...