numpy微实践 - 分形几何森林中的一片树叶

海洋饼干叔叔

发布于：2024年3月25日

Python

曲折的海岸线、起伏不平的山峦、大地上的河流、纵横交错的人体血管、树叶的脉络、甚至一棵花椰菜都展现出传统几何学难以描述的形状：其结构精细而不规则，整体和局部呈现自相似特性。数学上，称这种整体和局部呈自相似的图形特性为分形（fractal）。为了研究与分形相关的数学规律，人们甚至专门创立了一个新的数学分支 —— 分形几何学（fractal geometry）。

这里，我们撷取分形几何森林中的一片树叶，与读者共同领略数学之美。

版权声明

本文可以在互联网上自由转载，但必须：注明出处(作者：海洋饼干叔叔)并包含指向本页面的链接。

本文不可以以纸质出版为目的进行改编、摘抄。

本文节选自《Python编程基础及应用》，高等教育出版社，作者：陈波，刘慧君

表1 线性函数组

图1 分形蕨类树叶

迭代函数系统（iterated function system）是一种创建分形图案的简单算法。表1列出了的4个线性函数组，每个函数组均可以根据一个二维平面点的点坐标（x_i,y_i）计算得到一个新的点坐标（x_i+1,y_i+1）。从坐标原点(0,0)出发，反复利用上述线函数组进行迭代计算，即可凭“空”生成如图1所示的分形树叶。具体步骤如下：

① 将坐标原点(0,0)赋值给(x0,y0)；

② 按表10-x规定的选用概率随机选择一个线性函数组；

③ 应用该函数组对(x0,y0)进行计算，得到一个新点坐标(x1,y1)；

④ 重复②③步，从(x1,y1)计算得到(x2,y2)，从(x2,y2)计算得到(x3,y3)…… 总共计算得到10万个坐标点；

⑤ 将迭代计算而得的10万个点绘制成散点图，即得图1。

仔细观察图1，可见，如果切下这片树叶中的任何一片子叶，竖立起来并放大，会发现它跟整片树叶长得一模一样。这说明图案的整体跟局部是自相似的，这是一个分形图案。

在编程实现上述迭代函数系统之前，需要解决两个具体的技术问题。其一是如何以规定的选用概率在迭代过程中随机选择函数组。以表1中的函数组1为例，其选用概率为0.07，这意味着，10万个点中的7%，需要使用函数组1来计算生成。借助于numpy的searchsorted()函数以及随机数数组，可以便捷地解决这一问题。

IDLE Shell
>>>import numpy as np
>>>p = np.cumsum([0.01,0.07,0.07,0.85]) #对选用概率数组进行累积求和
>>>p
array([0.01, 0.08, 0.15, 1.  ])

首先，应用np.cumsum()函数对各函数组的选用概率“数组”进行累积求和。观察结果数组p，容易发现，p[2] = 0.15 = 0.01 + 0.07 + 0.07，它正好等于原始概率数组中下标0至下标2的3个元素之和。如果把结果数组p标记在数轴上，如图2所示，0与p[0]之间的“缝隙”宽度为0.01，p[0] ~ p[1]之间的“缝隙”宽度为0.07，p[2] ~ p[3]之间的“缝隙”宽度为0.85，正好等于相关函数组的选用概率。

图2 数轴上的概率累积数组p

IDLE Shell
>>>rands = np.random.rand(100000)
>>>rands
array([0.83528652, 0.53627275, 0.34296887, ..., 0.79190539, 0.49528487,
       0.55354825])

接下来，生成包含10万个元素的均匀分布随机数数组rands，随机数的值域为[0,1)。不难想象，如果把介于0和1之间的随机数“倾泻”在图2所示的数轴上，85%的数会落在0.15和1.0之间，1%的数会落在0和0.01之间…… 这正是下述np.searchsorted()函数的所做的工作。

IDLE Shell
>>>selects = np.searchsorted(p,rands)
>>>selects[:200]
array([3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 0, 3, 3, 3, 3, 3, 3, 3,
       ... #有节略
       3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 2, 3, 1, 3, 3,
       ... #有节略
       3, 3], dtype=int64)

np.searchsorted(p,rands)逐一取出rands数组中的随机数r，然后在数组p中进行搜索，找出使得r≤p[i]成立的最小下标i，并将其生成在结果数组selects中。以随机数rands[0] = 0.83528652为例，p[2] < rands[0] ≤ p[3]，故selects[0] = 3。

由searchsorted()函数生成的selects数组包括10万个整数，由于rands中的随机数在[0,1)区间内均匀分布，这些整数中，包含约1%的0、7%的1、7%的2以及85%的3。在后续程序中，selects[i]（取值0~3）即对应第i轮迭代时应选用的函数组的编号。

