牛顿分形

<image>

牛顿分形（Newton fractal）是将牛顿法应用于一给定多项式p(Z) ∈ ℂZ或超越函数而得到的复平面上的一个边界集。它是由牛顿法所定义的亚纯函数z ↦ z − p(z)的朱利亚集。当不存在吸引循环（阶数大于1）时，它将复平面划分为不同的区域Gk，每个区域与多项式的根ζk相关联，其中k = 1, …, deg(p)。此时牛顿分形类似于曼德博集合，并且与其他分形一样，它将简单的数学描述变成了非常繁复的图像。从数值分析的角度而言，牛顿分形表现出牛顿法在二次收敛区域之外对于初始点的选择非常敏感。

将复平面上的某一点作为牛顿法迭代zn + 1 := zn − p(zn)的初始点z0，可以通过迭代得到一个点序列z1, z2, …,。如果这一序列收敛于根ζk，则将z0划入区域Gk。如此便能将复平面上的这一点与多项式的某一个根相对应。不过值得注意的是，对于二次以上的多项式，都存在一些点会使得牛顿迭代无法收敛到任何根上，例如不同根的吸引域的边界。甚至存在一些多项式，某些开集中的任意初始点都无法收敛到任何根上。一个简单的例子是z3 − 2z + 2，某些点会被吸引到循环0、1、0、1……中，而不被任何根所吸引。

如果以一个开集中的任意点为初始点，迭代最终都收敛于某一根或循环，则该集合是这一牛顿迭代的法图集。一个法图集对应于一个根或循环。所有这些法图集的并集与朱利亚集为互补集。这一朱利亚集即是法图集的共同边界。因此，朱利亚集中的每个点都是每个法图集的一个聚点。正是由于这一性质导致了朱利亚集的分形结构（当多项式的次数大于2时）。

为了绘制一个牛顿分形图像，可以首先选择指定数量d的复点(ζ1, …, ζd)并计算多项式的系数(p1, …, pd)

$p(z)=z^d+p_1z^{d-1}+\cdots+p_{d-1}z+p_d:=(z-\zeta_1)(z-\zeta_2)\cdots(z-\zeta_d)$ .

于是，对于复平面上的一个矩形网格

$$z_{mn} = z_{00} + m \, \Delta x + in \, \Delta y; \quad m = 0, \ldots, M - 1; \quad n = 0, \ldots, N - 1$$

找到每个点(m,n)对应的根ζk(m,n)的编号k(m,n)，并通过为每一点分配一个颜色fk(m,n)来填充这一M × N网格。另外，颜色可以取决于距离D(m,n)。对于某一固定的小ε > 0，距离D可以定义为第一个使得zD − ζk(m,n) < ε成立的D值。

## 牛顿分形的推广

牛顿迭代的一种推广可表示为

$$z_{n+1}=z_n- a \frac{p(z_n)}{p'(z_n)}$$

其中a是任意复数。当a = 1时即对应于牛顿分形。当a位于以1为圆心、半径为1的圆盘以内时，该映射的不动点是稳定的。而当a位于这个圆盘之外时，不动点是局部不稳定的，不过该映射仍能表现出朱利亚集的分形结构。如果p是d次多项式，则当是a位于以d为圆心、半径为d的圆盘内时，序列zn是有界的。

更一般地，牛顿分形是朱利亚集的一个特例。

<image>

- p(z) = z3 − 1的牛顿分形，颜色表示需要的迭代步数

<image>

- p(z) = z3 − 1的牛顿分形，颜色表示最终收敛到的根

<image>

- p(z) = z3 − 2z + 2的牛顿分形，红色区域中的点不收敛到任何根

<image>

- 一个7次多项式的牛顿分形，颜色表示最终收敛到的根，深浅表示收敛速度

<image>

- p(z) = z8 + 15z4 − 16的牛顿分形

<image>

- p(z) = z5 − 3iz3 − (5 + 2i)z2 + 3z + 1的牛顿分形，颜色表示最终收敛到的根，深浅表示迭代步数

<image>

- p(z) = sin z的牛顿分形，颜色表示最终收敛到的根，深浅表示迭代步数

<image>

- p(z) = sin z的牛顿分形

<image>

- p(z) = z3 − 1的广义牛顿分形（a = −1）

<image>

- p(z) = z2 − 1的广义牛顿分形（a = 1 + i）

<image>

- p(z) = z3 − 1的广义牛顿分形（a = 2）

<image>

- p(z) = z4 + 3i − 1的广义牛顿分形（a = 2.1）

<image>

- p(z) = z6 + z3 − 1的牛顿分形

<image>

- p(z) = sin z − 1的牛顿分形

<image>

- p(z) = cosh z − 1的牛顿分形

### 新星分形

新星分形（Nova fractal）是由Paul Derbyshire于1990年代发明的一种分形。它也是牛顿分形的一种推广，即在每一步迭代时都增加了一个c值：

$$z_{n+1}=z_n- a \frac{p(z_n)}{p'(z_n)} + c = G(a, c, z)$$

新星分形的“朱利亚”版本令c为常数，并将像素坐标设为初始值z0。而新星分形的“曼德博”版本则用像素坐标来初始化c，并将z0设置为一临界点，满足

$$\frac{\partial}{\partial z} G(a, c, z) = 0.$$

例如，p(z) = (z − 1)3的临界点位于z = 1处。

## 实现

为了用计算机实现牛顿分形，需要有一个起始函数及其导函数：

$$\begin{align} f(z) &= z^3 - 1 \\ f'(z) &= 3z^2 \end{align}$$

该函数的三个根是

$$z = 1,\ -\tfrac12 + \tfrac{\sqrt 3}{2}i,\ -\tfrac12 - \tfrac{\sqrt 3}{2}i$$

该函数可以以伪代码表示如下：

```c
//z^3-1 
float2 Function (float2 z)
{
	return cpow(z, 3) - float2(1, 0); //cpow is an exponential function for complex numbers
}

//3*z^2
float2 Derivative (float2 z)
{
	return 3 * cmul(z, z); //cmul is a function that handles multiplication of complex numbers
}
```

之后只需用给定函数实现牛顿法即可：

```c
float2 roots[3] = //Roots (solutions) of the polynomial
{
	float2(1, 0), 
	float2(-.5, sqrt(3)/2), 
	float2(-.5, -sqrt(3)/2)
};
	
color colors[3] =  //Assign a color for each root
{
    red,
    green,
    blue
}

For each pixel (x, y) on the target, do:
{
	zx = scaled x coordinate of pixel (scaled to lie in the Mandelbrot X scale (-2.5, 1))
    zy = scaled y coordinate of pixel (scaled to lie in the Mandelbrot Y scale (-2, 1))

    float2 z = float2(zx, zy); //z is originally set to the pixel coordinates

	for (int iteration = 0;
	     iteration < maxIteration;
	     iteration++;)
	{
		z -= cdiv(Function(z), Derivative(z)); //cdiv is a function for dividing complex numbers

        float tolerance = 0.000001;
        
		for (int i = 0; i < roots.Length; i++)
		{
		    float2 difference = z - roots[i];
		    
			//If the current iteration is close enough to a root, color the pixel.
			if (abs(difference.x) < tolerance && abs (difference.y) < tolerance)
			{
				return colors[i]; //Return the color corresponding to the root
			}
		}
		
    }
    
    return black; //If no solution is found
}
```