I hate to bash on this article, but the number of commenters calling it "nice" makes me want to warn that the implementations are not really very good.
Except for the bisection method, all of these implementations take an argument specifying the number of iterations to run. In most cases, the only way to terminate in fewer iterations is by hitting an "exact" root, i.e., calculating the residual as exactly zero. This is poor practice for a number of reasons. First, in practice it's pretty rare for a method to find an exact zero. Second, once a method has converged to the numerical precision of the machine, making more iterations just wastes flops. So a much better approach is to specify a solution tolerance (as shown with the bisection method). Even better is to provide absolute and relative tolerances, and to choose those values based on either the domain requirements, or on the machine characteristics. Dennis & Schnabel's excellent "Numerical Methods for Unconstrained Optimization and Nonlinear Equations" has a good discussion on choosing convergence tolerances.
This dependence on iteration counts to terminate, by the way, is probably why the author equates low iteration counts with greater accuracy. But in fact these methods don't vary in their intrinsic accuracy, rather, they vary in their order of convergence.
Another example of poor practice is in the bisection method implementation. One generally should not bisect an interval using c = (a+b)/2, because the nature of finite-precision arithmetic means there is no guarantee that c will lie between a and b, even if the machine can represent numbers between a and b. A better approach is to ensure a < b, then to set c = a + (b-a)/2. This expression is much less subject to roundoff errors.
All good valid points. The article gives a high overview but lacks on the detail. All these points are valid and some mentioned in the article. The author mentioned a few times that the for-loops should be replaced with while-loops based on convergence tolerance. As for the bisection interval, I am not sure the basic c=(a+b)/2 is not enough. Can you provide an example where it fails?
As for the bisection interval, I am not sure the basic c=(a+b)/2 is not enough. Can you provide an example where it fails?
Sure. I've experienced this myself, for example, when doing the research cited here:
http://pubs.acs.org/doi/abs/10.1021/es202359j
Unfortunately, my Lab doesn't have our report version of the work up on the web yet (I believe it should be freely available at some point, but they're re-doing the library).
In addition, I knew about this particular effect before I started that work, so I must have run into it earlier.
In most cases, naive bisection works fine. That's especially true if you use pretty loose tolerances-- for example, you if only want to identify the zero to a relative tolerance of 1e-6 or something coarse like that. But this work demanded very tight tolerances (because the output gets used in a time-marching numerical simulation that may iterate on time steps, and hence has to be very repeatable no matter what the initial condition). In this work, the relative tolerance was about 2e-15, and the absolute tolerance about 1e-18 (on a machine with double-precision epsilon about 2e-16).
By the way, if you're testing for convergence according to a bracket on the value, then you want to know (b-a) anyhow. So there's no reason not to go ahead and use it to bisect the interval as well.
Another interesting implementation detail-- at very tight relative tolerances, you also need to calculate the terminating bracket width based on the larger of |b| and |a|-- due to floating point effects, if you find the terminating bracket width by applying the relative tolerance to the smaller value, you can occasionally run into a problem where adding that width to one of a or b fails to change the estimated value.
Wouldn't you consider that a very extreme situation? Also, I am not aware of a MAX_INT in python. As far as I know it is the limits of your memory. If you are trying to work with such big numbers you should really use something a little more than Python and Bisection.
I've been using sympy to automate some derivative computation and I loved it. I found the evalf[1] function to be really helpful as it helps bridge the gap between the symbolic and numeric world.
I was fairly impressed by this article until I noticed that the author was using Python 2.6. Isn't Python 2.7 the version usually used by those who aren't ready to make the leap to Python 3?
Then I noticed that there was no mention at all of the various ways of using the R language with Python (RPy, etc.).
The easy way to solve equations in Python (or any other language for that matter) would be to take full advantage of the work that other bright people have done.
I think the article would better be described as "Basic methods for numerically solving equations". It gives a nice introduction to a number of those methods, and is really not Python specific as the code can easily be ported to any other language. I would view it as a way to learn the methods, rather than how to actually solve an equation in a program.
If it was more focused on math I would agree. However, the focus is on the solution rather the technique. It is a nice refresher to those who have seen the information before.
The easy way to solve the equation is to send it to wolfram alpha and wait for a result. I think the intention was to not use anything beside standard library. Using the R language is equivalent to using NumPy, SciPy or any other tool. As for the Python version, I have no clue why he uses 2.6. I found that all the code runs perfectly well on 2.7. That is something you should point out to the author.
You missed the point of the article. It's not what version of Python he's using or what libraries already exist. It's about creating an understanding of simple methods for numerically solving equations.
> Isn't Python 2.7 the version usually used
Yes. But because the 2.x series is backward compatible, you will be able to run 2.6 code on 2.7. Honestly, if he was writing on 2.7 the code probably would look identical. For that matter, it would probably look identical on Python 3, except for print statements.
> using the R language
I'm not familiar with R, I always thought it was more for statistics like curve fitting or ANOVA. (I actually am not really sure what ANOVA is, but stats people seem to like that word, so I'll toss it out there.) Sage [1] would be the tool I'd turn to for symbolic equation solving.
I think anything above 2.6 is fair game. Then again, it depends on the system. I have been forced before to develop in python 2.4 and less for specific platforms.
I wrote this a while ago during uni to help cement the various numerical analysis processes covered in my classes in my head. It's not formal or well-written or covers all of the detail perfectly because it was largely for my own use. However, some of you may find it interesting I guess:
Except for the bisection method, all of these implementations take an argument specifying the number of iterations to run. In most cases, the only way to terminate in fewer iterations is by hitting an "exact" root, i.e., calculating the residual as exactly zero. This is poor practice for a number of reasons. First, in practice it's pretty rare for a method to find an exact zero. Second, once a method has converged to the numerical precision of the machine, making more iterations just wastes flops. So a much better approach is to specify a solution tolerance (as shown with the bisection method). Even better is to provide absolute and relative tolerances, and to choose those values based on either the domain requirements, or on the machine characteristics. Dennis & Schnabel's excellent "Numerical Methods for Unconstrained Optimization and Nonlinear Equations" has a good discussion on choosing convergence tolerances.
This dependence on iteration counts to terminate, by the way, is probably why the author equates low iteration counts with greater accuracy. But in fact these methods don't vary in their intrinsic accuracy, rather, they vary in their order of convergence.
Another example of poor practice is in the bisection method implementation. One generally should not bisect an interval using c = (a+b)/2, because the nature of finite-precision arithmetic means there is no guarantee that c will lie between a and b, even if the machine can represent numbers between a and b. A better approach is to ensure a < b, then to set c = a + (b-a)/2. This expression is much less subject to roundoff errors.