程序需要处理的另一个技术问题是应用函数组从点坐标(x_i,y_i)计算得到下一个点坐标(x_i+1,y_i+1)。线性代数课程中讨论的矩阵乘法是解决这一问题的有力工具。以表1中的函数组2为例，应用该函数组从(x_i,y_i)计算得到(x_i+1,y_i+1)的过程可用矩阵乘法实现。

IDLE Shell
>>>import numpy as np
>>>eq = np.array([[-0.14,  0.26, 0], [ 0.26, 0.25, 0.47]])
>>>eq                              #函数组的系数矩阵，2行3列
array([[-0.14,  0.26,  0.  ],
       [ 0.26,  0.25,  0.47]])
>>>pt = np.array([-0.8,6.7,1])
>>>pt                              #点坐标向量[xi,yi,1]，1行3列
array([-0.8,  6.7,  1. ])
>>>np.dot(eq,pt)                  #求两个矩阵的积
array([1.854, 1.937])

按照矩阵运算的数学规则，2行3列的函数组系数矩阵应与“3行1列”的点坐标矩阵pt相乘，结果为2行1列的新点坐标矩阵。但在上述IDLE交互过程中，数组pt事实上是包含3个元素的一维数组（向量），np.dot()函数在将矩阵eq和向量pt相乘前，会自动将向量pt转置成“3行1列”。

请读者注意，向量pt除了包含平面点的x和y坐标外，还包含了一个值为1的浮点数，以便处理函数组中的常数。

本实践的完整程序如下：

#ifsleaf.py
import numpy as np, matplotlib.pyplot as plt
eqs = [np.array([[ 0,     0,    0], [ 0,    0.15, 0   ]]),
       np.array([[ 0.21, -0.19, 0], [ 0.24, 0.27, 1.5 ]]),
       np.array([[-0.14,  0.26, 0], [ 0.26, 0.25, 0.47]]),
       np.array([[ 0.87,  0,    0], [-0.05, 0.84, 1.54]])]
p = np.cumsum([0.01,0.07,0.07,0.85])
N = 100000
rands = np.random.rand(N)
selects = np.searchsorted(p,rands)

pt = np.array([0,0,1],dtype=np.float64)
pts = np.zeros((N,2),dtype=np.float64)
for i in range(N):
    r = np.dot(eqs[selects[i]],pt)
    pts[i] = pt[:2] = r

plt.figure(figsize=(3,5))
plt.scatter(pts[:,0],pts[:,1],s=1,c='g',linewidths=0)
plt.axis("off")
plt.subplots_adjust(0,0,1,1)

plt.figure(figsize=(3,5))
plt.scatter(pts[:,0],pts[:,1],s=1,c=selects,linewidths=0)
plt.axis("off")
plt.subplots_adjust(0,0,1,1)

plt.show()

🚩第3 ~ 6行：eqs列表包含了4个2行3列的数组，对应4个函数组的系数矩阵。

🚩第7 ~ 10行：依稍前讨论的方法，生成包含10万个整数的函数组随机选用数组selects。其中，selects[i]对应生成第i个平面点时应使用的函数组编号。

🚩第12行：当前平面点坐标向量。整片树叶的迭代生成过程从坐标原点(0,0)出发。如前所述，pt向量尾部附加的1用于处理函数组中的常数。

🚩第13行：生成形状为(100000,2)的二维零值数组pts，用于保存全部平面点的坐标。第i轮迭代生成的平面点坐标保存于pts[i]。

🚩第14 ~ 16行：循环N(10万)次，迭代生成10万个平面点。

🚩第15行：以selects[i]为下标，选取对应的函数组系数矩阵，与当前平面点坐标向量pt相乘，得到新平面点坐标向量r。

🚩第16行：将新平面点坐标向量r存入平面点坐标数组pts；更新当前平面点坐标向量pt的坐标值，保持末尾元素1不变，作为下一轮迭代计算的输入。

🚩第18 ~ 21行：以绿色(c=’g’)绘制散点图至Figure 1，见图1。pts[:,0]，pts[:,1]是对二维数组的多维切片，分别对应全部平面点的x坐标和y坐标。关于多维切片语言的详细解释，请回顾本书10.1.3节。

🚩第23 ~ 26行：以selects向量为颜色绘制散点图至Figure 2。第24行中的参数c = selects表示在绘制第i个平面点时，其颜色由selects[i]的值决定。如图3所示，由不同函数组生成的平面点以不同颜色作了区分。选用概率为1%的函数组0生成了深红色的叶杆部分，选用概率为7%的函数组1、2生成了左右两片子叶（绿、蓝），选用概率为85%的函数组3则生成了树叶主体的黄色部分。

请读者观察函数组0，无论当前平面点的坐标(x_i,y_i)为何值，函数组0计算而得的下一个平面点的横坐标始终为0。这解释了图3中深红色叶杆部分为何左右居中（横坐标为0）且竖直向下。

本实践的代码并不难解读，但要真正理解IFS的迭代方程组是如何生成这片分形树叶却并不容易。这需要更多的线性代数以及二维仿射变换的知识，这超出了本书的范围。