Also shown in Fig.3,4 are the numerical results for the full self-consistent integral equation, (14). This we have solved iteratively after the fashion of Broyles : an initial guess input function is used in computing the right hand side of Eq.(14) giving an output function , and the next input function is chosen to be:
where is a constant chosen between 0 and 1 for optimal convergence of the process. For computational purposes, the function has to be specified on a two dimensional grid (assuming unbroken symmetry about the axis normal to the plates), with h being a fixed parameter during the iterations. The function is obtained by fitting bi-cubic splines between grid points. A fine grid means more grid points where evaluations of the function n has to be made, whereas a coarse one would not supply a satisfactory accuracy (which we define to be convergence of the force to within of the maximum value of for any given force curve at fixed ).
The bottleneck of the iterative process is the integral contained in , which is of high dimension, and makes precise evaluations very difficult. This problem is alleviated by splitting the integral into two more straightforward ones:
where , , and f and g are two dimensionless functions:
Substituting f and g into Eq.(34), together with the relations between 's, gives the original Eq.(15). Note that f(t,t') can be calculated and tabulated beforehand to speed up the iteration procedure, and is a simple integral with its integrand specifying a straight line segment in the plane.
It turns out the number of iterations needed to converge Eq.(33) is of the order 5-15 with between 0.1 and 0.5, and rather insensitive to the initial guess function , which is therefore always set to the constant bulk value . Not surprisingly, the process is a lot faster for smaller concentrations ( ), and it gets increasingly difficult when approaches 2, which we found to be the maximum value for which converged results could be obtained in reasonable time (5-10 hrs CPU per value on a DEC Alpha 600/5/266). In order to have the fastest convergence, a two dimensional adaptive grid was adopted which becomes finer wherever the function n varies most rapidly.
After numerically solving Eq.(14) for , and following Eq.(4), we can find:
We found that this leads, with a high degree of numerical accuracy and irrespective of the value of to:
so that the resulting pressure is once more identical to the second virial result [21, 18], which is known to be exact in the limit of which we have adopted.
The depletion force for general separation h is then given by equations (3, 37) as:
As in the perturbative treatments of Section 3, the force between spheres of radius R is then found by the Derjaguin integral, Eq.(1); the interaction energy between two such spheres can then be found by a further integration.