牛顿迭代通常用于求解实方程，但公式本身并不要求变量是实数。把复平面上的每个点都作为初值运行一次牛顿迭代，再按照最后收敛到的根着色，就能得到一张牛顿分形（Newton fractal）。

同一个多项式的几个根会像不同国家一样瓜分复平面。吸引域内部相当平静，边界却会不断分叉；在边界附近稍微改变初值，迭代就可能收敛到另一个根。

复平面上的牛顿迭代

对于多项式 $p(z)$，牛顿迭代为

$$
z_{k+1}=N_p(z_k)=z_k-\frac{p(z_k)}{p’(z_k)}.
$$

多项式的根是牛顿映射的不动点，根附近的初值会快速向它收敛。所有最终收敛到同一个根的初值组成该根的吸引域，不同吸引域的公共边界则是牛顿映射的 Julia 集。

牛顿分形例子

最经典的例子是 $p(z)=z^3-1$，它的三个根位于单位圆上。

在绘图中，三个根各自分配一种颜色，颜色亮度表示收敛速度：

  • 远离边界，尤其是在根附近时，通常只需几次迭代，颜色较亮；
  • 越接近边界，迭代路线越曲折，颜色也越暗。

对于 $p(z)=z^n-1$，它的 $n$ 个根均匀分布在单位圆上，根的旋转对称性也会反映在吸引域上。

对于二次方程，镜像对称结构比较简单：

四次方程得到四重对称结构：

五次方程则得到五个吸引域：

六次方程结果如下：

随着次数增加，原点附近会出现更多狭长的交错区域。
但是牛顿迭代不保证对所有初值都收敛。原点是牛顿映射的极点，Julia 集上还有不收敛到根的特殊初值。
但是对于 $z^n-1$ 这类开方问题，不收敛集合的平面测度为零，因此更准确的说法是“几乎处处收敛到某个根”，而不是全局收敛。

需要注意的是：牛顿迭代的长期行为不一定只有发散或者收敛到某个根，还可能出现周期轨道等更复杂的结构。

例如 $p(z)=z^3-2z+2$ 的导数为 $p’(z)=3z^2-2$，把 0 和 1 代入开头的牛顿公式可得

$$
N_p(0)=0-\frac{p(0)}{p’(0)}=1, \quad
N_p(1)=1-\frac{p(1)}{p’(1)}=0.
$$

因此 0 与 1 构成一个二周期轨道，并且这个周期会吸引附近的一片初值。下面图中的黑色圆形区域并不是绘图精度不足，而是其中的点没有收敛到多项式的三个根，大量点最终进入了这个二周期。

注：这里实际绘图中的黑色像素只能解释为“在给定次数内没有收敛到已知根”。它可能来自周期轨道、导数为零、数值发散，也可能只是迭代次数不够，不能直接把所有黑色像素都视为 Julia 集。

Python 实现

数值绘制只需要以下几步：

  1. 在指定矩形内生成复数网格，每个像素对应一个初值；
  2. 对所有仍然活跃的像素同时执行牛顿迭代；
  3. 检查当前值是否进入某个根的容差邻域；
  4. 用根的编号选择颜色，用迭代次数控制亮度；
  5. 超过最大迭代次数仍未收敛的点保留为黑色。

这里我们使用对数方式把迭代次数映射为亮度，再与对应根的基础 RGB 颜色相乘。相比线性映射，这种做法可以保留快速收敛区域中的明暗层次。

完整代码如下，可以增加 width 和 height 以绘制高清图片，但是脚本的计算会更慢。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
"""使用牛顿迭代绘制多项式在复平面上的吸引域。"""

import matplotlib.pyplot as plt
import numpy as np


PALETTE = np.array(
    [
        [255, 0, 82],
        [255, 212, 0],
        [0, 198, 141],
        [100, 80, 255],
        [255, 140, 0],
        [0, 180, 255],
        [230, 80, 200],
        [160, 220, 60],
    ],
    dtype=float,
) / 255.0


