1. 项目概述：当牛顿法遇上分形

如果你学过数值分析，一定对牛顿法不陌生。这个以艾萨克·牛顿命名的迭代方法，是寻找方程根的一把利器，以其超线性的收敛速度著称。教科书上通常会给你一个光滑的二次函数曲线，画上几条切线，然后告诉你：“看，经过几次迭代，我们就找到了根。”整个过程显得优雅、确定，甚至有些枯燥。但今天，我们要聊的恰恰是牛顿法“不优雅”、“不确定”甚至“混乱”的另一面——牛顿法分形。

这听起来可能有些矛盾。牛顿法不是用来找确定解的吗？怎么会和描绘无限复杂、自相似图案的分形扯上关系？秘密就藏在“初始值”的选择里。对于一个简单的多项式方程，比如 z^3 - 1 = 0 ，它在复数域上有三个根。如果你从复平面上的任意一点出发，应用牛顿法进行迭代，它最终会收敛到这三个根中的一个。但问题是， 从哪一点出发，会收敛到哪一个根？

直觉可能会告诉你，离哪个根近，就会收敛到哪个根。事实也大致如此，但边界绝非清晰明了。当你用计算机为复平面上的每一个像素点（代表一个初始复数）着色，标记其牛顿迭代最终收敛的根时，你会惊讶地发现，收敛域的边界并非光滑曲线，而是错综复杂、蜿蜒曲折，并且具有无限精细的结构——这就是牛顿法分形，也称为牛顿盆地（Newton Basins）。它直观地揭示了牛顿法对初始值的极端敏感性，以及确定性算法在复数域上可以产生何等惊人的混沌与美感。

这个项目适合任何对数学可视化、计算科学或创意编程感兴趣的人。你不需要是数学博士，只要对复数、迭代和绘图有基本了解，就能亲手创造出这些令人着迷的图案。通过它，你不仅能深入理解牛顿法的收敛行为，还能掌握将数学算法转化为艺术图像的核心技巧。

2. 核心数学原理与算法拆解

要生成牛顿法分形，我们必须先吃透两个核心：牛顿迭代公式在复数域上的形式，以及如何将迭代结果映射为图像上的颜色。

2.1 复数域上的牛顿迭代法

牛顿法的核心思想是用切线逼近。对于实函数 f(x) ，公式是 x_{n+1} = x_n - f(x_n)/f'(x_n) 。在复数域上，变量 x 变成了复数 z ，但公式形式完全一致： z_{n+1} = z_n - f(z_n) / f'(z_n)

我们以最经典的例子 f(z) = z^3 - 1 为例。它的导数 f'(z) = 3z^2 。因此，迭代公式为： z_{n+1} = z_n - (z_n^3 - 1) / (3z_n^2)

这个方程在复数域上有三个根：

  1. z = 1
  2. z = -1/2 + i√3/2
  3. z = -1/2 - i√3/2

这三个根均匀分布在复平面的单位圆上，夹角为120度。牛顿法的任务就是从任意初始点 z_0 出发，通过迭代逼近其中一个根。

2.2 收敛判定与着色逻辑

算法在实现时，有几个关键点需要确定：

  1. 迭代终止条件 ：我们不可能无限迭代下去。通常设置两个条件：

    • 最大迭代次数 ：例如 100 次或 200 次，防止不收敛的点陷入无限循环。
    • 收敛容差 ：当 |f(z_n)| 或 |z_{n+1} - z_n| 小于一个极小值（如 1e-6 ）时，认为已收敛到根。

  2. 根的匹配与着色 ：迭代终止后，我们需要判断当前的 z_n 最接近哪一个根。计算 z_n 到三个根的欧几里得距离，归属到距离最近的那个根。然后，为每个根分配一种颜色（例如，根1红色，根2绿色，根3蓝色）。该像素点的颜色就由它归属的根决定。

  3. 迭代次数着色（增强艺术效果） ：这是让分形图案产生丰富渐变和细节的关键。除了用根来分配基础色（色调），我们还可以用收敛所需的迭代次数来调制颜色的明度或饱和度。例如，迭代次数少的区域更亮，迭代次数多（在边界附近挣扎）的区域更暗。这能清晰地揭示出收敛域的复杂边界。

