Hello, Could we add "tolf" argument into the newton function signature? I am guessing this should match the tolf argument in IDL newton ( http://star.pst.qub.ac.uk/idl/NEWTON.html) TOLF: Set the convergence criterion on the function values. The default value is 1.0 x 10-4. def newton(func, x0, fprime=None, args=(), tol=1.48e-8, maxiter=50, tolf= 1.e-4) In second method part of the newton function (from https://github.com/scipy/scipy/blob/master/scipy/optimize/zeros.py) # Secant method p0 = x0 if x0 >= 0: p1 = x0*(1 + 1e-4) + 1e-4 else: p1 = x0*(1 + 1e-4) - 1e-4 without increasing the tolf or (1e-4) in these statements I can't get a proper root solution for my function. The reason that I am experiencing with newton is because fsolve seems slower comparing to the newton for scalar root finding for a given function. Consider this example run Sage v4.6.1 notebook: %cython cpdef double myfunc(double x): return x**3 + 2*x - 1 timeit('scipy.optimize.newton(myfunc,1)') 625 loops, best of 3: 22.1 µs per loop timeit('scipy.optimize.fsolve(myfunc,1)') 625 loops, best of 3: 86.5 µs per loop I am also experimenting to Cythonize the Newton secant method, which the Cython written version shows significant speed-ups comparing to the Python version. Without going any further, I would like to know if there is any Cythonized code around for newton or any other approach to make fsolve faster? Thanks. -- Gökhan
On Sat, May 21, 2011 at 4:24 PM, Gökhan Sever <gokhansever@gmail.com> wrote:
Hello,
Could we add "tolf" argument into the newton function signature?
I am guessing this should match the tolf argument in IDL newton ( http://star.pst.qub.ac.uk/idl/NEWTON.html)
TOLF: Set the convergence criterion on the function values. The default value is 1.0 x 10-4.
def newton(func, x0, fprime=None, args=(), tol=1.48e-8, maxiter=50, tolf= 1.e-4)
In second method part of the newton function (from https://github.com/scipy/scipy/blob/master/scipy/optimize/zeros.py)
# Secant method
p0 = x0
if x0 >= 0:
p1 = x0*(1 + 1e-4) + 1e-4
else:
p1 = x0*(1 + 1e-4) - 1e-4
without increasing the tolf or (1e-4) in these statements I can't get a proper root solution for my function.
The reason that I am experiencing with newton is because fsolve seems slower comparing to the newton for scalar root finding for a given function.
Consider this example run Sage v4.6.1 notebook:
%cython cpdef double myfunc(double x): return x**3 + 2*x - 1
timeit('scipy.optimize.newton(myfunc,1)')
625 loops, best of 3: 22.1 µs per loop
timeit('scipy.optimize.fsolve(myfunc,1)')
625 loops, best of 3: 86.5 µs per loop
I am also experimenting to Cythonize the Newton secant method, which the Cython written version shows significant speed-ups comparing to the Python version. Without going any further, I would like to know if there is any Cythonized code around for newton or any other approach to make fsolve faster?
Thanks.
You could probably adapt one of the other 1d zero finders, say ritter, just ignore all the fancy stuff for the bounding interval and such. I don't much like the stopping criterion in newton either and ftol would probably help, but it might be worth thinking about overstepping and looking for a sign change. Or something like that that would give more assurance that a zero was at hand. Note that there is also a pull request for using the second derivative as well as the first. Chuck
On Sat, May 21, 2011 at 5:25 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
You could probably adapt one of the other 1d zero finders, say ritter, just ignore all the fancy stuff for the bounding interval and such.
I have tried to use those solvers (listed under scalar functions at http://docs.scipy.org/doc/scipy/reference/optimize.html) but fail to make them work for my function. I am not sure how to automatically find an interval for the solution and make sure that f(a), f(b) will be opposite sign to each other.
I don't much like the stopping criterion in newton either and ftol would probably help, but it might be worth thinking about overstepping and looking for a sign change. Or something like that that would give more assurance that a zero was at hand. Note that there is also a pull request for using the second derivative as well as the first.
Now, I see an issue with the different versions of scipy used for newton. See my work at -> http://www.sagenb.org/home/pub/2801/ The version of the SciPy at sagenb.org is 0.7, my local sage v.4.6.1 uses v0.8 and my local SciPy is '0.10.0dev' -- a source build from a couple weeks ago. SciPy v.0.8 and above use these definitions in secant method: # Secant method if x0 >= 0: p1 = x0*(1 + 1e-4) + 1e-4 else: p1 = x0*(1 + 1e-4) - 1e-4 whereas v0.7 at sagenb.org p1 = x0*(1+1e-4) A full diff shows a bit more changes --style corrections and warnings added in 0.8 above. For the cythonized function petters_solve_for_rw "fsolve" can successfully finds the root, likewise the newton at sagenb using tol=1.e-10 argument. (newton being 5 times faster.) Later I copied the newton function from https://github.com/scipy/scipy/blob/master/scipy/optimize/zeros.py and just use the secant method parts to make it work locally. This goes to the idea of adding ftol argument and making appropriate change of 1.e-4 to 1.e-10 (at least this work well for my petters_solve_for_rw function case.) This seems slower comparing to the Scipy newton --might be due to the internals of Sage. In the next step, I tried to Cython compiled the pnewton function. This is about 10-15X faster comparing to the fsolve, about 8-9X faster than the python version (pnewton). For some reason newton in sagenb.org is fast --but 2-3 slower than Cython. One solution for myself is I can go ahead and define a local cythonized newton in my code library, unless I find another fast method for root solving from a scalar function. Indeed I use fsolve for finding two roots and passing an array rather than a scalar to make some estimations in another part of the code. However, these are a few times called calculations and don't add any overhead to the final computation. The newton that I am trying to accelerate is called about 1 to many millions times depends on the simulation case, thus the need for a fast root solver. I haven't tried using fprime option yet. How to approach this one? Use a symbolic package and estimate derivative of a function and call the solver with fprime set?
Chuck
_______________________________________________ SciPy-Dev mailing list SciPy-Dev@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-dev
-- Gökhan
On Sat, May 21, 2011 at 7:07 PM, Gökhan Sever <gokhansever@gmail.com> wrote:
On Sat, May 21, 2011 at 5:25 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
You could probably adapt one of the other 1d zero finders, say ritter,
just ignore all the fancy stuff for the bounding interval and such.
I have tried to use those solvers (listed under scalar functions at http://docs.scipy.org/doc/scipy/reference/optimize.html) but fail to make them work for my function. I am not sure how to automatically find an interval for the solution and make sure that f(a), f(b) will be opposite sign to each other.
Yeah, that left me thinking that we could really use a bracketing function. The brent optimizer must have something like that it uses to bracket minimums and perhaps that could be adapted. I know there are bracketing methods out there, IIRC there is one in NR.
I don't much like the stopping criterion in newton either and ftol would
probably help, but it might be worth thinking about overstepping and looking for a sign change. Or something like that that would give more assurance that a zero was at hand. Note that there is also a pull request for using the second derivative as well as the first.
Now, I see an issue with the different versions of scipy used for newton. See my work at -> http://www.sagenb.org/home/pub/2801/ The version of the SciPy at sagenb.org is 0.7, my local sage v.4.6.1 uses v0.8 and my local SciPy is '0.10.0dev' -- a source build from a couple weeks ago.
SciPy v.0.8 and above use these definitions in secant method:
# Secant method if x0 >= 0: p1 = x0*(1 + 1e-4) + 1e-4 else: p1 = x0*(1 + 1e-4) - 1e-4
whereas v0.7 at sagenb.org p1 = x0*(1+1e-4)
IIRC, the old version blew up when the root was at zero, the problem was posted on the list.
A full diff shows a bit more changes --style corrections and warnings added in 0.8 above.
For the cythonized function petters_solve_for_rw "fsolve" can successfully finds the root, likewise the newton at sagenb using tol=1.e-10 argument. (newton being 5 times faster.)
Later I copied the newton function from https://github.com/scipy/scipy/blob/master/scipy/optimize/zeros.py and just use the secant method parts to make it work locally. This goes to the idea of adding ftol argument and making appropriate change of 1.e-4 to 1.e-10 (at least this work well for my petters_solve_for_rw function case.) This seems slower comparing to the Scipy newton --might be due to the internals of Sage.
? I believe there is only one newton method in scipy, but we moved it into the zeros module and deprecated the version at the old location. It has since been removed.
In the next step, I tried to Cython compiled the pnewton function. This is about 10-15X faster comparing to the fsolve, about 8-9X faster than the python version (pnewton). For some reason newton in sagenb.org is fast --but 2-3 slower than Cython.
One solution for myself is I can go ahead and define a local cythonized newton in my code library, unless I find another fast method for root solving from a scalar function. Indeed I use fsolve for finding two roots and passing an array rather than a scalar to make some estimations in another part of the code. However, these are a few times called calculations and don't add any overhead to the final computation. The newton that I am trying to accelerate is called about 1 to many millions times depends on the simulation case, thus the need for a fast root solver.
I haven't tried using fprime option yet. How to approach this one? Use a symbolic package and estimate derivative of a function and call the solver with fprime set?
I would just stick to the secant method in the general case unless the derivative is easy to come by. Note that newton is one of the 'original' functions in scipy, the other 1d zero finders came later. So if you can make an improved cythonized version I don't see any reason not to use it. If you do make cythonized versions it might be worth implementing the Newton and secant parts separately and make the current newton a driver function. How many function evaluations are you seeing in the root finding? Chuck
On Sat, May 21, 2011 at 7:50 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
Yeah, that left me thinking that we could really use a bracketing function. The brent optimizer must have something like that it uses to bracket minimums and perhaps that could be adapted. I know there are bracketing methods out there, IIRC there is one in NR.
I have managed to get those bracketed solvers working after plotting my function. brenth seems converging the fastest but still providing the search interval is a problem for me. I am probably skipping these solvers.
IIRC, the old version blew up when the root was at zero, the problem was posted on the list.
? I believe there is only one newton method in scipy, but we moved it into the zeros module and deprecated the version at the old location. It has since been removed.
Yes, you are right. At the current source repository, there is only one newton in scipy.optimize which resides in zeros.py
I would just stick to the secant method in the general case unless the derivative is easy to come by. Note that newton is one of the 'original' functions in scipy, the other 1d zero finders came later. So if you can make an improved cythonized version I don't see any reason not to use it. If you do make cythonized versions it might be worth implementing the Newton and secant parts separately and make the current newton a driver function.
I have figured out the derivative option and tested fsolve and newton with fprime arg provided. Still slower comparing to the secant method. Most likely, that the derivative function requires a bit calculation to be evaluated. As you can see, the function is: 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)) but the derivative is quite complex comparing to the function: cpdef double pprime(double x, double rd, double rh): return -3*(rd**3 - x**3)*x**2*exp(kelvin/x)/((kappa - 1.0)*rd**3 +x**3)**2 - 3*x**2.*exp(kelvin/x)/((kappa - 1.0)*rd**3 + x**3)-\ (rd**3 - x**3)*kelvin*exp(kelvin/x)/(((kappa - 1.0)*rd**3 + x**3)*x**2) Skipping the fprime option, I focus on the secant method which works the fastest and quite robust for my case. Below you can see the latest version of the cythonized secant (probably I should update the name) that I use: cpdef double cnewton(func, double x0, args=(), double tol=1e-10, int maxiter=50): # Secant method p0 = x0 p1 = x0*(1 + 1e-10) + 1e-10 #p1 = x0*(1+1e-4) q0 = func(*((p0,) + args)) q1 = func(*((p1,) + args)) for iter in range(maxiter): p = p1 - q1*(p1 - p0)/(q1 - q0) if abs(p - p1) < tol: return p p0 = p1 q0 = q1 p1 = p q1 = func(*((p1,) + args)) I simplified this block if x0 >= 0: p1 = x0*(1 + 1e-4) + 1e-4 else: p1 = x0*(1 + 1e-4) - 1e-4 as just, one line since I am not interested with negative init point. In other words, reals drops can never be less than < 0 meters. p1 = x0*(1 + 1e-10) + 1e-10 #p1 = x0*(1+1e-4) I also skipped this block: if q1 == q0: if p1 != p0: msg = "Tolerance of %s reached" % (p1 - p0) warnings.warn(msg, RuntimeWarning) return (p1 + p0)/2.0 because, this part is almost never executed.
How many function evaluations are you seeing in the root finding?
With this version of the newton, I see about 9-10 function evaluations to converge to a meaningful and reasonable root.
Chuck
On Sun, May 22, 2011 at 1:55 AM, Gökhan Sever <gokhansever@gmail.com> wrote:
On Sat, May 21, 2011 at 7:50 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
Yeah, that left me thinking that we could really use a bracketing function. The brent optimizer must have something like that it uses to bracket minimums and perhaps that could be adapted. I know there are bracketing methods out there, IIRC there is one in NR.
I have managed to get those bracketed solvers working after plotting my function. brenth seems converging the fastest but still providing the search interval is a problem for me. I am probably skipping these solvers.
IIRC, the old version blew up when the root was at zero, the problem was posted on the list.
? I believe there is only one newton method in scipy, but we moved it into the zeros module and deprecated the version at the old location. It has since been removed.
Yes, you are right. At the current source repository, there is only one newton in scipy.optimize which resides in zeros.py
I would just stick to the secant method in the general case unless the derivative is easy to come by. Note that newton is one of the 'original' functions in scipy, the other 1d zero finders came later. So if you can make an improved cythonized version I don't see any reason not to use it. If you do make cythonized versions it might be worth implementing the Newton and secant parts separately and make the current newton a driver function.
I have figured out the derivative option and tested fsolve and newton with fprime arg provided. Still slower comparing to the secant method. Most likely, that the derivative function requires a bit calculation to be evaluated. As you can see, the function is:
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))
Wouldn't this be easier (for derivatives), and maybe be more stable, taking logs np.log(rh) - kelvin/x + np.log(..) ... ? (independently of any improvement to the solvers) Josef
On Sun, May 22, 2011 at 6:22 AM, <josef.pktd@gmail.com> wrote:
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))
Wouldn't this be easier (for derivatives), and maybe be more stable, taking logs
np.log(rh) - kelvin/x + np.log(..) ... ?
(independently of any improvement to the solvers)
Seems like this produces more terms in derivatives (tested below in Sage v.4.6.1 via notebook): myfunc (rd^3 - x^3)*e^(kelvin/x)/((kappa - 1.0)*rd^3 + x^3) + rh myfunc.derivative(x).simplify() -3*(rd^3 - x^3)*x^2*e^(kelvin/x)/((kappa - 1.0)*rd^3 + x^3)^2 - 3*x^2*e^(kelvin/x)/((kappa - 1.0)*rd^3 + x^3) - (rd^3 - x^3)*kelvin*e^(kelvin/x)/(((kappa - 1.0)*rd^3 + x^3)*x^2) p = myfunc.log() p.derivative(x).simplify() -(3*(rd^3 - x^3)*x^2*e^(kelvin/x)/((kappa - 1.0)*rd^3 + x^3)^2 + 3*x^2*e^(kelvin/x)/((kappa - 1.0)*rd^3 + x^3) + (rd^3 - x^3)*kelvin*e^(kelvin/x)/(((kappa - 1.0)*rd^3 + x^3)*x^2))/((rd^3 - x^3)*e^(kelvin/x)/((kappa - 1.0)*rd^3 + x^3) + rh) -- Gökhan
On Sun, May 22, 2011 at 12:59 PM, Gökhan Sever <gokhansever@gmail.com> wrote:
On Sun, May 22, 2011 at 6:22 AM, <josef.pktd@gmail.com> wrote:
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))
Wouldn't this be easier (for derivatives), and maybe be more stable, taking logs
np.log(rh) - kelvin/x + np.log(..) ... ?
(independently of any improvement to the solvers)
Seems like this produces more terms in derivatives (tested below in Sage v.4.6.1 via notebook):
myfunc (rd^3 - x^3)*e^(kelvin/x)/((kappa - 1.0)*rd^3 + x^3) + rh
myfunc.derivative(x).simplify() -3*(rd^3 - x^3)*x^2*e^(kelvin/x)/((kappa - 1.0)*rd^3 + x^3)^2 - 3*x^2*e^(kelvin/x)/((kappa - 1.0)*rd^3 + x^3) - (rd^3 - x^3)*kelvin*e^(kelvin/x)/(((kappa - 1.0)*rd^3 + x^3)*x^2)
p = myfunc.log()
I proposed taking logs of left hand side and right hand side separately, since you are just looking for a zero, with myfunc.log(), it is not simplified (I don't have a quick way to do the symbolic derivative, but there shouldn'd be any exp left in the expression) Josef
p.derivative(x).simplify() -(3*(rd^3 - x^3)*x^2*e^(kelvin/x)/((kappa - 1.0)*rd^3 + x^3)^2 + 3*x^2*e^(kelvin/x)/((kappa - 1.0)*rd^3 + x^3) + (rd^3 - x^3)*kelvin*e^(kelvin/x)/(((kappa - 1.0)*rd^3 + x^3)*x^2))/((rd^3 - x^3)*e^(kelvin/x)/((kappa - 1.0)*rd^3 + x^3) + rh)
-- Gökhan _______________________________________________ SciPy-Dev mailing list SciPy-Dev@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-dev
On Sun, May 22, 2011 at 11:12 AM, <josef.pktd@gmail.com> wrote:
I proposed taking logs of left hand side and right hand side separately, since you are just looking for a zero, with myfunc.log(), it is not simplified
OK, I have gotten this right this time: myfunc = rh - exp(kelvin/x) * (x**3 - rd**3) / (x**3 - rd**3 * (1.0 - kappa)) myfunc2 = log(myfunc) = log(rh) - (kelvin/x) + log(x**3 - rd**3) - log(x**3 - rd**3*(1.0 - kappa)) myfunc2_prime = -3*x**2/(rd**3 - x**3) - 3*x**2/((kappa - 1.0)*rd**3 + x**3) + kelvin/x**2 How can I proceed this point onwards?
(I don't have a quick way to do the symbolic derivative, but there shouldn'd be any exp left in the expression)
On Sun, May 22, 2011 at 4:24 PM, Gökhan Sever <gokhansever@gmail.com> wrote:
On Sun, May 22, 2011 at 11:12 AM, <josef.pktd@gmail.com> wrote:
I proposed taking logs of left hand side and right hand side separately, since you are just looking for a zero, with myfunc.log(), it is not simplified
OK, I have gotten this right this time:
myfunc = rh - exp(kelvin/x) * (x**3 - rd**3) / (x**3 - rd**3 * (1.0 - kappa)) myfunc2 = log(myfunc) = log(rh) - (kelvin/x) + log(x**3 - rd**3) - log(x**3 - rd**3*(1.0 - kappa))
myfunc2_prime = -3*x**2/(rd**3 - x**3) - 3*x**2/((kappa - 1.0)*rd**3 + x**3) + kelvin/x**2
How can I proceed this point onwards?
try newton with fprime. My initial suggestion was in response to your statement that newton with fprime is too slow because the expression for the derivative is too complicated and slow. Trying to get the function in a nicer form might help quite a bit, but it won't be a solution if you have a large set of functions that might show up in different simulations. It's just an aside for the main topic of the thread, improving the solvers. Josef
(I don't have a quick way to do the symbolic derivative, but there shouldn'd be any exp left in the expression)
_______________________________________________ SciPy-Dev mailing list SciPy-Dev@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-dev
On Sun, May 22, 2011 at 2:50 PM, <josef.pktd@gmail.com> wrote:
try newton with fprime. My initial suggestion was in response to your statement that newton with fprime is too slow because the expression for the derivative is too complicated and slow.
This seems to be not working for my case: cpdef double myfunc2(double x, double rd, double rh): return log(rh) - (kelvin/x) + log(x**3 - rd**3) - log(x**3 - rd**3*(1.0 - kappa)) cpdef double myfunc2_prime(double x, double rd, double rh): -3*x**2/(rd**3 - x**3) - 3*x**2/((kappa - 1.0)*rd**3 + x**3) + kelvin/x**2 rd = 5.75e-08; rh = 0.95 I[4]: newton(myfunc2, rd, args=(rd, rh), fprime=myfunc2_prime) /usr/lib64/python2.7/site-packages/scipy/optimize/zeros.py:106: RuntimeWarning: derivative was zero. warnings.warn(msg, RuntimeWarning) O[4]: 5.75e-08 # eliminating the zero derivative. I[5]: newton(myfunc2, rd, args=(rd, rh), fprime=myfunc2_prime, tol=1.e-10) O[5]: 5.75e-08 The correct result is: Setting tol to different accuracies makes a different. In this case tol=1.e-20 yields exact solution. I[7]: cnewton(petters_solve_for_rw, rd, args=(rd, rh), tol=1.e-20) O[7]: 1.4972782377152967e-07
Trying to get the function in a nicer form might help quite a bit,
Yes, I can confirm this from Charles Harris' suggestion, bit as I said the speed gain isn't that significant in this case and readability still counts. but
it won't be a solution if you have a large set of functions that might show up in different simulations.
It's just an aside for the main topic of the thread, improving the solvers.
Josef
-- Gökhan
On Sat, May 21, 2011 at 11:55 PM, Gökhan Sever <gokhansever@gmail.com>wrote:
On Sat, May 21, 2011 at 7:50 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
Yeah, that left me thinking that we could really use a bracketing function. The brent optimizer must have something like that it uses to bracket minimums and perhaps that could be adapted. I know there are bracketing methods out there, IIRC there is one in NR.
I have managed to get those bracketed solvers working after plotting my function. brenth seems converging the fastest but still providing the search interval is a problem for me. I am probably skipping these solvers.
IIRC, the old version blew up when the root was at zero, the problem was posted on the list.
? I believe there is only one newton method in scipy, but we moved it
into
the zeros module and deprecated the version at the old location. It has since been removed.
Yes, you are right. At the current source repository, there is only one newton in scipy.optimize which resides in zeros.py
I would just stick to the secant method in the general case unless the derivative is easy to come by. Note that newton is one of the 'original' functions in scipy, the other 1d zero finders came later. So if you can
make
an improved cythonized version I don't see any reason not to use it. If you do make cythonized versions it might be worth implementing the Newton and secant parts separately and make the current newton a driver function.
I have figured out the derivative option and tested fsolve and newton with fprime arg provided. Still slower comparing to the secant method. Most likely, that the derivative function requires a bit calculation to be evaluated. As you can see, the function is:
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))
You could also try rewriting this in various ways. For instance cpdef double petters_solve_for_rw(double x, double rd, double rh): return rh*( x**3 - rd**3 * (1.0 - kappa)) - exp(kelvin/x) * (x**3 - rd**3) or cpdef double petters_solve_for_rw(double x, double rd, double rh): return rh*(1 - (1 - kappa)*y**3) - exp(y*kelvin/rd) * (1 - y**3) or cpdef double petters_solve_for_rw(double x, double rd, double rh): return rh*kappa*y**3 - (exp(y*kelvin/rd) - rh) * (1 - y**3) where x = rd/y. The last might allow you to bracket things fairly easily, i.e., (exp(y*kelvin/rd) - rh) * (1 - y**3) has to be >0 if you expect y>0 <snip> Chuck
On Sun, May 22, 2011 at 12:00 PM, Charles R Harris <charlesr.harris@gmail.com>
cpdef double petters_solve_for_rw(double x, double rd, double rh): return rh*kappa*y**3 - (exp(y*kelvin/rd) - rh) * (1 - y**3)
This last modification is converging faster than the original version, but readability of the function is reduced now.
where x = rd/y. The last might allow you to bracket things fairly easily, i.e., (exp(y*kelvin/rd) - rh) * (1 - y**3) has to be >0 if you expect y>0
"rh" also plays role in determining the sign of the right portion of this equation. Throughout the model rh usually goes from 0.95 and pass beyond 1.0. This causes a sign change.
<snip>
Chuck
_______________________________________________ SciPy-Dev mailing list SciPy-Dev@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-dev
-- Gökhan
On Sun, May 22, 2011 at 2:46 PM, Gökhan Sever <gokhansever@gmail.com> wrote:
On Sun, May 22, 2011 at 12:00 PM, Charles R Harris < charlesr.harris@gmail.com>
cpdef double petters_solve_for_rw(double x, double rd, double rh): return rh*kappa*y**3 - (exp(y*kelvin/rd) - rh) * (1 - y**3)
This last modification is converging faster than the original version, but readability of the function is reduced now.
where x = rd/y. The last might allow you to bracket things fairly easily, i.e., (exp(y*kelvin/rd) - rh) * (1 - y**3) has to be >0 if you expect y>0
"rh" also plays role in determining the sign of the right portion of this equation. Throughout the model rh usually goes from 0.95 and pass beyond 1.0. This causes a sign change.
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? Chuck
On Sun, May 22, 2011 at 2:51 PM, Charles R Harris <charlesr.harris@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.
On Sun, May 22, 2011 at 3:21 PM, Gökhan Sever <gokhansever@gmail.com> wrote:
On Sun, May 22, 2011 at 2:51 PM, Charles R Harris <charlesr.harris@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
On Sun, May 22, 2011 at 6:04 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
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.
rh = 1 is a special case, right when the supersaturation is reached within the parcel model. I highly suspect that we get an exact rh=1 throughout the simulations. This value is usually rh=1+-small number. However I might need to verify this further. Soon, I will work on separating the model thermodynamics for rh<1 and rh>1 cases which I will have to estimate the closest rh=1 point where the saturation starts occurring.
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].
How did you get this one for rh<1? What happened to the exponential term? Even if so, how I am going to ensure that for the [rd, b] interval f(rd) and f(b) will result opposite sign results?
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.
x is >= 0. I don't use negative x as an initial estimator. Neither the function or the root solver should yield a negative result.
Using the reciprical of x can help as it makes the exponential continuous for x = +/- inf.
Chuck
_______________________________________________ SciPy-Dev mailing list SciPy-Dev@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-dev
-- Gökhan
On Sun, May 22, 2011 at 8:05 PM, Gökhan Sever <gokhansever@gmail.com> wrote:
On Sun, May 22, 2011 at 6:04 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
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.
rh = 1 is a special case, right when the supersaturation is reached within the parcel model. I highly suspect that we get an exact rh=1 throughout the simulations. This value is usually rh=1+-small number. However I might need to verify this further. Soon, I will work on separating the model thermodynamics for rh<1 and rh>1 cases which I will have to estimate the closest rh=1 point where the saturation starts occurring.
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].
How did you get this one for rh<1? What happened to the exponential term? Even if so, how I am going to ensure that for the [rd, b] interval f(rd) and f(b) will result opposite sign results?
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.
x is >= 0. I don't use negative x as an initial estimator. Neither the function or the root solver should yield a negative result
I'm not so confidant about the rh >= 1 case, but I've attached an example for rh = .95, rd=1e-8. The light blue line is the lhs from above, the labled lines are for the rhs and different values of kappa. The heavy horizontal line is the rh and the bracket I was suggesting was between the zero of the rhs at x=rd and its crossing with the rh line. There are corner cases here depending on the parameter values, so this probably needs more exploration, there might be cases with two zeros. Also same thing with rd=1.5e-6. Note that the zero is very near the upper limit. I suspect there will none or two zeros in the supersaturated case. Chuck
On Sun, May 22, 2011 at 9:13 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
I'm not so confidant about the rh >= 1 case, but I've attached an example for rh = .95, rd=1e-8. The light blue line is the lhs from above, the labled lines are for the rhs and different values of kappa. The heavy horizontal line is the rh and the bracket I was suggesting was between the zero of the rhs at x=rd and its crossing with the rh line. There are corner cases here depending on the parameter values, so this probably needs more exploration, there might be cases with two zeros.
Also same thing with rd=1.5e-6. Note that the zero is very near the upper limit.
I suspect there will none or two zeros in the supersaturated case.
Chuck
Thanks for spending your time and producing those plots. Could you please provide the code that you used to create the figures? This might help me to better understand some of the points you have made in your latest reply. Your comment on having no or two zero worries me a bit. I can't easily see how apart these two zeros from each if they ever exist. One of them could be unrealistic to easily disregard, but yet to verify this claim. I went ahead and tested the secant and fsolve solvers to see if they produce any significant differences in terms of the root they return. I assume fsolve is a more robust solver. You can see this comparison in the attached figure. I use a tolerance value of 1.e-20 for both solvers. Again this comparison is based on the estimate of about 20k different values. Most of the difference is zero. I focused in to a more interesting part of the figure. Still the difference is about 1.e-17 which is quite insignificant.
On Mon, May 23, 2011 at 12:21 AM, Gökhan Sever <gokhansever@gmail.com>wrote:
I'm not so confidant about the rh >= 1 case, but I've attached an example for rh = .95, rd=1e-8. The light blue line is the lhs from above, the labled lines are for the rhs and different values of kappa. The heavy horizontal line is the rh and the bracket I was suggesting was between the zero of
On Sun, May 22, 2011 at 9:13 PM, Charles R Harris <charlesr.harris@gmail.com> wrote: the
rhs at x=rd and its crossing with the rh line. There are corner cases here depending on the parameter values, so this probably needs more exploration, there might be cases with two zeros.
Also same thing with rd=1.5e-6. Note that the zero is very near the upper limit.
I suspect there will none or two zeros in the supersaturated case.
Chuck
Thanks for spending your time and producing those plots. Could you please provide the code that you used to create the figures? This might help me to better understand some of the points you have made in your latest reply.
I've attached the module with the lhs, rhs functions. I hope I got them right ;) The plots were done using x = linspace(small number > 0, 5*rd, 500) and a loop for the values of kappa. They might actually look better as a semilogx plot.
Your comment on having no or two zero worries me a bit. I can't easily see how apart these two zeros from each if they ever exist. One of them could be unrealistic to easily disregard, but yet to verify this claim.
The reason I think there will be two zeros is that the upper branch of the hyperbolic rhs is concave up while the lhs is concave down, so if they intersect it will be at two points or a tangent (double zero). At least the supersaturated case shows up as being a bit squirrelly, which is probably a good sign for the model ;) Also note that the zeros go off to +inf as rh -> 1, which might be a good argument for using rd/x as the independent variable.
I went ahead and tested the secant and fsolve solvers to see if they produce any significant differences in terms of the root they return. I assume fsolve is a more robust solver. You can see this comparison in the attached figure. I use a tolerance value of 1.e-20 for both solvers. Again this comparison is based on the estimate of about 20k different values. Most of the difference is zero. I focused in to a more interesting part of the figure. Still the difference is about 1.e-17 which is quite insignificant.
I think you are right that the secant solver needs a user input for the initial step size. Although working near singularities is always going to be a problem, hence variable changes. In fact, the whole newton thing could probably use a think through. Chuck
On Sun, May 22, 2011 at 6:04 PM, Charles R Harris <charlesr.harris@gmail.com
wrote:
On Sun, May 22, 2011 at 3:21 PM, Gökhan Sever <gokhansever@gmail.com>wrote:
On Sun, May 22, 2011 at 2:51 PM, Charles R Harris <charlesr.harris@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.
This isn't quite right it seems.
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
participants (3)
-
Charles R Harris -
Gökhan Sever -
josef.pktd＠gmail.com