<br><br><div class="gmail_quote">On Sun, May 22, 2011 at 3:21 PM, Gökhan Sever <span dir="ltr"><<a href="mailto:gokhansever@gmail.com">gokhansever@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin: 0pt 0pt 0pt 0.8ex; border-left: 1px solid rgb(204, 204, 204); padding-left: 1ex;">
<div class="im">On Sun, May 22, 2011 at 2:51 PM, Charles R Harris<br>
<<a href="mailto:charlesr.harris@gmail.com">charlesr.harris@gmail.com</a>> wrote:<br>
> I think the zeros of this function can be bracketed by inspection. What sort<br>
> of values do rd, rh, and kappa have? What is kelvin?<br>
<br>
</div>This is the original function:<br>
<div class="im"><br>
cpdef double petters_solve_for_rw(double x, double rd, double rh):<br>
</div><div class="im"> return rh - exp(kelvin/x) * (x**3 - rd**3) / (x**3 - rd**3 * (1.0 - kappa))<br>
<br>
</div>"kelvin" is a constant: 1.04962912337e-09 and stays constant<br>
throughout all of the simulations.<br>
<br>
"kappa" is a constant, but its value set before the simulation.<br>
Default is 1, but can range from 0.001 to 2 depend on the simulation<br>
<br>
"rd" is initialized differently. For one simulation about 20k element<br>
rd array created --this number changes depends on the simulation<br>
--solving a set of five ODE equations. For one case:<br>
<br>
I[3]: rd.min()<br>
O[3]: 1.1926858018899999e-08<br>
<br>
I[4]: rd.max()<br>
O[4]: 1.3455000000000001e-06<br>
<br>
"rh" is the relative humidity. Starts at rh=0.95, and evolves like<br>
"rd" within the simulation, and differs from simulation to simulation<br>
depends on the initial conditions.<br>
<br>
I[9]: rh.max()<br>
O[9]: 1.0050122345200001<br>
<br>
I[10]: rh.min()<br>
O[10]: 0.95017287164200004<br>
<br>
With these numbers, I still think it is hard to bracket this function<br>
within which a root is searched.</blockquote><div><br>Solve<br><br>rh*exp(-kelvin/x) = (x**3 - rd**3) / (x**3 - rd**3 * (1.0 - kappa))<br><br>The lhs increases from 0 to rh as x -> inf, hence is <= rh.<br>The rhs looks sort like a hyperbola with a horizontal asymptote at y = 1, and a vertical asymptote at rd*(1 - kappa)**1/3.<br>
<br>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. <br><br>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].<br>
<br>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.<br>
<br>Using the reciprical of x can help as it makes the exponential continuous for x = +/- inf.<br><br>Chuck <br></div></div>