核心逻辑伪代码 ： 

对于图像中的每个像素点 (x, y):
    // 将像素坐标映射到复平面区域，如 [-2, 2] x [-2, 2]
    z0 = complex(x_scaled, y_scaled)
    z = z0
    iteration = 0
    while iteration < max_iterations:
        如果 |f(z)| < tolerance: 跳出循环
        z = z - f(z) / f'(z) // 牛顿迭代
        iteration += 1

    // 判断收敛到哪个根
    root_index = 找出离z最近的根的索引
    // 根据 root_index 决定基础色 (hue)
    // 根据 iteration 决定亮度/饱和度 (value/saturation)
    像素颜色 = HSV_to_RGB(hue, saturation, value)

3. 实战：从零生成你的第一幅牛顿分形图

理论说得再多，不如动手一试。我们将使用 Python 的 NumPy 和 Matplotlib 库来实现，这是科学计算和可视化的黄金组合。

3.1 环境准备与依赖安装

确保你的 Python 环境已就绪。通过 pip 安装必要的库： 

pip install numpy matplotlib

NumPy 将允许我们以向量化方式对整个像素网格进行高效计算，这比用 for 循环遍历每个像素快几个数量级。 Matplotlib 则用于将数值数组渲染为图像。

3.2 代码实现步步解析

下面是一个完整、可运行的代码示例，我们将逐块解析。 

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import hsv_to_rgb

# 1. 定义图像参数和复平面区域
width, height = 1600, 1600  # 图像分辨率，越高细节越丰富
x_min, x_max = -2.0, 2.0
y_min, y_max = -2.0, 2.0

# 创建代表复平面网格的数组
x = np.linspace(x_min, x_max, width)
y = np.linspace(y_min, y_max, height)
X, Y = np.meshgrid(x, y)
Z = X + 1j * Y  # 关键的复数网格，包含了所有初始点

# 2. 定义目标函数 f(z) = z^3 - 1 及其导数
def f(z):
    return z**3 - 1

def df(z):
    return 3 * z**2

# 3. 定义三个根
roots = np.array([1, -0.5 + 0.5j*np.sqrt(3), -0.5 - 0.5j*np.sqrt(3)])
# 为每个根分配一个色调 (Hue)，范围在 [0, 1)
root_hues = np.array([0.0, 0.33, 0.67])  # 红，绿，蓝的HSV色调近似值

# 4. 牛顿迭代参数
max_iter = 50
tolerance = 1e-6

# 初始化迭代计数和根归属数组
iter_count = np.zeros(Z.shape, dtype=int)
root_index = np.full(Z.shape, -1, dtype=int)  # -1 表示未收敛

# 复制初始网格用于迭代
z = Z.copy()

# 5. 向量化牛顿迭代（核心循环）
for i in range(max_iter):
    # 计算当前点的函数值和导数值
    fz = f(z)
    dz = df(z)
    
    # 牛顿迭代步： z = z - f(z)/f'(z)
    # 注意：防止除以零，给导数加上一个极小值
    z_next = z - fz / (dz + 1e-15)
    
    # 判断哪些点已经收敛 (|f(z)| < tolerance)
    converged = np.abs(fz) < tolerance
    # 对于本次迭代新收敛的点，记录其迭代次数和最终的z值
    newly_converged = converged & (root_index == -1)
    iter_count[newly_converged] = i
    # 暂时将收敛点的最终z值存到root_index数组中（后续再匹配具体根）
    # 这里巧用root_index存储复数，需要先转换为浮点复数类型
    if newly_converged.any():
        # 我们需要另一个数组来存储收敛值，这里用一个小技巧
        # 实际上，我们应该把收敛的z值存到另一个变量。为了简化，我们稍后统一匹配。
        pass 
    
    # 更新z值
    z = z_next
    
    # 提前终止：如果所有点都收敛了，就跳出循环
    if np.all(converged):
        break

