6  Finite Difference Method

At the dose f, the qi(x,f) functions are deemed to be sufficiently smooth in the region [0,xP(f)]. In this region the governing PDEs can be satisfactorily solved by a finite difference, finite element or a related method. The use of such methods takes the form of using the solution q(x,f*) for x Î [0,xP] and at a particular dose level f* (and perhaps previous dose levels) to compute the solution q(x,f*+d) at the following dose level. Traditionally, for such methods, the term dose f is substituted by time t and the application of such time-stepping methods to parabolic and hyperbolic PDEs are considered in many texts including Richtmyer and Morton [10] and Mitchell and Griffiths [11].
In the original IMPETUS code the most straightforward explicit FDM was employed. Although this method provided a satisfactory numerical solution, the severe restriction it placed on the size of the dose steps was felt to lead to a relatively inefficient method and that a more complex (implicit) method could reduce the execution time of the program considerably.
In this section the solution of the PDE of equations (1)-(3) by the FDM is examined. In section 5 it was demonstrated that there was a transition from parabolic to hyperbolic nature in the PDE. This also affects the development of a computational method. The particular method of solution employed by the IMPETUS II code is described in full.

6.1  The Straightforward Explicit Method

The most straigthforward FDM to apply to the problem to be obtained through writing

v(x,f)
f
 |(x,f)=(x*,f*) » 1
d
 [ v( x*,f*+D)-v(x*,f) ],
(28)

v(x,f)
x
 |(x,f)=(x*,f*) » 1
2 D
 [ v(x*+D,f*)-v(x*-D,f*) ] ,
(29)

2 v(x,f)
x2
 |(x,f)=(x*,f*) » 1
D2
 [ v(x*+D,f*)-2 v(x*,f*) +v(x*-D,f*) ] ,
(30)
The molecule for the above method is given in the figure 7.


Figure 7. The molecule for the explicit method.
The method is explicit; the solution at any particular point on the next dose level can be obtained independent of the other points at that dose level. The method is the easiest to implement and it is the method employed in the original IMPETUS code. However, the method has some unfortunate properties. For the model equation [(v(x,f))/(f)] = U [(v(x,f))/(x)] the Courant - Friedrichs - Lewy condition requires that d £ [(D)/U] and for the model equation [(v(x,f))/(f)] = D [(2 v(x,f))/(x2)] the condition d £ [(D2)/2 D] must be satisfied if the resulting method is to be stable. The application of this method to the solution of the partial differential equation therefore places a severe restriction on the dose step and hence can result in a relatively inefficient numerical method.

6.2  The Linearisation Technique

For the IMPETUS II code the above finite difference method was replaced by a implicit method. In general implicit methods do not place severe resrictions on d, the dose step length. However the methods are more difficult to implement since they require that the entire solution is obtained simultaneously at the f* + d dose level. For non-linear PDEs, such as the ones considered in this paper, the development of the method is further complicated by the inter-dependence of the functions in the equations. For example for the PDEs (1), with the expressions (2), (3) substituted for the qi, the collective current q and the diffusivities Di at the dose level f+ d depend on the (unknown) material distribution of the species qi(x,f*+d) (for i=1,2,..n) at that dose level. The application of an implicit method to such a non-linear problem is only possible either by a solution of a non-linear set of equations (by Newton's method, for example) or by linearising the PDE. Linearisation leads to a linear banded system of equations to be solved at each dose step and the technique is employed to obtain the method used in IMPETUS II.
The functions q and Di are related to the material distribution. Hence in order to obtain a solution q(x,f+d) at the following dose level their values must be fixed over the dose interval in the method. Note that the currents qi are not linearised, their values at a particular depth x may change over the dose interval.

6.3  The Finite Difference Method of IMPETUS II

The FDM employed in IMPETUS II is a two-step, three-level scheme. The advantages of using two-step methods in convection-diffusion equations are considered in Chen [12]. The method employed can be considered to be of predictor-corrector type. The first step advances the solution from the initial dose level f* to f* +d, the functions q and Di being set on at dose f*. The second step of the scheme involves using the solution at the dose levels f* and f*+d to obtain a solution at f*+ 2 d. In the second step the functions q and Di are linearised in respect of their values at the central dose level, hence it is expected that the second step 'corrects' the result of the first step to a higher accuracy.
The first step, advancing the solution from the f* dose level to f*+d, is carried out by the Crank-Nicolson method. In this method the non-linearised terms in the PDE are approximated as follows:

v(x,f)
f
 |(x,f)=(x*,f*) » 1
d
 [ v( x*,f*+d)-v(x*,f) ],
(31)

v(x,f)
x
 |(x,f)=(x*,f*) » 1
4 D
 [ v(x*+D,f*+d)-v(x*-D,f*+d) +v(x*+D,f*)-v(x*-D,f*) ] ,
(32)

2 v(x,f)
x2
 |(x,f)=(x*,f*) » 1
D2
 [ v(x*+D,f*+d)-2 v(x*,f*+d) +v(x*-D,f*+d)

              + v(x*+D,f*)-2 v(x*,f*) +v(x*-D,f*) ] ,
(33)
and there is no stability restriction on the dose step d.
The molecule for the method is given in figure 8.

Figure 8. The molecule for the Crank-Nicolson method.

In the second step the solution at the initial (first) and second dose levels are used to obtain the solution at the third level. In this method the non-linearised terms in the PDE are approximated as follows:
v(x,f)
f
 |(x,f)=(x*,f*) » 1
2 d
 [ 3 { v( x*,f* + 2 d) - v(x*,f* + d) } - { v(x*,f*+d) - v(x*,f) } ],
(34)

v(x,f)
x
 |(x,f)=(x*,f*) » 1
4 D
 [ v(x*+D,f*+d)-v(x*-D,f*+d) ] ,
(35)

2 v(x,f)
x2
 |(x,f)=(x*,f*) » 1
D2
 [ v(x*+D,f*+2 d)-2 v(x*,f*+2 d) +v(x*-D,f*+ 2 d) ] .
(36)
The molecule for the method is given in figure 9.


Figure 9. The molecule for the two-step method.
Further information the two methods is given in table 8.1 of Richtmyer and Morton [10] and in Chen [12].

6.4  Assembling and Solving the finite difference equations

Since the two finite difference steps are implicit methods, the complete solution over [0,xP] for each qi must be computed simultaneously at the dose level f*+d and then f*+2d. In the IMPETUS II program the FDMs require that the region [0,xP] is divided into a uniform grid of N nodes so that (N-1) D = xP.
By enforcing a conservation of material condition on the interval [0,D] and determining a boundary condition at xR, the finite difference equations can be assembled as a linear system
A
^
q
 
i 
 = b
(37)
where [^(q)]i = [[^(q)]i1,[^(q)]i2, ... , [^(q)]iN] with [^(q)]ij » qi((j-1)D, f) and f is the solution dose level f*+d or f*+2 d. In each step the matrix A and the vector b are determined by the finite difference approximations given in the previous subsection. In both steps the matrix A is tridiagonal and the system is solved using a method similar to the one described in section 2.5 of Mitchell and Griffiths [11] at a cost of O(N) floating point operations.
In subsection 5.1 it was shown that the PDEs were parabolic near to the surface but were predominantly hyperbolic at greater depths. Let the PDEs be regarded as predominantly parabolic in the region [0,xT). In the remaining domain of the FDM ( that is the region [xT,xP]) the hyperbolic nature of the equation can be removed by removing the u [(qi)/(x)] term, solving the remaining PDE then translating the solution by u d at each dose step. Solving the equations in this way is more error efficient than solving them directly and this technique is used in the IMPETUS II code.