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 [22]: 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.

Tue Sep 17 16:40:38 BST 1996