# 6. 迭代结束后，为所有点匹配最近的根
# 计算网格上每个点（最终z值）到三个根的距离
# 注意：对于未收敛的点，使用的是最后一次迭代的z值
distances = np.abs(z[:, :, np.newaxis] - roots)  # 利用广播机制
root_index = np.argmin(distances, axis=2)  # 每个点对应最近根的索引

# 7. 基于根和迭代次数创建HSV颜色图像
# H（色调）由根决定
H = root_hues[root_index]
# S（饱和度）可以设为固定值或与迭代次数相关，这里设为1（全饱和）
S = np.ones_like(H)
# V（明度）由迭代次数决定：迭代越多，越暗
# 将迭代次数归一化到 [0, 1] 区间，然后取反（次数多则V小，变暗）
V = 1.0 - (iter_count / max_iter)
# 对于未达到最大迭代就收敛的点，V值会较高（较亮）
V = np.clip(V, 0.1, 1.0)  # 限制明度范围，避免全黑

# 将HSV转换为RGB
HSV_image = np.stack([H, S, V], axis=-1)
RGB_image = hsv_to_rgb(HSV_image)

# 8. 绘制图像
plt.figure(figsize=(10, 10))
plt.imshow(RGB_image, extent=[x_min, x_max, y_min, y_max])
plt.axis('off')  # 不显示坐标轴
plt.title("Newton Fractal for f(z) = z^3 - 1", fontsize=16)
plt.tight_layout()
plt.savefig('newton_fractal_z3-1.png', dpi=300, bbox_inches='tight')
plt.show()

关键操作解析：

  1. 网格创建 ( np.meshgrid ) ：这是向量化计算的基础。我们一次性生成了包含所有初始复数的数组 Z ，避免了低效的 Python 层循环。
  2. 向量化迭代 ：整个 for 循环是对整个复数网格进行操作。 f(z) 和 df(z) 的调用会返回整个数组。这是性能的关键，使得计算1600x1600的网格也只需秒级完成。
  3. 收敛判断 ： np.abs(fz) < tolerance 会产生一个布尔数组，标记所有已收敛的点。我们利用 root_index == -1 来确保只记录第一次收敛时的迭代次数。
  4. 根匹配的广播技巧 ： distances = np.abs(z[:, :, np.newaxis] - roots) 这行代码是 NumPy 广播机制的经典应用。它将形状为 (H, W) 的 z 数组增加一个维度变成 (H, W, 1) ，然后减去形状为 (3,) 的 roots 数组，NumPy 会自动将 roots 广播到 (H, W, 3) ，并逐元素相减，最终得到每个点到三个根的距离数组 (H, W, 3) 。 np.argmin(axis=2) 则沿着第三个维度（根索引的维度）找到最小距离的索引。
  5. HSV色彩空间 ：使用色调（H）、饱和度（S）、明度（V）模型比直接使用RGB更直观。H 决定属于哪个根（颜色族），V 由迭代次数决定，能完美展现边界处的细腻结构。最后通过 hsv_to_rgb 转换。

运行这段代码，你就能得到一幅绚丽而复杂的牛顿分形图。在平滑的区域，颜色均匀，表示牛顿法快速收敛；而在那些蜿蜒曲折、如同闪电或海藻的边界区域，颜色深浅交织，揭示了初始值微小的变动就可能导致收敛到不同根的混沌特性。

4. 高级探索与参数调优

生成基础分形只是开始。通过调整参数和函数，你可以创造出无限多样的图案。

4.1 探索不同的多项式函数

z^3 - 1 只是冰山一角。尝试不同的多项式，收敛域的图案会截然不同。

  • z^4 - 1 ：有四个根（1, -1, i, -i）。你会看到四个主要的收敛盆地，它们的边界交织成更复杂的四重对称图案。 
    def f(z):
    return z**4 - 1
def df(z):
    return 4 * z**3
