Iterative Cunnigham correction in Millikans oil-drop experiment

Question

I tried to write an analysis program in python for the Millikan experiment based on a paper from Jones "The Millikan oil‐drop experiment: Making it worthwhile"

However, I cannot reproduce the calculated values of this paper for the radius of the oil drops and the charges when using an iterative Cunningham correction approach. If I iterate 6 times as indicated in the paper I get completely different results and if I take enough iteration steps, $r$ and $q$ seem to converge to zero, which doesn't make sense.

Here is a quote from the paper (with some minor changes in notation to match my program).

From Stokes' law, the radius, $r$, of the drop is given by $$ r = \sqrt{\frac{9\eta v_\text{Fall}}{2g\rho}} $$

where $\eta$ is the viscosity of the air, $v_{\text{Fall}}$ is the terminal speed of the drop falling under gravity, $g$ is the acceleration of gravity and $\rho$ is the mass density of the oil measured in air.

Stoke's law must be corrected for spheres that are comparable in size to the mean-free path of the air molecules by multiplying the viscosity by a correction factor:

$$ \eta_{\text{eff}} = \eta \left (1 + \frac{b}{pr} \right)^{-1} $$

where $\eta_{\text{eff}}$ is the effective viscosity of air $p$ is the absolute barometric pressure in $\text{cm Hg}$ and $b = 6.17\cdot 10^{-6}\,\text{m cm Hg}$ is a constant. This correction factor depends on the drop radius, and the calculated drop radius depends on this factor via the equation above. It seems to be widely assumed that that iteration between those equations is not necessary. In his original work Millikan iterated twice so that this correction converged sufficiently. Teaching versions of this experiment use lower voltages and thus often work with smaller drops than did Millikan. [...] In this study the computer iterated six times for each timing of a drop. This produced complete convergence in single-precision computer arithmetics.

Do you have any idea what may be wrong in my program:

• Is there an error in the implementation?
• Did I understand wrongly what is meant be "iteration" in the paper? (Just see my source code, of how I interpreted it)
• Why does this iteration even make sense? I don't see how to justify it mathematically. However it seems to be used by Millikan in his work which lead to a Nobel prize and also by lot's of other works about Millikan's oil drop experiment.

%pylab inline
# ipython3 notebook

def mil(v_Fall,v_Rise,U,d,rho=859.9,g=9.81,eta=1.81e-5,N=0,verbose=False):
    pi = math.pi
    b = 6.17e-6 # m cm Hg
    p = 76 # cm Hg

    # Analytic approach
    r_ana = (9*eta*v_Fall/(2*rho*g) + (b/(2*p))**2)**(0.5) - b/(2*p)
    q_ana = (6*pi*d/U)*(9*eta**3/(2*rho*g))**(0.5)* \
            (1+(b/(p*r_ana)))**(-3/2)*(v_Rise+v_Fall)*v_Fall**(0.5)

    # Iterative Approach
    r= sqrt(9*eta*v_Fall/(2*g*rho))
    for i in range(0,N):
        eta = eta*(1 + b/(p*r))**(-1)
        if verbose==True:
            print(r,eta) # it doesn't seem to converge to a value > 0
        r= sqrt(9*eta*v_Fall/(2*g*rho))
    q = 4*pi*d/(3*U)*(g*rho*r**3 + 9/2*eta*r*v_Rise)
    return array([r*1e6,r_ana*1e6,q/1.5793e-19,q_ana/1.5793e-19])

def evalMil(s_Fall,s_Rise,t_Fall,t_Rise,U,d,**kwargs):
    v_Fall = s_Fall/t_Fall
    v_Rise = s_Rise/t_Rise
    return mil(v_Fall,v_Rise,U,d,**kwargs)

# Read data
M = np.loadtxt('Millikan_Jones_test.dat', skiprows=1)
t_Fall = M[:,1]
t_Rise= M[:,0]
d = 0.0044012 # m
U = 405 # V
s = 0.0015945 # m

# eta value was not given in the paper and a bit different for each value
# (depending on the temperature)
A = evalMil(s,s,t_Fall,t_Rise,U,d,g=9.79637,N=6,eta=1.79237e-5) 
# Set N = 100 and you get q = 0 everywhere