# 修改这个字典即可绘制不同多项式和视口。
CONFIG = {
    "function": lambda z: z**3 - 1,
    "derivative": lambda z: 3 * z**2,
    "roots": np.exp(2j * np.pi * np.arange(3) / 3),
    "title": r"$p(z)=z^3-1$",
    # 复平面显示范围。
    "xlim": (-2.0, 2.0),
    "ylim": (-2.0, 2.0),
    # width、height 决定分形本身的细节数量。
    "width": 1800,
    "height": 1800,
    "max_iter": 150,
    "tolerance": 1e-8,
    "derivative_epsilon": 1e-14,
    "divergence_limit": 1e10,
    # 迭代越慢颜色越暗，取值应在 0 到 1 之间。
    "darkest_brightness": 0.30,
    "palette": PALETTE,
    # Figure 参数只影响带坐标轴版本；raw_path 严格保持 width × height。
    "figure_size": (9.0, 9.0),
    "display_dpi": 120,
    "save_dpi": 300,
    "interpolation": "nearest",
    "show_axes": True,
    "show_title": True,
    "figure_path": "newton-fractal.png",
    "raw_path": None,
    "show": True,
}


def compute_newton_fractal(config=None):
    """计算每个像素收敛到的根、迭代次数和 RGB 图像。"""
    cfg = CONFIG.copy()
    if config:
        cfg.update(config)

    function = cfg["function"]
    derivative_function = cfg["derivative"]
    roots = np.asarray(cfg["roots"], dtype=np.complex128)
    if not callable(function) or not callable(derivative_function):
        raise TypeError("function 和 derivative 必须是可调用对象")
    if roots.ndim != 1 or len(roots) == 0 or not np.isfinite(roots).all():
        raise ValueError("roots 必须是一维、非空且只包含有限复数")
    width = int(cfg["width"])
    height = int(cfg["height"])
    x = np.linspace(*cfg["xlim"], width)
    y = np.linspace(*cfg["ylim"], height)
    z = x[None, :] + 1j * y[:, None]

    root_index = np.full(z.shape, -1, dtype=np.int16)
    iteration_count = np.full(z.shape, cfg["max_iter"], dtype=np.int16)
    active = np.ones(z.shape, dtype=bool)
    tolerance_sq = cfg["tolerance"] ** 2

    for iteration in range(cfg["max_iter"] + 1):
        # 逐根比较比构造 (height, width, root_count) 数组更节省内存。
        nearest = np.zeros(z.shape, dtype=np.int16)
        minimum_distance_sq = np.full(z.shape, np.inf)
        for index, root in enumerate(roots):
            distance_sq = (z.real - root.real) ** 2 + (z.imag - root.imag) ** 2
            closer = distance_sq < minimum_distance_sq
            minimum_distance_sq[closer] = distance_sq[closer]
            nearest[closer] = index

        converged = active & (minimum_distance_sq < tolerance_sq)
        root_index[converged] = nearest[converged]
        iteration_count[converged] = iteration
        active[converged] = False

        if iteration == cfg["max_iter"] or not active.any():
            break

        derivative = derivative_function(z)
        safe = active & np.isfinite(derivative) & (
            np.abs(derivative) > cfg["derivative_epsilon"]
        )
        with np.errstate(over="ignore", divide="ignore", invalid="ignore"):
            z[safe] -= function(z[safe]) / derivative[safe]

        failed = active & (
            ~np.isfinite(z) | (np.abs(z) > cfg["divergence_limit"]) | ~safe
        )
        active[failed] = False

    palette = np.asarray(cfg["palette"], dtype=float)
    if len(palette) < len(roots):
        raise ValueError("palette 中的颜色数量不能少于多项式的根数")

    image = np.zeros((*z.shape, 3), dtype=float)
    known = root_index >= 0
    level = np.log1p(iteration_count[known]) / np.log1p(cfg["max_iter"])
    brightness = 1 - (1 - cfg["darkest_brightness"]) * np.minimum(level, 1)
    image[known] = palette[root_index[known]] * brightness[:, None]
    return image, roots, root_index, iteration_count