roots = np.array([1, -1, 1j, -1j])
root_hues = np.array([0.0, 0.25, 0.5, 0.75]) # 分配四种色调
  • z^5 - 1 ：五个根在单位圆上呈五边形分布，产生具有五重旋转对称性的分形。
  • z^3 - 2z + 2 ：这个函数有复数根。它的牛顿分形没有中心对称性，图案更加抽象和非对称，挑战你的视觉预期。

注意 ：当多项式的根不是简单的单位根时，你需要准确计算或定义这些根，用于后续的着色匹配。可以使用 np.roots([系数列表]) 来数值求解多项式的根。

4.2 调色方案的艺术化处理

着色逻辑是艺术创作的画笔。

  1. 迭代次数映射 ：上面我们用 V = 1.0 - (iter_count / max_iter) 进行线性映射。你可以尝试非线性映射来增强对比度，例如使用对数函数： 
    V = 1.0 - np.log1p(iter_count) / np.log1p(max_iter)
  2. 饱和度映射 ：除了明度，饱和度也可以变化。例如，让边界区域（迭代次数中等）饱和度更高，颜色更鲜艳；而快速收敛区（迭代次数少）和慢速收敛区（迭代次数多）饱和度降低，接近白色或灰色。 
    # 一个示例：让中间迭代次数的区域饱和度最高
normalized_iter = iter_count / max_iter
S = 1.0 - 4.0 * (normalized_iter - 0.5)**2  # 一个开口向下的抛物线
S = np.clip(S, 0.3, 1.0)
  3. 使用更丰富的色彩映射 ：不为每个根分配单一色调，而是分配一个小范围的色调区间，或者使用预定义的色彩映射（如 viridis , plasma ），再根据迭代次数进行插值，可以得到类似星云或地质图的效果。

4.3 性能优化与大规模渲染

当你想渲染4K甚至更高分辨率，或者探索极精细的细节时，性能成为瓶颈。

  1. 使用 numba 加速 ：如果算法逻辑复杂，难以完全向量化，可以使用 numba 的 @jit 装饰器编译单个像素的计算函数，然后结合多进程并行处理图像块。 
    from numba import jit, prange
@jit(nopython=True, parallel=True)
def compute_pixel(x, y, max_iter, tolerance, roots):
    # ... 单个像素的牛顿迭代和根匹配逻辑 ...
    return root_idx, iter_cnt
# 然后在循环中调用，或使用 prange 进行并行循环
  2. 自适应采样 ：分形图像大部分区域是平滑的，只有边界处极其复杂。可以采用自适应采样算法：先低分辨率渲染，检测颜色梯度大的区域（边界），然后只对这些区域进行高分辨率重采样，能极大节省计算量。
  3. GPU 计算 ：利用 cupy (对应 NVIDIA GPU) 或 jax 库，将整个网格计算移植到GPU上，可以获得数百倍的加速，实现实时缩放和探索。

5. 常见问题、调试技巧与避坑指南

在实际编码和渲染过程中，你肯定会遇到一些典型问题。以下是我踩过坑后总结的经验。

5.1 图像出现“斑点”或“条纹”伪影

  • 问题描述 ：生成的图像在某些区域有不自然的、块状或线状的错误颜色。
  • 可能原因与解决 ：
    1. 除以零 ：在牛顿迭代公式 z - f(z)/f'(z) 中，如果 f'(z) 为零，会导致计算溢出。虽然 z^3-1 的导数 3z^2 只在 z=0 为零，但数值计算中接近零也会产生极大值。
      • 解决 ：给分母加上一个极小值，如 dz + 1e-15 。
    2. 迭代溢出 ：某些初始点可能不收敛， z 值会变得非常大（趋于无穷），导致后续计算溢出。
      • 解决 ：在迭代循环中增加一个判断，如果 abs(z) 超过一个很大的阈值（如 1e10 ），则终止迭代，并将其标记为“不收敛”，可以赋予一个特殊的颜色（如黑色）。
    3. 根匹配容差问题 ：迭代结束后， z 值可能并未精确收敛到任何一个根，导致距离判断出错。
      • 解决 ：在计算到根的距离前，可以尝试对 z 进行一步“净化”操作，例如将其设置为距离最近的根（如果距离小于某个阈值）。或者，使用更宽松的匹配策略。

