[SciPy-Dev] Comments on optimize.newton function

[Charles R Harris]

On Sun, May 22, 2011 at 3:21 PM, Gökhan Sever <gokhansever at gmail.com> wrote:

> On Sun, May 22, 2011 at 2:51 PM, Charles R Harris
> <charlesr.harris at gmail.com> wrote:
> > I think the zeros of this function can be bracketed by inspection. What
> sort
> > of values do rd, rh, and kappa have? What is kelvin?
>
> This is the original function:
>
> 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))
>
> "kelvin" is a constant: 1.04962912337e-09 and stays constant
> throughout all of the simulations.
>
> "kappa" is a constant, but its value set before the simulation.
> Default is 1, but can range from 0.001 to 2 depend on the simulation
>
> "rd" is initialized differently. For one simulation about 20k element
> rd array created --this number changes depends on the simulation
> --solving a set of five ODE equations. For one case:
>
> I[3]: rd.min()
> O[3]: 1.1926858018899999e-08
>
> I[4]: rd.max()
> O[4]: 1.3455000000000001e-06
>
> "rh" is the relative humidity. Starts at rh=0.95, and evolves like
> "rd" within the simulation, and differs from simulation to simulation
> depends on the initial conditions.
>
> I[9]: rh.max()
> O[9]: 1.0050122345200001
>
> I[10]: rh.min()
> O[10]: 0.95017287164200004
>
> With these numbers, I still think it is hard to bracket this function
> within which a root is searched.


Solve

rh*exp(-kelvin/x) = (x**3 - rd**3) / (x**3 - rd**3 * (1.0 - kappa))

The lhs increases from 0 to rh as x -> inf, hence is <= rh.
The rhs looks sort like a hyperbola with a horizontal asymptote at y = 1,
and a vertical asymptote at rd*(1 - kappa)**1/3.

If rh = 1, there is no solution unless kappa = 0 and x = +/- inf. I suspect
that might be a problem and a hint that the model might be a bit off.

If rh < 1, solve rh = (b**3 - rd**3) / (b**3 - rd**3 * (1.0 - kappa)) for b,
which you can do algebraically, and the root will lie in the interval [rd,
b].

If rh > 1, things are a mess, but x < 0 and also to the left of the vertical
asymptote, and to the right of b solved for previously from above. Is the
negative x a problem? There is no (real) solution in this case if 1/(1 -
kappa)  < rh, and more generally, if the bracket doesn't contain any values.

Using the reciprical of x can help as it makes the exponential continuous
for x = +/- inf.

Chuck
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/scipy-dev/attachments/20110522/1ab842bb/attachment.html>