[math-fun] Cubic eqn with 3 real zeros, again
Has anyone seen the following treatment of the real cubic? A cubic equation with 3 real zeros can be reduced to either (your choice): 4*c^3-3*c=c3 or 4*s^3-3*s=-s3 where c,c3 reminds us of cos(a),cos(3*a), and s,s3 reminds us of sin(b),sin(3*b), respectively. Now Newton's method _can_ be used to compute one of the zeros, but an unfortunate choice of starting point will attempt to converge on a double (or almost double) root, thereby seriously slowing the convergence. For example, Newton on the cosine version gives c' = (c3+8*c^3)/[3*(2*c-1)*(2*c+1)] which leads to very unfortunate behavior for c near +- 1/2. But if we recognize that we are attempting to solve the 1-dimensional projection of a 2-dimensional problem, then we can avoid this convergence problem. Specifically, if we are given 4*c^3-3*c=c3, then we first compute s3=sqrt(1-c3^2), and then instead solve the problem z^3 = z3 = = c3+i*s3 = (c+i*s)^3 = c3+i*s3 = c^3 - 3*c*s^2 - i*(s^3 - 3*s*c^2) = c^3 - 3*c*(1-c^2) - i*(s^3 - 3*s*(1-s^2)) = c^3 - 3*c +3*c^3 - i*(s^3 - 3*s + 3*s^3) = 4*c^3 - 3*c - i*(4*s^3 - 3*s) = c3 - i*(-s3) = c3 + i*s3 But Newton's method on the complex problem z^3=z3 is much more robust, so long as we stay away from z=0: z' = (z3+2*z^3)/(3*z^2) (Throw away s=imagpart(z) & keep c=realpart(z) for the answer. But before doing so, we can trivially produce the other two roots with complex multiplications: w*z and w^2*z, where w is cube root of unity: w^3=1.) But we know that |z|^2 must be 1, so we can easily avoid the region near z=0. On modern hardware with multiple floating point units, complex arithmetic may have no more latency than real arithmetic, so this method may be no slower, and (due to faster convergence) may actually be faster.
This suggests a meta-problem: If you find your numerical algorithm is converging slowly, can you introduce an extra dimension, perhaps imaginary, to speed things up? Rich ------------- Quoting Henry Baker <hbaker1@pipeline.com>:
Has anyone seen the following treatment of the real cubic?
A cubic equation with 3 real zeros can be reduced to either (your choice):
4*c^3-3*c=c3 or 4*s^3-3*s=-s3
where c,c3 reminds us of cos(a),cos(3*a), and s,s3 reminds us of sin(b),sin(3*b), respectively.
Now Newton's method _can_ be used to compute one of the zeros, but an unfortunate choice of starting point will attempt to converge on a double (or almost double) root, thereby seriously slowing the convergence.
For example, Newton on the cosine version gives
c' = (c3+8*c^3)/[3*(2*c-1)*(2*c+1)]
which leads to very unfortunate behavior for c near +- 1/2.
But if we recognize that we are attempting to solve the 1-dimensional projection of a 2-dimensional problem, then we can avoid this convergence problem.
Specifically, if we are given 4*c^3-3*c=c3, then we first compute s3=sqrt(1-c3^2), and then instead solve the problem
z^3 = z3 = = c3+i*s3 = (c+i*s)^3 = c3+i*s3 = c^3 - 3*c*s^2 - i*(s^3 - 3*s*c^2) = c^3 - 3*c*(1-c^2) - i*(s^3 - 3*s*(1-s^2)) = c^3 - 3*c +3*c^3 - i*(s^3 - 3*s + 3*s^3) = 4*c^3 - 3*c - i*(4*s^3 - 3*s) = c3 - i*(-s3) = c3 + i*s3
But Newton's method on the complex problem z^3=z3 is much more robust, so long as we stay away from z=0:
z' = (z3+2*z^3)/(3*z^2)
(Throw away s=imagpart(z) & keep c=realpart(z) for the answer. But before doing so, we can trivially produce the other two roots with complex multiplications: w*z and w^2*z, where w is cube root of unity: w^3=1.)
But we know that |z|^2 must be 1, so we can easily avoid the region near z=0.
On modern hardware with multiple floating point units, complex arithmetic may have no more latency than real arithmetic, so this method may be no slower, and (due to faster convergence) may actually be faster.
_______________________________________________ math-fun mailing list math-fun@mailman.xmission.com http://mailman.xmission.com/cgi-bin/mailman/listinfo/math-fun
Henry's post got me thinking in a different direction, maybe not the one he had intended. I was thinking about Newton's method and the problem that it is vulnerable to instabilities; even when it works perfectly and converges fast, it gets preoccupied with one root and ignores the others unless you start it over. So I started wondering whether anybody had played with a sort of adaptive version of Newton's algorithm. While one "thread" is converging on a root, a "kibitzer" thread divides the original function by x-r, as if r were the actual root, and gets a modified function that can be the host for another round of Newton's algorithm that would be scared away from the first root and be more likely to converge on the others. Perhaps I'm barking up a squirrel-less tree, but I was imagining different threads converging on different roots, and dividing out by each other's approximations in order to "isolate" their own target root. Has anyone heard of an approach like this? On Wed, Sep 15, 2010 at 6:14 PM, <rcs@xmission.com> wrote:
This suggests a meta-problem: If you find your numerical algorithm is converging slowly, can you introduce an extra dimension, perhaps imaginary, to speed things up?
Rich
-------------
Quoting Henry Baker <hbaker1@pipeline.com>:
Has anyone seen the following treatment of the real cubic?
A cubic equation with 3 real zeros can be reduced to either (your choice):
4*c^3-3*c=c3 or 4*s^3-3*s=-s3
where c,c3 reminds us of cos(a),cos(3*a), and s,s3 reminds us of sin(b),sin(3*b), respectively.
Now Newton's method _can_ be used to compute one of the zeros, but an unfortunate choice of starting point will attempt to converge on a double (or almost double) root, thereby seriously slowing the convergence.
For example, Newton on the cosine version gives
c' = (c3+8*c^3)/[3*(2*c-1)*(2*c+1)]
which leads to very unfortunate behavior for c near +- 1/2.
But if we recognize that we are attempting to solve the 1-dimensional projection of a 2-dimensional problem, then we can avoid this convergence problem.
Specifically, if we are given 4*c^3-3*c=c3, then we first compute s3=sqrt(1-c3^2), and then instead solve the problem
z^3 = z3 = = c3+i*s3 = (c+i*s)^3 = c3+i*s3 = c^3 - 3*c*s^2 - i*(s^3 - 3*s*c^2) = c^3 - 3*c*(1-c^2) - i*(s^3 - 3*s*(1-s^2)) = c^3 - 3*c +3*c^3 - i*(s^3 - 3*s + 3*s^3) = 4*c^3 - 3*c - i*(4*s^3 - 3*s) = c3 - i*(-s3) = c3 + i*s3
But Newton's method on the complex problem z^3=z3 is much more robust, so long as we stay away from z=0:
z' = (z3+2*z^3)/(3*z^2)
(Throw away s=imagpart(z) & keep c=realpart(z) for the answer. But before doing so, we can trivially produce the other two roots with complex multiplications: w*z and w^2*z, where w is cube root of unity: w^3=1.)
But we know that |z|^2 must be 1, so we can easily avoid the region near z=0.
On modern hardware with multiple floating point units, complex arithmetic may have no more latency than real arithmetic, so this method may be no slower, and (due to faster convergence) may actually be faster.
_______________________________________________ math-fun mailing list math-fun@mailman.xmission.com http://mailman.xmission.com/cgi-bin/mailman/listinfo/math-fun
_______________________________________________ math-fun mailing list math-fun@mailman.xmission.com http://mailman.xmission.com/cgi-bin/mailman/listinfo/math-fun
* Allan Wechsler <acwacw@gmail.com> [Sep 16. 2010 08:09]:
Henry's post got me thinking in a different direction, maybe not the one he had intended. I was thinking about Newton's method and the problem that it is vulnerable to instabilities;
due to the nature of the basins of attraction as soon as we got at least three roots. Some people would refer to Newton's method as "root polishing" rather than "root finding". (Yes there are methods to mitigate the problem like that damping factor multiplied to the amount of change, but I assume one can always cook up an example killing any particular root finder).
even when it works perfectly and converges fast, it gets preoccupied with one root and ignores the others unless you start it over.
So I started wondering whether anybody had played with a sort of adaptive version of Newton's algorithm. While one "thread" is converging on a root, a "kibitzer" thread divides the original function by x-r, as if r were the actual root, and gets a modified function that can be the host for another round of Newton's algorithm that would be scared away from the first root and be more likely to converge on the others. Perhaps I'm barking up a squirrel-less tree, but I was imagining different threads converging on different roots, and dividing out by each other's approximations in order to "isolate" their own target root.
Has anyone heard of an approach like this?
This may not be what you have in mind but there are methods to find _all_ roots simultaneously, see e.g. http://en.wikipedia.org/wiki/Durand-Kerner_method and the references at the bottom of the page
On Wed, Sep 15, 2010 at 6:14 PM, <rcs@xmission.com> wrote:
This suggests a meta-problem: If you find your numerical algorithm is converging slowly, can you introduce an extra dimension, perhaps imaginary, to speed things up?
Rich
Ignoring the "extra dimension": there is a method (essentially delta-squared) to obtain a superlinear iteration from a linear one. See page 165 (Thm.4.3.3) of Alston S.\ Householder: {The Numerical Treatment of a Single Nonlinear Equation}, McGraw-Hill, (1970) In simple terms, if \Phi(x) is an (only linear) iteration ( i.e. you iterate x_{n+1} = \Phi(x_{n}) ), then \Phi^* defined by 2 * [\Phi(\Phi(x)) - \Phi(x)] \Phi (x) = \Phi(\Phi(x)) - ----------------------------- \Phi(\Phi(x)) - 2 \Phi(x) + x will have superlinear convergence. (Properly typeset as relation (30.6-7) on page 598 of the fxtbook, followed by a simple example).
[...]
cheers, jj
On 9/16/10, Joerg Arndt <arndt@jjj.de> wrote:
... This may not be what you have in mind but there are methods to find _all_ roots simultaneously, see e.g. http://en.wikipedia.org/wiki/Durand-Kerner_method and the references at the bottom of the page
A very elegant method which computes all roots simultaneously, yielding circles in which they are guaranteed to lie, is buried (I don't remember exactly where --- anybody else know?) in the strangely neglected Henrici P., Applied and Computational Complex Analysis (Wiley). [Three volumes: 1974, 1977, 1986.] It's interesting that these methods actually gain in speed and stability by simultaneity, as opposed to attempting to exclude the other roots. WFL
How is it related to Schonhage's splitting circle method? http://en.wikipedia.org/wiki/Splitting_circle_method Victor On Thu, Sep 16, 2010 at 1:29 PM, Fred lunnon <fred.lunnon@gmail.com> wrote:
On 9/16/10, Joerg Arndt <arndt@jjj.de> wrote:
... This may not be what you have in mind but there are methods to find _all_ roots simultaneously, see e.g. http://en.wikipedia.org/wiki/Durand-Kerner_method and the references at the bottom of the page
A very elegant method which computes all roots simultaneously, yielding circles in which they are guaranteed to lie, is buried (I don't remember exactly where --- anybody else know?) in the strangely neglected
Henrici P., Applied and Computational Complex Analysis (Wiley). [Three volumes: 1974, 1977, 1986.]
It's interesting that these methods actually gain in speed and stability by simultaneity, as opposed to attempting to exclude the other roots.
WFL
_______________________________________________ math-fun mailing list math-fun@mailman.xmission.com http://mailman.xmission.com/cgi-bin/mailman/listinfo/math-fun
That's one sophisticated root-finding algorithm! Henry has evidently researched Henrici's algorithm more thoroughly than I --- but my recollection is that there's no connection. WFL On 9/16/10, Victor Miller <victorsmiller@gmail.com> wrote:
How is it related to Schonhage's splitting circle method? http://en.wikipedia.org/wiki/Splitting_circle_method
Victor
On Thu, Sep 16, 2010 at 1:29 PM, Fred lunnon <fred.lunnon@gmail.com> wrote:
On 9/16/10, Joerg Arndt <arndt@jjj.de> wrote:
... This may not be what you have in mind but there are methods to find _all_ roots simultaneously, see e.g. http://en.wikipedia.org/wiki/Durand-Kerner_method and the references at the bottom of the page
A very elegant method which computes all roots simultaneously, yielding circles in which they are guaranteed to lie, is buried (I don't remember exactly where --- anybody else know?) in the strangely neglected
Henrici P., Applied and Computational Complex Analysis (Wiley). [Three volumes: 1974, 1977, 1986.]
It's interesting that these methods actually gain in speed and stability by simultaneity, as opposed to attempting to exclude the other roots.
WFL
_______________________________________________ math-fun mailing list math-fun@mailman.xmission.com http://mailman.xmission.com/cgi-bin/mailman/listinfo/math-fun
_______________________________________________ math-fun mailing list math-fun@mailman.xmission.com http://mailman.xmission.com/cgi-bin/mailman/listinfo/math-fun
Considering those huge lists of surprises, I suspect that the really greatest surprise is that mathematics works at all. However, while following the discussion of cubic equations, it occurred to me that the coefficients of a (monic) polynomial just come from the hierarchy of surface areas of the rectangular figure formed by the roots (due care being given to signs). So the constant term is a volume, the coefficient of x the sum of surface areas, the coefficient of x^2 the sum of the edges, and so on. Probably everybody knew this; just never said so. -hvm ------------------------------------------------- www.correo.unam.mx UNAMonos Comunicándonos
I haven't been following this thread, but I want to mention that Peter Doyle and Curt McMullen made a global analysis of all possible root-finding algorithms for polynomials that involve iterating an complex analytic function of one variable with parameters depending on the coefficients of the polynomial. (Newton's method is one such iteration, but it gets trapped in unproductive loops with positive probability already for cubics. They showed that algorithms can be constructed that converge a.s. for polynomials of degrees 1 through 5, but no such method converges a.s. to a root for polynomials of degree 6 or higher. The proof is based on the fact that the simple parts of Galois groups through degree 5 are all automorphism groups of polyhedra, but this is false for higher degrees. The construction for quintics is ingenious, making essential use of the theory Klein developed describing invariants for the icosahedral group, expressing them in terms of the rational functions that describes the map of CP^1 to its symmetry quotient by the icosahedral group. I think these methods are all super-exponential once they're in the attracting basin, but I don't remember Bill On Sep 17, 2010, at 7:41 AM, Fred lunnon wrote:
That's one sophisticated root-finding algorithm!
Henry has evidently researched Henrici's algorithm more thoroughly than I --- but my recollection is that there's no connection. WFL
On 9/16/10, Victor Miller <victorsmiller@gmail.com> wrote:
How is it related to Schonhage's splitting circle method? http://en.wikipedia.org/wiki/Splitting_circle_method
Victor
On Thu, Sep 16, 2010 at 1:29 PM, Fred lunnon <fred.lunnon@gmail.com> wrote:
On 9/16/10, Joerg Arndt <arndt@jjj.de> wrote:
... This may not be what you have in mind but there are methods to find _all_ roots simultaneously, see e.g. http://en.wikipedia.org/wiki/Durand-Kerner_method and the references at the bottom of the page
A very elegant method which computes all roots simultaneously, yielding circles in which they are guaranteed to lie, is buried (I don't remember exactly where --- anybody else know?) in the strangely neglected
Henrici P., Applied and Computational Complex Analysis (Wiley). [Three volumes: 1974, 1977, 1986.]
It's interesting that these methods actually gain in speed and stability by simultaneity, as opposed to attempting to exclude the other roots.
WFL
_______________________________________________ math-fun mailing list math-fun@mailman.xmission.com http://mailman.xmission.com/cgi-bin/mailman/listinfo/math-fun
_______________________________________________ math-fun mailing list math-fun@mailman.xmission.com http://mailman.xmission.com/cgi-bin/mailman/listinfo/math-fun
_______________________________________________ math-fun mailing list math-fun@mailman.xmission.com http://mailman.xmission.com/cgi-bin/mailman/listinfo/math-fun
Thanks for the reference! Here's the Doyle-McMullen paper: http://math.dartmouth.edu/~doyle/docs/icos/icos.pdf http://www.math.dartmouth.edu/~doyle/docs/icos/icos/icos.html At 04:12 AM 9/20/2010, Bill Thurston wrote:
I haven't been following this thread, but I want to mention that Peter Doyle and Curt McMullen made a global analysis of all possible root-finding algorithms for polynomials that involve iterating an complex analytic function of one variable with parameters depending on the coefficients of the polynomial. (Newton's method is one such iteration, but it gets trapped in unproductive loops with positive probability already for cubics. They showed that algorithms can be constructed that converge a.s. for polynomials of degrees 1 through 5, but no such method converges a.s. to a root for polynomials of degree 6 or higher.
The proof is based on the fact that the simple parts of Galois groups through degree 5 are all automorphism groups of polyhedra, but this is false for higher degrees. The construction for quintics is ingenious, making essential use of the theory Klein developed describing invariants for the icosahedral group, expressing them in terms of the rational functions that describes the map of CP^1 to its symmetry quotient by the icosahedral group.
I think these methods are all super-exponential once they're in the attracting basin, but I don't remember
Bill
* Bill Thurston <wpt4@cornell.edu> [Sep 20. 2010 14:46]:
I haven't been following this thread, but I want to mention that Peter Doyle and Curt McMullen made a global analysis of [...]
This one(?): Peter Doyle, Curt McMullen: {Solving the quintic by iteration}, Acta Mathematica, vol.163, no.1, pp.151–180, (1989). http://www.math.harvard.edu/~ctm/papers/home/text/papers/icos/icos.pdf (Sadly much over my head). This is ref'd form iteration recipe at end of: http://en.wikipedia.org/wiki/Bring_radical
I think these methods are all super-exponential once they're in the attracting basin, but I don't remember
Bill
(Assuming the answer is 'yes', which I think it is): Still if you are very close to a boundary the method will _seem_ to fail. jj Possibly related: Scott Crass: "Solving the quintic by iteration in three dimensions" http://arxiv.org/abs/math/9903054
Well, moving into the complex plane helps in this instance because the roots become more widely separated in the complex plane. Some of the same intuition that we already have about series converging to the nearest singularity also works here. I've been musing about whether there is an analogous extension to the quaternions for solving quartics??? At 03:14 PM 9/15/2010, rcs@xmission.com wrote:
This suggests a meta-problem: If you find your numerical algorithm is converging slowly, can you introduce an extra dimension, perhaps imaginary, to speed things up?
Rich
------------- Quoting Henry Baker <hbaker1@pipeline.com>:
Has anyone seen the following treatment of the real cubic?
A cubic equation with 3 real zeros can be reduced to either (your choice):
4*c^3-3*c=c3 or 4*s^3-3*s=-s3
where c,c3 reminds us of cos(a),cos(3*a), and s,s3 reminds us of sin(b),sin(3*b), respectively.
Now Newton's method _can_ be used to compute one of the zeros, but an unfortunate choice of starting point will attempt to converge on a double (or almost double) root, thereby seriously slowing the convergence.
For example, Newton on the cosine version gives
c' = (c3+8*c^3)/[3*(2*c-1)*(2*c+1)]
which leads to very unfortunate behavior for c near +- 1/2.
But if we recognize that we are attempting to solve the 1-dimensional projection of a 2-dimensional problem, then we can avoid this convergence problem.
Specifically, if we are given 4*c^3-3*c=c3, then we first compute s3=sqrt(1-c3^2), and then instead solve the problem
z^3 = z3 = = c3+i*s3 = (c+i*s)^3 = c3+i*s3 = c^3 - 3*c*s^2 - i*(s^3 - 3*s*c^2) = c^3 - 3*c*(1-c^2) - i*(s^3 - 3*s*(1-s^2)) = c^3 - 3*c +3*c^3 - i*(s^3 - 3*s + 3*s^3) = 4*c^3 - 3*c - i*(4*s^3 - 3*s) = c3 - i*(-s3) = c3 + i*s3
But Newton's method on the complex problem z^3=z3 is much more robust, so long as we stay away from z=0:
z' = (z3+2*z^3)/(3*z^2)
(Throw away s=imagpart(z) & keep c=realpart(z) for the answer. But before doing so, we can trivially produce the other two roots with complex multiplications: w*z and w^2*z, where w is cube root of unity: w^3=1.)
But we know that |z|^2 must be 1, so we can easily avoid the region near z=0.
On modern hardware with multiple floating point units, complex arithmetic may have no more latency than real arithmetic, so this method may be no slower, and (due to faster convergence) may actually be faster.
participants (8)
-
Allan Wechsler -
Bill Thurston -
Fred lunnon -
Henry Baker -
Joerg Arndt -
mcintosh＠servidor.unam.mx -
rcs＠xmission.com -
Victor Miller