5.2 边界模糊，缺乏分形的精细细节

  • 问题描述 ：分形边界看起来糊成一团，没有那种“无限放大都有结构”的锐利感。
  • 可能原因与解决 ：
    1. 迭代次数不足 ： max_iter 设置得太低。在边界区域，迭代需要很多步才能满足收敛条件。增加 max_iter （如到100或200）可以捕获更精细的收敛过程。
    2. 收敛容差 tolerance 设置不当 ： tolerance 太大，迭代过早停止，导致很多点被错误地归类。通常 1e-6 是个不错的起点，对于高精度渲染，可以尝试 1e-8 或更小。
    3. 色彩映射过于平淡 ：仅用根来着色，没有充分利用迭代次数信息。 务必引入迭代次数来调制明度或饱和度 ，这是凸显边界细节的最有效方法。尝试不同的映射曲线（线性、对数、指数），找到最能突出细节的那一个。

5.3 计算速度太慢

  • 问题描述 ：渲染一张大图需要几分钟甚至更久。
  • 优化策略 ：
    1. 确保使用向量化 ：检查你的代码，确保所有对 f(z) 、 df(z) 和 z 的更新操作都是基于整个 NumPy 数组进行的，绝对避免对每个像素点写 Python 层的 for 循环。
    2. 降低初始分辨率 ：在调试和探索阶段，先将 width 和 height 设为较小的值（如 800x800）。
    3. 使用 float32 而非 float64 ：对于可视化，单精度浮点数通常足够，且能提升速度和减少内存占用。在创建数组时指定 dtype=np.complex64 。
    4. 裁剪感兴趣区域 ：不需要总是渲染 [-2,2]x[-2,2] 的全景。找到图案有趣的区域（如某个边界），缩小 x_min, x_max, y_min, y_max 的范围进行局部高清渲染。

5.4 探索特定区域：实现放大与漫游

最迷人的体验莫过于像使用地图软件一样，在分形世界中放大漫游。这需要将视图参数动态化。 

# 定义一个函数，给定视图中心和范围，生成分形
def generate_newton_fractal(center_re, center_im, range_re, range_im, width=800, height=800):
    x_min = center_re - range_re/2
    x_max = center_re + range_re/2
    y_min = center_im - range_im/2
    y_max = center_im + range_im/2
    # ... 后续使用这些边界生成网格和计算 ...
    return RGB_image

# 示例：放大查看 (0.3, 0.6) 点附近的有趣结构
img = generate_newton_fractal(center_re=0.3, center_im=0.6, range_re=0.5, range_im=0.5, width=1600, height=1600)

你可以将这个过程与一个简单的 GUI（如 matplotlib 的交互式点击事件）结合，实现点击哪里就放大哪里的效果，亲自感受“无限细节”的分形特性。

通过这个项目，你收获的远不止一幅漂亮的图片。你亲身体验了确定性算法中蕴含的混沌之美，实践了将数学公式转化为可视化代码的完整流程，并掌握了处理复杂计算和图像生成的实用技巧。牛顿法分形是一个完美的起点，从这里出发，你可以探索其他迭代函数系统（如曼德博集、朱利亚集），或者深入研究更复杂的动力系统。每一次参数的调整，每一次函数的改变，都可能打开一扇通往全新视觉宇宙的大门。

Logo
码道开发者社区

码道开发者社区，聚焦华为云码道 CodeArts 代码智能体，沉淀 Agent、Skill、鸿蒙开发实战内容，供开发者查阅资料、交流技术、分享工程实践

更多推荐

cover

不想先学一堆 React？用 Codex 帮你完成第一支 Remotion 动效视频

avatar 码道开发者社区
cover

突发！英伟达限制员工使用Claude；程序员一条代码删掉89TB数据获刑：被判五年十个月；曝iPhone Duo有望支持微信双开 | 极客头条

avatar 码道开发者社区
cover

每日热评｜学术牛马的救星还是幻觉放大器？47k星学术研究Agent技能包深度评测

avatar 码道开发者社区