在Mathematica中实现Newton-Raphson迭代的收敛时间算法

对于给定的比较简单的三次多项式，可以使用FixedPoint直接给出不动点，进而在一定程度上回答Cayley问题。但是这里不能够给出收敛速度。作为类比，这里的收敛速度，类似于绘制Mandelbrot集合的逃逸时间，可以用来对图形进行染色。

收敛时间算法：对于给定的一个点，将计数以此点为初始点的迭代收敛速度，并将相应的计数划分为从0到255的三个等分区间中。

NewtonCounter = Compile[{{z, _Complex}},Module[{counter = 0, zold = N[z] + 1.0, 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 = 2*zold/3 + 1/(3*zold^2))];Which[Abs[znew - 1] < 10^(-4), counter,Abs[znew + 0.5 - 0.866025 I] < 10^(-4), 85 + counter,Abs[znew + 0.5 + 0.866025 I] < 10^(-4), 170 + counter,True, 255]]];

使用函数NewtonCounter，对于平面中的某个区域，可以构造一个数表：

NewtonArray1[{{remin_, remax_}, {immin_, immax_}}, steps_] :=Table[NewtonCounter[x + y I],{y, immin, immax, (immax - immin)/steps},{x, remin, remax, (remax - remin)/steps}]

这样，就可使用ListDenstiyPlot函数进行绘图了：

region1 = NewtonArray1[{{-2, 2}, {-2, 2}}, 1000]; ListDensityPlot[region1]