使用历史上的三次多项式$x^3 - 2x - 5$来测试一些寻零算法。

内容

的脚注<一个name="07636ff6-6cc5-4796-886a-8b1261de4b46">

我的<一个href="//www.tatmou.com/blogs/cleve/2015/12/21/a-historic-cubic/">以前的文章是关于惠特克和罗宾逊的经典数值分析文本中的一个脚注，这是奥古斯都·德·摩根的妙语。“我之所以把$x^3-2x-5=0$称为一个著名的方程，是因为沃利斯第一次发表牛顿的方法时，碰巧在这个方程上展示了牛顿的方法，因此每个数值求解者都觉得有责任把它作为他的例子之一。发明了一种数值方法，却忽视了它是如何在这个方程上工作的，你就是一个不从小门上进来的朝圣者。见j·班扬)。”让我们尝试一些数值解法来解这个方程。

牛顿法<一个name="3e4aef46-7453-44ea-838d-76ee4e780934">

我必须先试试牛顿法。$$ x_{n+1} = x_n - \frac{f(x_n)} {f'(x_n)} $$我知道导数，我有一个很好的开始猜想。

F = @(x) x.^3-2*x-5;Fprime = @(x) 3*x.^2-2;X = 2;Disp (x) z = 0;而Abs (x-z) > eps(x) z = x;X = X - f(X)/fprime(X);disp (x)结束

2 2.100000000000000 2.094568121104185 2.094551481698199 2.094551481542327 2.094551481542327

这很快就得到了解决方案。只有五步，第五步证实了第四步已经精确无误。但是牛顿的方法需要导数的知识和一个很好的起点。这些都是通用的大缺点。

fzero<一个name="afa9423c-2ed5-466f-8eb6-46dbd120af78">

让我们看看MATLAB中的零点查找器是如何工作的。不久前，我写了一系列关于它的帖子;<一个href="//www.tatmou.com/blogs/cleve/2015/10/12/zeroin-part-1-dekkers-algorithm/">part1，<一个href="//www.tatmou.com/blogs/cleve/2015/10/26/zeroin-part-2-brents-version/">第二部分，<一个href="//www.tatmou.com/blogs/cleve/2015/11/09/zeroin-part-3-matlab-zero-finder-fzero/">part3．设置一个请求中间结果的显示参数。

% opt = optimset('display'，'iter');

假设我们对这个函数一无所知，然后从$x = 0$开始搜索。

% fzero (f, 0,选择)

每次我们得到一行输出一个减少和b是增加了。所有的值f (a)的所有值都是负的f (b)也都是负的，直到B = 2.56当到达第一个符号变化时。

搜索0周围包含符号变化的区间:% Func-count a f(a) b f(b)步骤% 1 0 -5 0 -5初始间隔% 3 -0.0282843 -4.94345 0.0282843 -5.05655搜索% 5 -0.04 -4.92006 0.04 -5.07994搜索% 7 -0.0565685 -4.88704 0.0565685 -5.11296搜索% 9 -0.08 -4.84051 0.08 -5.15949搜索% 11 -0.113137 -4.77517 0.113137 -5.22483搜索% 13 -0.16 -4.6841 0.16 -5.3159搜索% 15 -0.226274 -4.55904 0.226274 -5.44096搜索% 17 -0.32 -4.39277 0.32 -5.60723搜索% 19 -0.452548 -4.18759 0.452548 -5.81241搜索% 21 -0.64 -3.98214 0.64 -6.01786搜索% 23 -0.905097 -3.93126 -0.905097 -6.06874搜索% 25 -1.28 -4.53715 1.28 -5.46285搜索% 27 -1.81019 -7.31125 -1.81019 -2.68875搜索% 29 -2.56 -16.6572 2.56 6.65722搜索

现在是经典的zeroin算法可以快速快速地找到零点。要查看详细信息，请运行fzerogui从<一个href="//www.tatmou.com/blogs/matlabcentral/fileexchange/37976-numerical-computing-with-matlab">MATLAB数值计算然后点击汽车按钮。用割线法和逆二次插值法求零点。

% start -2.5600000000000001% start 2.5600000000000001% sec 1.0980323260716793% secant 1.7832168816106038% iqi 2.2478393639958036% sec 2.0660057758331045% secant 2.0922079131171945iqi 2.0945566700001779% secant 2.0945514746903111% secant 2.0945514815423065iqi 2.0945514815423265%小2.0945514815423274

所以，即使最初的猜测很糟糕，fzero通过了德摩根测试。

定点迭代<一个name="e070bbec-6025-483d-be18-5709beecf91a">

我喜欢称其为“世界上最简单”的“WS”算法。试着找到映射$G(x)$的一个固定点。$$ x = G(x) $$迭代是$$ x_{n+1} = G(x_n) $$我应该为$G$选择什么?最明显的选择是方程$$ x = \frac{1}{2} (x^3-5) $$但在接近零的地方，$G(x)$的斜率太大，WS迭代发散。另一种选择是$$ x = \sqrt[3]{2x+5} $$这是可行的。当我们在做的时候，也要制作一个情节。

G = @ (x) (2 * x + 5) ^ (1/3) ezplot (G, 1[3])线(3 [1],[1 3],“颜色”，“k”) dkgreen = [0 2/3 0];X = 1.4;Z = 0;disp (x)而Abs (x-z) > eps(x) z = x;x = G(x);Line ([z z]，[z x]，“颜色”，dkgreen) line([z x]，[x x]，“颜色”dkgreen) disp (x)结束

G = @(x)(2*x+5).^(1/3) 1.400000000000000 1.983192482680775 2.077490885128178 2.091955753470501 2.094515696278154228 2.09454550097140211 2.094551271169830 2.09455144951430152 2.094551481525281 2.094551481539736 2.094551481542317 2.09455148152325 2.09455148152326 2.0945514815232327

这需要很长时间。它只是线性收敛。但我们还是挤过了德摩根的大门。

根平方<一个name="2a5091a4-09f6-4821-a0c2-35929c9389b3">

Graffe平方根法是由Germinal Pierre Dandelin在1826年，Nikolai Lobachevsky在1834年，Karl Heinrich Graffe在1837年独立发明的。下面引用的Alston Householder的一篇文章详细介绍了谁发明了什么。其思想是将一个多项式的系数转换为另一个多项式的零点是原多项式的平方。如果0在大小上被很好地分开，重复这个过程将最终产生一个高次多项式，它的前两项提供了一个与主0很好的近似。该方法在手工计算中很受欢迎，但当涉及到可靠的自动实现时，存在严重的困难。重复的零、等大小的复杂零，特别是为防止浮点溢出和下溢而进行的缩放都是障碍。但它适用于历史三次幂，它的主子是0。假设$p(x)$是一个三次多项式。$$ p(x) = p_1 x^3 + p_2 x^2 + p_3 x + p_4 $$重新排列方程$p(x) = 0$，使$x$的奇数次项在等号的一边，偶数次项在另一边。$ p_1 x^3 + p_3 x = - (p_2 x^2 + p_4)对方程两边平方，收集项，形成一个新的多项式，$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.

函数Q = graffe(p);% grafff平方根。Q = graffe(p)。% p和q是四向量，系数为三次。%根(q) =根(p) ^2。Q = 0(大小(p));Q (1) = p(1)^2;Q (2) = -(p(2)²- 2*p(1)*p(3));Q (3) = p(3)²- 2*p(2)*p(4);Q (4) = -p(4)^2;结束

我们试一下历史多项式x^3-2x-5。系数的浮点指数很快开始在每一步翻倍。幸运的是，我们在八个步骤中达到了我的停止标准。如果必须在不重新缩放的情况下尝试更多步骤，则会溢出。但此时，我们求的0的256次幂完全压倒了其他根的幂，所以第一个系数之比的256次方根提供了一个精确到完全精确的结果。

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;而Max (abs(p)) <√(realmax) m = 2*m;P = graffe(P);Z = (-p(2)/p(1)).^(1/m);流(fmt, m, p, z);结束

11 0 -2 -5 0.000000000000000 21 -4 4 -25 2.000000000000000 41 -8 -184 -625 1.681792830507429 81 -1.3891e+05 2.135184796196703 161 -1.8833e+10 1.125e+16 -2.3283e+22 2.094553553759747412 64 1 -3.5467e+20 -7.5043e+32 -5.421e+44 2.094551481347086 128 1 -1.2579e+41 1.7861e+65 - 2.94551481542327 256 1 -1.5824e+82 -4.2032e+130 -8.6362e+178 2.094551481542327

这个基本的平方根实现通过了De Morgan的闸门。

倒数第二余数简化<一个name="8d205b79-eb44-4629-aed6-a196eb1ba0c7">

大约两年前，我写了一篇关于<一个href="//www.tatmou.com/blogs/cleve/2013/01/21/reduced-penultimate-remainder">简化倒数余数算法．这是一个本科生的个人研究项目，我现在还满怀深情地记得，它是我见过的最晦涩、最不切实际的算法。该算法试图找到一个低次多项式，它能精确地除除我们所寻找的零的多项式，没有任何余数。结果除数的0就是被除数的一些0。在那篇文章中，我试图找到$x^3-2x-5$的线性因子。这就得到了我在这篇文章中计算的0。但无论给出什么起始因子，算法都不能收敛。那计算二次因子呢?如果成功，这将产生两个零。这是只需要一步的函数。

函数Q = rpr(p, Q)% RPR减少倒数第二余数。% r = rpr(p,q)对于多项式p和q。假设p(1) = 1, q(1) = 1。%多项式长除法。当r的度= q的度时停止。N =长度(p);M =长度(q);R = p(1:m);为K = m+1:n r = r - r(1)*q;R = [R (2:end) p(k)];结束Q = r/r(1);结束

起始因子$q(x)$可以是任何虚数为零的二次元。我们试试q(x) = x^2+x+1。

% p = [1 0 -2 -5];% q = [1 1 1]% k = 1;% r = 0*q;当max(abs(q-r)) > 10*eps(max(abs(q)))% r = q;% q = rpr(p,q)% k = k+1;%结束% k

以下是114步中的前10步和后4步。我删掉了中间的100个步骤。

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.1814671814671821.000000000000000 2.119373776908023 2.534246575342465 . . . . . . . . .1.000000000000000 2.094551481542332 2.387145908831160 1.000000000000000 2.094551481542323 2.387145908831149 1.000000000000000 2.09455149088312327 2.387145908831159 1.000000000000000 2.094551481542328 2.387145908831155 k = 114

现在我可以用反褶积来做多项式除法p(x)/q(x)。这就产生了一个线性商和期望的零点。

% r = deconv(p,q)% x1 = -r(1)

R = 1.000000000000000 -2.094551481542328 x1 = 2.094551481542328

参考文献<一个name="8eed4963-6fe8-4357-9327-ea639cf15cb0">

阿尔斯通·s·豪斯霍尔德《丹德林，洛巴切夫斯基，或格拉夫》美国数学月刊， vol. 66, 1959, pp. 464-466。<<一个href="http://www.jstor.org/stable/2310626">http://www.jstor.org/stable/2310626 http://www.jstor.org/stable/2310626>

获取MATLAB代码

MATLAB®R2015a发布