# Analytic
rA = A[1]
rI = A[0]

# Iterative
qA = A[3]
qI = A[2]


figsize(15,10)
xlim(0,1.5)
ylim(0,10)
plot(rA,qA,'x')
plot(rI,qI,'x',color='red')
for k in range(0,35):
    axhline(y=k,color='gray')



# Output as table: 
class ListTable(list):
"""
Overridden list class which takes a 2-dimensional list of 
the form [[1,2,3],[4,5,6]], and renders an HTML Table in 
IPython Notebook.
"""

    def _repr_html_(self):
        html = ["<table>"]
        for row in self:
            html.append("<tr>")

            for col in row:
                html.append("<td>{0}</td>".format(col))

            html.append("</tr>")
        html.append("</table>")
        return ''.join(html)

tabular = ListTable()
tabular.append(['$r$ in $10^{-6}\mathrm{m}$ (iterative)',
                '$r$ in $10^{-6}\mathrm{m}$ (analytic)',
                '$q$ in $1.5793\cdot 10^{-19}\,\mathrm{C}$ (iterative)',
                '$q$ in $1.5793\cdot 10^{-19}\,\mathrm{C}$ (analytic)'])
for i in range(0,len(A[1])):
    tabular.append(A[:,i].tolist())
tabular

Here are my data as taken from the paper:

t_Rise t_Fall
9.181   17.490
22.436  52.182
35.473  11.766
30.967  12.875
29.823  12.955
9.454   7.260
4.148   5.554
23.698  17.446
28.967  7.602
7.269   12.229
11.285  7.252
7.417   15.982
7.198   17.095
15.170  34.503
28.977  45.775
26.561  46.392
10.375  47.318
10.588  2.939
6.675   10.481
14.269  11.695
24.748  50.816
23.838  51.498
4.848   7.251
8.753   15.073
18.803  28.796
88.926  30.198
36.984  12.453  

Answer 1

When solving a nonlinear equation $$ x = f(x), \qquad x\in\mathbb{R}^m, $$ the iteration method consists of generating successive guesses $$ x_1 = f(x_0), \quad x_2 = f(x_1), \qquad x_k = f(x_{k-1}). $$

If this converges to a limit $x_*$, that limit will satisfy the equation $x_* = f(x_*)$. Convergence is guaranteed by Taylor's theorem so long as $|f'|<1$ in a neighbourhood of $x_*$ and the initial guess is sufficiently close. This is discussed in every introductory numerical analysis textbook.

It is easy to explain why your solution converges to an incorrect value: iterating the equation $$ \eta = \eta \times (\cdots) $$ has only one fixed point, $\eta=0$, as the multiplier is less than 1. So this is solving the wrong equation, not the equations you wrote down.

In reality you have just one equation in one unknown $r$, but you write it as two equations. So write the two equations as $$ r = g(r, \eta_e), \qquad \eta_e = h(\eta, r). $$ Note that the two $\eta$'s are distinct: $\eta_e$ is an unknown, $\eta$ is a constant.

The fixed point method is to iterate $$ r_{k+1} = g(r_k, \eta_{e,k}), \qquad \eta_{e,k+1} = h(\eta, r_k). $$ Your code instead iterates $$ r_{k+1} = g(r_k, \eta_{e,k}), \qquad \eta_{e,k+1} = h(\eta_{e,k}, r_k), $$ which replaces the original equation for $\eta_e$ with $$ \eta_e = \eta_e / (1 + p/br). $$

Note: I haven't run your code and haven't checked for any other errors.

So your loop should read:

## Iterative Approach
r =  sqrt(9*eta*v_Fall/(2*g*rho))
for i in range(0,N):
    eta_e = eta*(1 + b/(p*r))**(-1)
    if verbose==True:
        print(r,eta_e) # now it works
    r = sqrt(9*eta_e*v_Fall/(2*g*rho))
q = 4*pi*d/(3*U)*(g*rho*r**3 + 9/2*eta_e*r*v_Rise)
return array([r*1e6, r_ana*1e6, q/1.5793e-19, q_ana/1.5793e-19])