Testing Zero Finders
Use the historic cubic polynomial $x^3 - 2x - 5$ to test a few zero-finding algorithms.
This takes a long time. It's just linear convergence. But we squeeze through de Morgan's wicket anyway.
Contents
The Footnote
My previous post was about a footnote in a classic numerical analysis text by Whittaker and Robinson that was a quip by Augustus de Morgan. "The reason I call $x^3-2x-5=0$ a celebrated equation is because it was the one on which Wallis chanced to exhibit Newton's method when he first published it, in consequence of which every numerical solver has felt bound in duty to make it one of his examples. Invent a numerical method, neglect to show how it works on this equation, and you are a pilgrim who does not come in at the little wicket (vide J. Bunyan)." So, let's try a few numerical solvers on this equation.Newton's method
I must first try Newton's method. $$ x_{n+1} = x_n - \frac{f(x_n)} {f'(x_n)} $$ I know the derivative and I have a good starting guess.f = @(x) x.^3-2*x-5; fprime = @(x) 3*x.^2-2; x = 2; disp(x) z = 0; while abs(x-z) > eps(x) z = x; x = x - f(x)/fprime(x); disp(x) end
2 2.100000000000000 2.094568121104185 2.094551481698199 2.094551481542327 2.094551481542327This reaches the solution very quickly. Only five steps, with the fifth confirming that the fourth was already correct to full precision. But Newton's method requires knowledge of the derivative and a really good starting point. These are big drawbacks for general use.
fzero
Let's see how the zero finder available in MATLAB works. I wrote a series of postings about it not too long ago; part1, part2, part3. Set a display parameter asking for intermediate results.% opt = optimset('display','iter');
Pretend we don't know anything about this function and start the search at $x = 0$.
% fzero(f,0,opt)
We get one line of output each time a is decreased and b is increased. All the values of f(a) are negative and all the values of f(b) are also negative, until b = 2.56 when the first sign change is reached.
% Search for an interval around 0 containing a sign change: % Func-count a f(a) b f(b) Procedure % 1 0 -5 0 -5 initial interval % 3 -0.0282843 -4.94345 0.0282843 -5.05655 search % 5 -0.04 -4.92006 0.04 -5.07994 search % 7 -0.0565685 -4.88704 0.0565685 -5.11296 search % 9 -0.08 -4.84051 0.08 -5.15949 search % 11 -0.113137 -4.77517 0.113137 -5.22483 search % 13 -0.16 -4.6841 0.16 -5.3159 search % 15 -0.226274 -4.55904 0.226274 -5.44096 search % 17 -0.32 -4.39277 0.32 -5.60723 search % 19 -0.452548 -4.18759 0.452548 -5.81241 search % 21 -0.64 -3.98214 0.64 -6.01786 search % 23 -0.905097 -3.93126 0.905097 -6.06874 search % 25 -1.28 -4.53715 1.28 -5.46285 search % 27 -1.81019 -7.31125 1.81019 -2.68875 search % 29 -2.56 -16.6572 2.56 6.65722 searchNow the classic zeroin algorithm can quickly and rapidly find the zero. To see the details, run fzerogui from Numerical Computing with MATLAB and click on the auto button. Both the secant method and inverse quadratic interpolation are used to find the zero.
% start -2.5600000000000001 % start 2.5600000000000001 % secant 1.0980323260716793 % secant 1.7832168816106038 % iqi 2.2478393639958036 % secant 2.0660057758331045 % secant 2.0922079131171945 % iqi 2.0945566700001779 % secant 2.0945514746903111 % secant 2.0945514815423065 % iqi 2.0945514815423265 % small 2.0945514815423274So, even when it is started at a lousy initial guess, fzero passes de Morgan's test.
Fixed point iteration
I like to call this the "WS" algorithm for "World's Simplest". Try to find a fixed point of a mapping $G(x)$. $$ x = G(x) $$ The iteration is $$ x_{n+1} = G(x_n) $$ What should I choose for $G$ ? The most obvious choice is the equation $$ x = \frac{1}{2} (x^3-5) $$ But near the zero the slope of $G(x)$ is too large and the WS iteration diverges. Another choice is $$ x = \sqrt[3]{2x+5} $$ This will work. While we're at it, also produce a plot.G = @(x) (2*x+5).^(1/3) ezplot(G,[1 3]) line([1 3],[1 3],'color','k') dkgreen = [0 2/3 0]; x = 1.4; z = 0; disp(x) while abs(x-z) > eps(x) z = x; x = G(x); line([z z],[z x],'color',dkgreen) line([z x],[x x],'color',dkgreen) disp(x) end
G =
@(x)(2*x+5).^(1/3)
1.400000000000000
1.983192482680775
2.077490885128178
2.091955753470501
2.094156962781544
2.094491529117444
2.094542371187228
2.094550097140211
2.094551271169830
2.094551449574315
2.094551476684497
2.094551480804135
2.094551481430152
2.094551481525281
2.094551481539736
2.094551481541933
2.094551481542267
2.094551481542317
2.094551481542325
2.094551481542326
2.094551481542327
This takes a long time. It's just linear convergence. But we squeeze through de Morgan's wicket anyway.
Root squaring
The "Graffe" root-squaring method was invented independently by Germinal Pierre Dandelin in 1826, Nikolai Lobachevsky in 1834, and Karl Heinrich Graffe in 1837. An article by Alston Householder referenced below goes into detail about who invented what. The idea is to transform the coefficients of one polynomial into another whose zeros are the squares of the original. If the zeros are well separated in magnitude, repeating the process will eventually produce a high degree polynomial whose first two terms provide a good approximation to the dominant zero. The method was popular for hand computation, but there are serious difficulties when it comes to making a reliable automatic implementation. Repeated zeros, complex zeros of equal magnitude, and especially scaling to prevent floating point overflow and underflow are all obstacles. But it works OK on our historic cubic, which has a simple dominant zero. Suppose $p(x)$ is a cubic polynomial. $$ p(x) = p_1 x^3 + p_2 x^2 + p_3 x + p_4 $$ Rearrange the equation $p(x) = 0$ so that the terms with odd powers of $x$ are on one side of the equals sign and the even powers are on the other. $$ p_1 x^3 + p_3 x = - (p_2 x^2 + p_4) $$ Square both sides of this equation and collect terms to form a new polynomial, $q(x)$. $$ q(x) = (p_1 x^3 + p_3 x)^2 - (p_2 x^2 + p_4)^2 $$ $$ = p_1^2 x^6 - (p_2^2 - 2 p_1 p_3) x^4 + (p_3^2 - 2 p_2 p_4) x^2 - p_4^2 $$ The zeros of $q$ are the squares of the zeros of $p$. Here is a function that does the job for cubics.function q = graffe(p); % Graffe root squaring. q = graffe(p). % p & q are 4-vectors with coefficients of cubics. % roots(q) = roots(p).^2. q = zeros(size(p)); q(1) = p(1)^2; q(2) = -(p(2)^2 - 2*p(1)*p(3)); q(3) = p(3)^2 - 2*p(2)*p(4); q(4) = -p(4)^2; endLet's try it on our historic polynomial $x^3-2x-5$. The floating point exponents of the coefficients soon start doubling at each step. Luckily, we reach my stopping criterion in eight steps. If we had to try more steps without rescaling, we would overflow. But by this time, the 256-th power of the zero we're after completely dominates the powers of the other roots, and so the 256-th root of the ratio of the just first coefficients provides a result that is accurate to full precision.
p = [1 0 -2 -5]; fmt = '%5.0f %5.0g %13.5g %13.5g %13.5g %20.15f\n'; fprintf(fmt,1,p,0) m = 1; while max(abs(p)) < sqrt(realmax) m = 2*m; p = graffe(p); z = (-p(2)/p(1)).^(1/m); fprintf(fmt,m,p,z); end
1 1 0 -2 -5 0.000000000000000
2 1 -4 4 -25 2.000000000000000
4 1 -8 -184 -625 1.681792830507429
8 1 -432 23856 -3.9063e+05 2.135184796196703
16 1 -1.3891e+05 2.3161e+08 -1.5259e+11 2.096144583759898
32 1 -1.8833e+10 1.125e+16 -2.3283e+22 2.094553557477412
64 1 -3.5467e+20 -7.5043e+32 -5.421e+44 2.094551481347086
128 1 -1.2579e+41 1.7861e+65 -2.9387e+89 2.094551481542327
256 1 -1.5824e+82 -4.2032e+130 -8.6362e+178 2.094551481542327
This basic implementation of root-squaring makes it through De Morgan's wicket.
Reduced Penultimate Remainder
Almost two years ago, I blogged about the Reduced Penultimate Remainder algorithm. It was an undergrad individual research project that I now fondly remember as the most obscure, impractical algorithm I have even seen. The algorithm tries to find a lower degree polynomial that divides exactly, without any remainder, into the polynomial whose zeros we seek. The zeros of the resulting divisor are then some of the zeros of the dividend. As an example in that post, I tried to find a linear factor of $x^3-2x-5$. That would have given me the zero that I've been computing in this post. But the algorithm fails to converge, no matter what starting factor is provided. How about computing a quadratic factor instead? If successful, that will yield two zeros. Here is the function that takes one step.function q = rpr(p,q) % RPR Reduced Penultimate Remainder. % r = rpr(p,q), for polynomials p and q. % Assume p(1) = 1, q(1) = 1. % Polynomial long division. % Stop when degree of r = degree of q. n = length(p); m = length(q); r = p(1:m); for k = m+1:n r = r - r(1)*q; r = [r(2:end) p(k)]; end q = r/r(1); endThe starting factor $q(x)$ can be any quadratic that has imaginary zeros. Let's try $q(x) = x^2+x+1$.
% p = [1 0 -2 -5]; % q = [1 1 1] % k = 1; % r = 0*q; % while max(abs(q-r)) > 10*eps(max(abs(q))) % r = q; % q = rpr(p,q) % k = k+1; % end % kHere are the first ten and the last four of 114 steps. I've edited out 100 steps in between.
1.000000000000000 1.000000000000000 1.000000000000000
1.000000000000000 3.000000000000000 5.000000000000000
1.000000000000000 2.333333333333334 1.666666666666667
1.000000000000000 1.571428571428571 2.142857142857143
1.000000000000000 2.636363636363636 3.181818181818182
1.000000000000000 1.965517241379310 1.896551724137931
1.000000000000000 1.982456140350877 2.543859649122807
1.000000000000000 2.292035398230088 2.522123893805309
1.000000000000000 1.972972972972974 2.181467181467182
1.000000000000000 2.119373776908023 2.534246575342465
. . .
. . .
. . .
1.000000000000000 2.094551481542332 2.387145908831160
1.000000000000000 2.094551481542323 2.387145908831149
1.000000000000000 2.094551481542327 2.387145908831159
1.000000000000000 2.094551481542328 2.387145908831155
k =
114
I can now use deconvolution to do the polynomial division $p(x)/q(x)$. This produces a linear quotient and the desired zero.
% r = deconv(p,q) % x1 = -r(1)r = 1.000000000000000 -2.094551481542328 x1 = 2.094551481542328
评论
要发表评论，请点击 此处 登录到您的 MathWorks 帐户或创建一个新帐户。