2801.sagews

%cython
cimport cython

cdef extern from "math.h":
    double exp(double)

cdef double gascon = 8.3144
cdef double mw_h2o = 0.018015
cdef double gascon_h2o = gascon / mw_h2o
cdef double sigma_at_tk = 0.072
cdef double rho_h2o = 997.0
cdef double tk_aerosol = 298.15
cdef double kappa = 1.0
cdef double kelvin = 2.0 * sigma_at_tk / (gascon_h2o * rho_h2o * tk_aerosol)
                                            
@cython.cdivision(True)
cpdef double petters_solve_for_rw(double x, double rd, double rh):
    return rh - exp(kelvin/x) * (x**3 - rd**3) / (x**3 - rd**3 * (1.0 - kappa))

from scipy import optimize
# These are input to the petters_solve_for_rw
rd = 5.75e-08
rh = 0.95

# First, find a solution
optimize.fsolve(petters_solve_for_rw, rd, args=(rd, rh))

1.4972782377152967e-07

timeit('optimize.fsolve(petters_solve_for_rw, rd, args=(rd, rh))')

625 loops, best of 3: 133 µs per loop

# my local scipy v0.10.dev can't produce this result
optimize.newton(petters_solve_for_rw, rd, args=(rd, rh), tol=1.e-10)

1.49726676560525e-7

timeit('optimize.newton(petters_solve_for_rw, rd, args=(rd, rh))')

625 loops, best of 3: 26.2 µs per loop

#SAGE scipy version
from scipy import __version__
print __version__
#newton function seems changed in between versions.

0.7.0

# Slightly modified newton -- 1.e-4 changed to 1.e-10
# Newton --using the secant method
def pnewton(func, x0, args=(), tol=1.e-10, maxiter=50):
    """Find a zero using the secant method.

    Find a zero of the function `func` given a nearby starting point `x0`.

    Parameters
    ----------
    func : function
        The function whose zero is wanted. It must be a function of a
        single variable of the form f(x,a,b,c...), where a,b,c... are extra
        arguments that can be passed in the `args` parameter.
    x0 : float
        An initial estimate of the zero that should be somewhere near the
        actual zero.
    args : tuple, optional
        Extra arguments to be used in the function call.
    tol : float, optional
        The allowable error of the zero value.
    maxiter : int, optional
        Maximum number of iterations.

    Returns
    -------
    zero : float
        Estimated location where function is zero.
    """
    # Secant method
    p0 = x0
    if x0 >= 0:
        p1 = x0*(1 + 1e-10) + 1e-10
    else:
        p1 = x0*(1 + 1e-10) - 1e-10
    q0 = func(*((p0,) + args))
    q1 = func(*((p1,) + args))
    for iter in range(maxiter):
        if q1 == q0:
            if p1 != p0:
                return (p1 + p0)/2.0
        else:
            p = p1 - q1*(p1 - p0)/(q1 - q0)
        if abs(p - p1) < tol:
            return p
        p0 = p1
        q0 = q1
        p1 = p
        q1 = func(*((p1,) + args))

# Now the result is correct.
pnewton(petters_solve_for_rw, rd, args=(rd, rh))

1.49726708556492e-7

timeit('pnewton(petters_solve_for_rw, rd, args=(rd, rh))')

625 loops, best of 3: 79.3 µs per loop

%cython
# An effort to cythonize the pnewton()
cpdef double cnewton(func, double x0, args=(), double tol=1.e-10, int maxiter=50):
    """Find a zero using the secant method.

    Find a zero of the function `func` given a nearby starting point `x0`.

    Parameters
    ----------
    func : function
        The function whose zero is wanted. It must be a function of a
        single variable of the form f(x,a,b,c...), where a,b,c... are extra
        arguments that can be passed in the `args` parameter.
    x0 : float
        An initial estimate of the zero that should be somewhere near the
        actual zero.
    args : tuple, optional
        Extra arguments to be used in the function call.
    tol : float, optional
        The allowable error of the zero value.
    maxiter : int, optional
        Maximum number of iterations.

    Returns
    -------
    zero : float
        Estimated location where function is zero.
    """
    # Secant method
    p0 = x0
    if x0 >= 0:
        p1 = x0*(1 + 1e-10) + 1e-10
    else:
        p1 = x0*(1 + 1e-10) - 1e-10
    q0 = func(*((p0,) + args))
    q1 = func(*((p1,) + args))
    for iter in range(maxiter):
        if q1 == q0:
            if p1 != p0:
                return (p1 + p0)/2.0
        else:
            p = p1 - q1*(p1 - p0)/(q1 - q0)
        if abs(p - p1) < tol:
            return p
        p0 = p1
        q0 = q1
        p1 = p
        q1 = func(*((p1,) + args))

cnewton(petters_solve_for_rw, rd, args=(rd, rh))

1.4972670855649154e-07

timeit('cnewton(petters_solve_for_rw, rd, args=(rd, rh))')

625 loops, best of 3: 9.05 µs per loop