def demo(config):
    """按照配置字典绘制、保存并返回牛顿分形。"""
    cfg = CONFIG.copy()
    cfg.update(config)

    image, roots, root_index, iteration_count = compute_newton_fractal(cfg)
    fig, ax = plt.subplots(
        figsize=cfg["figure_size"],
        dpi=cfg["display_dpi"],
    )
    ax.imshow(
        image,
        extent=(*cfg["xlim"], *cfg["ylim"]),
        origin="lower",
        interpolation=cfg["interpolation"],
        resample=False,
    )
    ax.set_aspect("equal")

    if cfg["show_axes"]:
        ax.set_xlabel(r"$\operatorname{Re}(z)$")
        ax.set_ylabel(r"$\operatorname{Im}(z)$")
    else:
        ax.set_axis_off()
    if cfg["show_title"]:
        ax.set_title(cfg["title"])
    fig.tight_layout()

    if cfg["figure_path"]:
        fig.savefig(
            cfg["figure_path"],
            dpi=cfg["save_dpi"],
            bbox_inches="tight",
            pad_inches=0.05,
        )
    if cfg["raw_path"]:
        plt.imsave(cfg["raw_path"], image, origin="lower")

    if cfg["show"]:
        plt.show()
    else:
        plt.close(fig)
    return image, roots, root_index, iteration_count

使用不同配置即可绘制牛顿分形。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
cfg1 = {
    "function": lambda z: z**3 - 1,
    "derivative": lambda z: 3 * z**2,
    "roots": np.exp(2j * np.pi * np.arange(3) / 3),
    "title": r"$p(z)=z^3-1$",
    "xlim": (-2.0, 2.0),
    "ylim": (-2.0, 2.0),
    "figure_path": "newton-z3-minus-1.png",
}

demo(cfg1)

cfg2 = {
    "function": lambda z: z**4 - 1,
    "derivative": lambda z: 4 * z**3,
    "roots": np.exp(2j * np.pi * np.arange(4) / 4),
    "title": r"$p(z)=z^4-1$",
    "xlim": (-2.0, 2.0),
    "ylim": (-2.0, 2.0),
    "figure_path": "newton-z4-minus-1.png",
}

demo(cfg2)

cfg3 = {
    "function": lambda z: z**5 - 1,
    "derivative": lambda z: 5 * z**4,
    "roots": np.exp(2j * np.pi * np.arange(5) / 5),
    "title": r"$p(z)=z^5-1$",
    "xlim": (-2.0, 2.0),
    "ylim": (-2.0, 2.0),
    "figure_path": "newton-z5-minus-1.png",
}

demo(cfg3)

cfg4 = {
    "function": lambda z: z**3 - 2 * z + 2,
    "derivative": lambda z: 3 * z**2 - 2,
    "roots": np.array(
        [
            -1.7692923542386314,
            0.8846461771193155 + 0.5897428050222055j,
            0.8846461771193155 - 0.5897428050222055j,
        ]
    ),
    "title": r"$p(z)=z^3-2z+2$",
    "xlim": (-2.5, 2.5),
    "ylim": (-2.5, 2.5),
    "figure_path": "newton-z3-minus-2z-plus-2.png",
}

demo(cfg4)

cfg5 = {
    "function": lambda z: z**2 - 1,
    "derivative": lambda z: 2 * z,
    "roots": np.exp(2j * np.pi * np.arange(2) / 2),
    "title": r"$p(z)=z^2-1$",
    "xlim": (-2.0, 2.0),
    "ylim": (-2.0, 2.0),
    "figure_path": "newton-z2-minus-1.png",
}

demo(cfg5)

cfg6 = {
    "function": lambda z: z**6 - 1,
    "derivative": lambda z: 6 * z ** 5,
    "roots": np.exp(2j * np.pi * np.arange(6) / 6),
    "title": r"$p(z)=z^6-1$",
    "xlim": (-2.0, 2.0),
    "ylim": (-2.0, 2.0),
    "figure_path": "newton-z6-minus-1.png",
}

demo(cfg6)

参考：