greg@skippy.berkeley.edu (Greg) (03/17/88)
In article <7741@alice.UUCP> td@alice.UUCP writes: >Restating, this is > p0(x)=p(x) > p1(x)=p'(x) > p[n+2](x)=p[n](x) mod p[n-1](x) >where a(x) mod b(x) is the remainder left after `synthetic division.' > >This is incorrect. Sturm chains have > > p[n+2](x)= -(p[n](x) mod p[n-1](x)) Given that the purpose of computing the chain is to compute S(a)-S(b), where S(x) is the number of sign changes in the chain evaluated at x, your modification has a trivial effect on the algorithm. The new number of sign changes is S'(x) = k - S(x), where k is the length of the chain, and therefore S'(b) - S'(a) = S(a) - S(b). Perhaps the point of your sign convention is that S'(b)-S'(a) resembles F(b)-F(a) in the fundamental theorem of calculus. >>Then the number of roots in the interval [a,b] is S(a)-S(b). >This is wrong for two reasons: first, the interval must be open at >the bottom (i.e. (a,b], not [a,b]). You are correct. However, this bug is mild for a floating-point algorithm. Also, in my case, the best plan of attack if either endpoint is a root is to abort the Sturm chain computation and go for the roots directly. We know at this point there are at most two roots in [a,b]. If p(a) = p(b) = 0, we are done. If p(a) = 0, we can check if there is another root by comparing the sign of p'(a) with that of p(b). p'(a) can't be zero because double roots of p have been factored out. >Second the method of counting >sign changes is subtly wrong. It is possible for p[k](a) to be zero. >Any zeros in the sequence must be deleted before counting sign changes. You have a point, but again the bug is mild. Computationally the fastest solution is to declare zero to be either positive or negative, because the two terms on either side of a zero in a Sturm chain always have opposite sign. >Another criticism of the method is that his numerical method for converging >to the root once it's isolated is over-elaborate. Since the root is bracketed, >methods like regula falsi are guaranteed to converge, albeit slowly. The whole point of the algorithm is that the root is bracketed in such a way that Newton's method will work beautifully. As a general rule, Newton's method either crashes, stumbles in the dark, or converges much more quickly than an uninspired method such as regula falsi. In this case, there are no inflection points to interfere. There may be two roots, which would make regula falsi choke, but Newton steps from the two endpoints will get you both roots. If there is one root, you pick the starting endpoint based on concavity and again Newton's method can't fail. For those of you who haven't seen Newton's method, the formula for it is: x[n+1] = x[n] - f(x[n])/f'(x[n]). If the initial guess x[0] is good, each guess will have twice the precision of the previous as long as f'(x) is continuous and non-zero near the root. -- Greg
td@alice.UUCP (03/17/88)
greg@skippy claims that the sequence p0(x)=p(x) p1(x)=p'(x) p[n+2](x)=p[n](x) mod p[n-1](x) where a(x) mod b(x) is the remainder left after `synthetic division.' is an adequate sequence for use in applying Sturm's theorem. I claimed that it was not, but that p0(x)=p(x) p1(x)=p'(x) p[n+2](x)=-p[n](x) mod p[n-1](x) greg@skippy followed this up, claiming: >Given that the purpose of computing the chain is to compute S(a)-S(b), >where S(x) is the number of sign changes in the chain evaluated at x, >your modification has a trivial effect on the algorithm. The new >number of sign changes is S'(x) = k - S(x), where k is the length of >the chain, and therefore S'(b) - S'(a) = S(a) - S(b). This is untrue for two reasons. First, I do not claim that S'(b)-S'(a) is the correct formula. The correct formula is S'(a)-S'(b). Second, the modification does not have the claimed trivial effect. Rather it corrects an incorrect method. I will give an example: Let p(x)=x^2-x. This has roots 0 and 1. According to greg@skippy, a Sturm sequence for this polynomial and its values at a=-1 and b=2 are: p[i](x) p[i](-1) p[i](2) x^2-x 2 2 2*x-1 -3 3 -1/4 -1/4 -1/4 So, S(a)=1 and S(b)=1, yielding S(a)-S(b)=0. But p(x) has two roots on the interval (-1,2]. According to me, a Sturm sequence for this polynomial and its values at a=-1 and b=2 are: p[i](x) p[i](-1) p[i](2) x^2-x 2 2 2*x-1 -3 3 1/4 1/4 1/4 In this case, S'(a)=2 and S'(b)=0, so S'(a)-S'(b)=2, as required. My previous criticism of greg@skippy's numerical methods was based on not-too-careful reading of his posting. If you are interested in finding all real roots of a polynomial, then his method is endorsable. In a ray-tracing situation, we are interested in finding the smallest positive root, in which case I believe his method to be overkill. I have not implemented it, so that is an opinion.