在Mathematica中实现Newton-Raphson迭代的收敛时间算法（一般三次多项式）

考察一般的三次多项式，以r为参数：

p[z_, r_] := z^3 + (r - 1) z - r;
roots[r_] := z /. Solve[p[z, r] == 0, z]；

此多项式的根为：

尽管看起来这个多项式是特殊的，其实一般的三次多项式都是可以通过线性变换化为这个形式的。此时的Newton-Raphson映射为：

此时，可以编写带有参数r的收敛速度计数函数：

NewNewtonCounter = Compile[{{z, _Complex}, {r, _Real}, {otherroot, _Complex}},Module[{counter = 0, zold = N[z] + 1, znew = N[z]},If[Abs[znew] < 10^(-9), znew = 10^(-9) + 0.0 I,znew = znew];For[counter = 0,(Abs[zold - znew] > 10^(-6)) && (counter < 85), counter++,(zold = znew; znew = (r + 2*zold^3)/(-1 + r + 3*zold^2))];Which[Abs[znew - 1] < 10^(-4), counter,Abs[znew - otherroot] < 10^(-4), 85 + counter,Abs[znew - Conjugate[otherroot]] < 10^(-4), 170 + counter,True, 255]]];(*返回给定区域中的每个点相应的收敛计数*)  NewNewtonArray[r_, {{remin_, remax_}, {immin_, immax_}}, steps_] :=Module[{croot = -N[(1 + Sqrt[1 - 4 r])/2]},Table[NewNewtonCounter[x + y I, r, croot],{y, immin, immax, (immax - immin)/steps},{x, remin, remax, (remax - remin)/steps}]]

region = NewNewtonArray[2, {{-2, 2}, {-2, 2}}, 1000];
NewtonPlot1[{{-2, 2}, {-2, 2}}, region, NewtonColorRGB]

region = NewNewtonArray[0.05, {{-2, 2}, {-2, 2}}, 1000];
NewtonPlot1[{{-2, 2}, {-2, 2}}, region, NewtonColorRGB]