使用历史悠久的三次多项式x ^ 3 - 2 x - 5美元来测试几zero-finding算法。

内容

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

我的<一个href="//www.tatmou.com/blogs/cleve/2015/12/21/a-historic-cubic/">以前的文章是关于一个脚注中一个经典的数值分析文本,惠塔克和罗宾逊是奥古斯都•德•摩根的妙语。“我之所以叫$ x ^ 3-2x-5 = 0美元方程是著名的因为它是沃利斯偶然表现出牛顿法首次出版时,由于每一个义务中的数值解算器都觉得他的一个例子。发明了一种数值方法,忽视对这个方程显示它是如何工作的,你是一位朝圣者不进来的小wicket (见j·班扬)。”So, let's try a few numerical solvers on this equation.

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

我必须首先尝试牛顿法。$ $间{n + 1} = x_n - \压裂{f (x_n)} {f ' (x_n)} $ $我知道导数和我有一个好开始猜测。

f = @ x (x)。^ 3 - 2 * * 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

这个达到很快的解决方案。只有5个步骤,第五第四已经确认正确完整的精度。但牛顿法需要导数的知识和一个很好的起点。这些都是普遍使用的缺点。

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

让我们看看零finder中可用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。设置显示参数要求中间结果。

%选择= optimset(“显示”、“iter”);

假装我们不知道任何关于这个函数,并开始搜索$ x = 0 $。

% fzero (f, 0,选择)

每次我们得到一行输出一个却降低了,b是增加了。所有的值f (a)是负面的和所有的值f (b)也是消极的,直到b = 2.56当第一个信号变化。

% 0周围寻找一个区间包含一个信号变化:% Func-count f b (a) (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搜索-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数值计算,然后点击汽车按钮。割线法和逆二次插值法用于发现零。

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

所以,即使是开始一个糟糕的初始猜测,fzero通过德摩根的测试。

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

我喜欢称之为“WS”为“世界最简单的”算法。试图找到一个固定的点映射G (x)美元。$ $ x = G (x) $ $ $ $的迭代间{n + 1} = G (x_n) $ $ $ G $我应该选择什么?最显而易见的选择是方程$ $ x = \压裂{1}{2}(x ^ 3 - 5) $ $但接近零的斜率G (x)太大美元和WS迭代发散的。另一个选择是$ $ x = \ sqrt [3] {2 x + 5} $ $这将工作。而我们,也产生一个阴谋。

G = @ (x) (2 * x + 5) ^ (1/3) ezplot (G, 1[3])线(3 [1],[1 3],“颜色”,“k”)dkgreen = [0 0] 2/3;x = 1.4;z = 0;disp (x)而abs (x z) > eps (x) z = x;x = G (x);x线([z z]、[z],“颜色”dkgreen)线(x [z], [x],“颜色”dkgreen) disp (x)结束

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

这需要很长时间。它只是线性收敛。但无论如何我们挤过德摩根的wicket。

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

“Graffe”root-squaring方法由胚皮埃尔Dandelin独立发明于1826年,1834年尼古拉Lobachevsky,卡尔·海因里希Graffe在1837年。一篇文章由阿尔斯通户主引用下面进入细节谁发明了什么。这个想法是一个多项式的系数转化为另一个零的原始的广场。如果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) $ $广场方程两边并收集形成新的多项式,问(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);% Graffe根平方。q = graffe (p)。% p和q 4-vectors三次曲线的系数。%根(q) =根(p) ^ 2。q = 0(大小(p));q (1) = p (1) ^ 2;问(2)= - (p (2) ^ 2 - 2 * p (1) * p (3));问(3)= p (3) ^ 2 - 2 * p (2) * (4);问(4)= - p (4) ^ 2;结束

让我们试穿历史性多项式$ x ^ 3-2x-5美元。的浮点指数系数每一步很快就翻倍。幸运的是,我们在八个步骤达到我的停止准则。如果我们尝试更多的步骤没有重新调节,我们会溢出。但此时,零的256次方我们完全主宰的权力其他根后,所以的第256根的比例只是第一个系数提供了一个结果,是准确完整的精度。

p = [1 0 2 5];fmt =' % 5.0 f % 5.0克% 13.5 g % 13.5克% 13.5 g % 20.15 f \ n ';流(fmt, 1 p, 0) m = 1;而马克斯(abs (p) <√(最大浮点数)m = 2 * m;p = graffe (p);z = (- p(2) /(1)页)。^ (1 / m);流(fmt, m, p, z);结束

1 1 0 2 5 0.000000000000000 - 2 1 4 4 -25 1.681792830507429 2.000000000000000 - 4 1 8 -184 -625 8 -432 23856 -3.9063 e + 05年16 1 -1.3891 2.135184796196703 2.3161 e + e + 05年08年-1.5259 e + 11 2.096144583759898 32 1 -2.3283 -1.8833 1.125 e + e + 10 16 e + 22 64 1 -3.5467 2.094553557477412 -5.421 e + 20 -7.5043 e + 32 e + 44 128 1 -1.2579 2.094551481347086 -2.9387 e + e e + 1.7861 + 41 65 89 256 2.094551481542327 82 -4.2032 -1.5824 e + e + e + 178 -8.6362 130 2.094551481542327

这个基本的实现root-squaring使它通过德摩根的wicket。

减少剩余倒数第二<一个name="8d205b79-eb44-4629-aed6-a196eb1ba0c7">

大约两年前,我在博客中写道<一个href="//www.tatmou.com/blogs/cleve/2013/01/21/reduced-penultimate-remainder">减少了倒数第二部分算法。这是一个大学生的个人研究项目,现在我天真地记得最晦涩难懂的,我甚至见过不切实际的算法。算法试图找到一个多项式分低,没有剩余,我们所寻求的多项式的0。产生的零因子的一些零股息。在那篇文章为例,我试图找到一个美元的线性因子x ^ 3-2x-5美元。会给我零,我已经计算在这篇文章中。但是算法不收敛,不管什么因素开始。如何计算一个二次因子呢?如果成功,这将产生两个零。这是函数需要一个步骤。

函数q = rpr (p, q)% RPR减少剩余倒数第二。% r = rpr (p, q), p和q为多项式。%认为p (1) = 1, (1) = 1。%多项式长除法。%停止当r =度程度的问。n =长度(p);m =长度(q);r = p (1: m);为k = m + 1: n r = r - r (1) * q;r = [r(结束)2:p (k)];结束q = r / r (1);结束

起始因子q (x)美元可以是任何二次虚构的0。让我们尝试问美元(x) = x ^ 2 + x + 1美元。

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

这是第十和114年最后四个步骤。我编辑了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.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

我现在可以使用反褶积多项式除法$ p (x) / (x)美元。这产生一个线性系数和所需的零。

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

r = 1.000000000000000 - -2.094551481542328 x1 = 2.094551481542328

引用<一个name="8eed4963-6fe8-4357-9327-ea639cf15cb0">

阿尔斯通户主,“Dandelin Lobacevskii,或Graeffe”,美国数学月刊66卷,1959年,页464 - 466。<<一个href="http://www.jstor.org/stable/2310626">http://www.jstor.org/stable/2310626 http://www.jstor.org/stable/2310626>

得到了MATLAB代码

发表与MATLAB®R2015a