Menu
▾
▴
#38 num_solve does not work for simple 2 equations
load_package numeric;
num_solve({3=c*a^4,3=c*a^3},{c,a});
***** Newton method does not converge
as far as I tested, it works for 3=c*a^3, 3=c*a^2.
I also try higher degrees like 3=c*a^5, 3=c*a^4, in this case it is even worse: the method converged but gave wrong answer.
You can try "on trroot,trnumeric;"; to get some fragments of trace output. I believe that num_solve is JUST trying to find a solution numerically, so it can not see the algebraic structure that we as humans can see and that makes us believe that these cases may be simple! I have not yet looked hard at the code that is being used, but I do know that both Newton Raphson is only certain to behave well if somehow it gets a reasonably good initial guess for a solution, and also that in cases where Newton Raphson has not got a close initial starting point deciding how many steps to try to take in the hope that it will converge is hard. I hope that when we look harder at this particular case (ie writing out the Jacobian and looking at the pattern of the iteration in detail) we will find some hack that improves things, but even if we fix this one I will always expect that there are equations that do not yield reliable numerical solutions!
Arthur
OK, A brief further look. num_solve, as you see in the mainual, can be called as in
num_solve()3-c*a^4, 3-c*a^5}, {c=2, a=2.5});
or some such to provide Newton Raphson with starting approximations (I have passed 2 and 1.5 here. If you do not provide starting values it uses a random starting "guess" and a random value is hardly even a guess. In this case if I go "in trnumeric;" to get a bit of trace output and run the same test several times in succession (within the same reduce) I can see the effects of having different random initial "guesses" for c and a, and they generally lead to big residuals and no comfortable convergence. If I provide {2,1.5} it still takes a few steps to home in but it then converges well.
The documentation also explains that the algorithms used in the numeric package are not suited to ill conditioned problems, and that you can give anther parameter to control the number of iterations used. I think I would like to update the documentation to make it clear that unless reasonably accurate starting estimates are provided there can be no certainty of respectable behaviour here, and that multi-variable cases are liable to be harder than things in just one variable. I suspect strongly that error reporting could be improved, but generally finding numerical solutions to arbitrary sets of nonlinear simultaneous equations is nasty so there is not liable to be a quick generic fix,,,,
Arthur
I read the chapters in the book "Numerical Recipe" for C++ and understood that Newton-Raphson method occasionally can't converge; there is a quick remedy, i.e. "tracking backward" that may help a bit.
Also, have you tried {3=c*a^5, 3=c*a^4} ? The problem is that it didn't return with non-convergence, but actually returned some result which was wrong... would there be a final "test" that the result is really a solution? otherwise it may raise some confusion.
Thanks for checking some more. If you grab the full set of Reduce sources from sourceforge you can find the c ode in trunk/packages/numeric/numsolve.red and newton.red. I am about to go away on a summer break so can certainly not look at it for a while, but even though the code is written in Reduce at5 a level you may not be very used to you may find that adjusting the algorithmic details is in fact not too bad. You can not just use code from Numeric Recipes because of their license terms, and stuff in Reduce wants to be under the BSD license. But you can pick up ideas from their explanation (and from professional books on Numerical Analysis). If you could produce a version that was much better behaved we might adopt it. Collecting a set of further examples that are hard (and easy too?) to go in an additional test file or at the end of the current numeric.tst script would also help.
I will see if you have got anywhere when I get back from my break!
I do not know how feasible you will find this, but the file newton.red is only 240 lines long so to my mind can not be TOO bad.
Arthur
Bug was fixed a while ago, but not recorded as such.
I added a regression test file for it.