说实话,南邮《数学实验》这门课,前几周大家还能对着PPT敲敲代码应付过去,到了模块二“函数的迭代”基本就开始两极分化了。有的同学照抄了老师给的示例,跑出图来却看不懂在干什么;有的同学代码报错连迭代是收敛还是发散都判断不了;还有的卡在期末报告那道“分析迭代收敛性与参数关系”的题上。借这篇文章,我把当时整理的一份完整参考思路和配套代码全部放出来,从原理、代码到测量结果,一条线讲透。
1. 这个模块到底在做什么
1.1 先看懂“迭代”这个动作
“函数的迭代”听起来是个很高大上的概念,其实说白了就是你手里有一个函数 (f(x)),任意给定一个初始值 (x_0),然后反复执行同一个操作:(x_1 = f(x_0)),(x_2 = f(x_1)),一直到 (x_{n+1} = f(x_n))。这样产生的序列
(x_0, x_1, x_2, \dots, x_n, \dots)
就是我们说的迭代序列。
很多同学一开始觉得这有什么好实验的?不就是套公式循环算吗?真正到了实验课你才会发现,同一个函数、不同参数、不同初值,跑出来的序列可能完全不同。有的序列会稳定在一个数附近,有的会在几个数之间来回跳,有的看起来完全没有规律。这门实验的核心任务,就是通过数值实验观察、分类、解释这些“奇怪”的现象。
1.2 这个模块在整门课里的定位
如果你把《数学实验》这门课看作一次“从公式到现象”的思维训练,那模块二就是第一个真正让你脱离解析推导的限制,靠计算机去“看”数学行为的模块。前面的线性代数实验大多还能手算验证,到了函数迭代这里,很多现象你根本没法用笔算出来,比如混沌状态下的迭代序列,解析上只能证明它有界,但具体怎么跑,必须上机器。
我在做这个实验时最大的感受是:这门课的重点并非“编程”,而是“用程序做数学观察”。老师最终想考察的,是你能否通过迭代观察到不动点、周期点、混沌等动力学行为,并且把观察结果用图形和数值清晰地表达出来。所以,参考答案给代码只是第一步,更重要的是理解每张图、每个数值到底对应哪个数学概念。
1.3 最终要交出什么样的报告
南邮这个实验模块的作业通常要求包含:实验目的、迭代原理说明、程序代码、运行结果图、结果分析、思考题回答。其中“结果分析”是最拉分的部分。同一个班很多人的代码是从学长那边拷贝的,图都差不多,但你能不能在几幅图里准确说出“这是一个2周期点”“分岔点大约在 r=3.449 附近”“Feigenbaum常数测量值约为4.66”,这决定了报告是60分还是90分。下面的章节我直接按这个要求来拆。
2. 主要实验内容与核心函数设计
2.1 实验一:用二次函数观察收敛与发散
第一个必做实验通常是考察函数
(f(x) = x^2 + c)
在不同参数 (c) 下的迭代行为。实际做的时候,你会调整 (c) 的取值,观察序列 (x_n) 的长期行为。
我当时在MATLAB里写了这样一个简单的脚本:
% 函数迭代实验:f(x) = x^2 + c clear; clf; x0 = 0.2; % 初值 N = 100; % 迭代步数 c_list = [0.5, -0.5, -1, -1.3, -2]; % 不同参数 for i = 1:length(c_list) c = c_list(i); x = zeros(1, N); x(1) = x0; for n = 1:N-1 x(n+1) = x(n)^2 + c; end subplot(2, 3, i); plot(1:N, x, '.-', 'MarkerSize', 6); title(['c = ', num2str(c)]); xlabel('迭代步数 n'); ylabel('x_n'); grid on; end跑出来的现象大致是这样的:
| 参数 c | 表现 | 说明 |
|---|---|---|
| 0.5 | 序列迅速增大,趋于无穷 | 发散 |
| -0.5 | 序列从0.2开始衰减,最后稳定在约0.366 | 收敛到不动点 |
| -1 | 序列在0和-1之间反复跳动 | 2周期点 |
| -1.3 | 序列在4个值之间循环 | 4周期点 |
| -2 | 序列长期在[-2,2]内波动,但似乎不重复 | 混沌或高周期 |
这个实验的核心收获是:一个形式上极其简单的二次函数,仅仅改变一个常数,就能产生如此丰富的动态行为。很多人跑完图没有感觉,我建议你额外做一件事:把 (c=-1) 情况和 (c=-2) 情况打印出后面20步的具体数值,瞪大眼睛比较“周期循环”和“看似乱跳”的区别,这个对比比图上的线条更有冲击力。
2.2 实验二:蛛网图可视化迭代过程
第二个实验通常是画蛛网图,也叫科布韦布图。它的画法非常直观:在坐标平面里画出 (y=f(x)) 曲线和 (y=x) 对角线,然后从 x 轴上的初始点开始,做竖线到曲线上得到 f(x_0),再做横线到对角线上,得到下一个 x_1,再重复竖线、横线。
我当时写了一个通用函数来画蛛网图,这样后面多个实验都能复用:
function cobweb_plot(f, x0, n, xlim_range) % f: 函数句柄 % x0: 初始值 % n: 迭代次数 % xlim_range: x显示范围 [xmin, xmax] t = linspace(xlim_range(1), xlim_range(2), 400); plot(t, f(t), 'b-', 'LineWidth', 1.5); hold on; plot(t, t, 'k--', 'LineWidth', 1.0); x = x0; for i = 1:n y = f(x); plot([x, x], [x, y], 'r-', 'LineWidth', 0.8); plot([x, y], [y, y], 'r-', 'LineWidth', 0.8); x = y; end xlabel('x_n'); ylabel('x_{n+1}'); grid on; end调用方式比如:
f = @(x) x.^2 - 1; cobweb_plot(f, 0.2, 50, [-1.5, 1.5]);蛛网图的真正价值不是好看,而是把“序列收敛/周期/混沌”变成了几何信息。收敛时蛛网线绕成一个“螺旋”粘在曲线与对角线的交点上;周期时蛛网会形成一个闭合的矩形回路;混沌时蛛网在某一区域内杂乱缠绕,把整个区域越画越满。
2.3 实验三:Logistic 映射与倍周期分岔图
如果你只做前两个实验,期末肯定不够。模块二的重头戏,是研究经典的Logistic映射:
(x_{n+1} = r , x_n (1 - x_n))
这个模型本来是生态学里用来描述种群数量变化的:r代表增长率,x代表当前种群数量占环境容纳量的比例。但它后来成为混沌理论的教科书级例子,原因就是参数 r 从 2 增加到 4 的过程中,系统的最终行为呈现出极其规整的从周期走向混沌的路径。
绘制分岔图的程序,我当时是这样写的:
% Logistic映射分岔图 clear; clf; r_list = 2.5:0.005:4.0; % 参数范围 x0 = 0.3; % 统一初值 N_transient = 200; % 丢弃的前瞬态点数 N_plot = 200; % 保留的画图点数 for i = 1:length(r_list) r = r_list(i); x = x0; for n = 1:N_transient x = r * x * (1 - x); end for n = 1:N_plot x = r * x * (1 - x); plot(r, x, '.', 'MarkerSize', 1, 'Color', [0, 0.2, 0.6]); hold on; end end xlabel('参数 r'); ylabel('迭代最终状态 x'); title('Logistic映射分岔图');这张图一出来,几乎每个第一次跑出分岔图的同学都会惊叹:在 r < 3 时只有一个值,所以图上是一条线;在 r 超过 3 之后裂成两支,这就是2周期;后续不断裂变,形成4、8、16……直到某个临界值之后突然出现一片密密麻麻的点,那就是混沌区。更神奇的是,混沌区里还偶尔会突然出现几条干净的“白色竖线”,那是周期窗口,比如 r 约等于 3.83 附近有一个明显的3周期窗口。
2.4 实验四:测量倍周期分岔点与 Feigenbaum 常数
这个模块的高级题目通常会要求你定量分析分岔点。你需要找到四个连续分岔点的位置:(r_1, r_2, r_3, r_4),然后用公式
(\delta = \frac{r_3 - r_2}{r_4 - r_3})
估算 Feigenbaum 常数,理论上它的值是 4.669201609...。
因为直接肉眼从分岔图读分岔点误差太大,我当时采用了一个比较聪明的扫描方法。原理是:在周期区域内,长时期迭代后的序列值集合是有限个点;从2周期变4周期时,序列值的种类会翻倍。于是我可以对每个 r,迭代足够多步,记录最后若干个不同数值的个数,然后观察个数突变的位置。
这里有一个很小的技巧:判断两个数值是否“不同”不要用严格相等,要用容差判断,否则因为计算机浮点误差,哪怕理论上是同一个点,实际算出两个相差 1e-15 的数也会被误算成两个周期点。
我当时用的是一段粗扫加细扫结合的代码:
function r_bif = find_bifurcation(r_start, r_end, period_target) % 找到从指定周期翻倍的参数近似位置 tol = 1e-6; options = optimset('TolX', 1e-10); f = @(r) period_count(r) - period_target; r_bif = fzero(f, [r_start, r_end], options); end function num = period_count(r) x = 0.3; N = 500; % 迭代次数 for n = 1:N x = r * x * (1 - x); end vals = zeros(1, 512); vals(1) = x; for n = 2:512 x = r * x * (1 - x); vals(n) = x; end num = length(unique(round(vals, 6))); % 四舍五入到1e-6去重 end这套代码在精度要求不夸张的情况下,实际测出来的分岔点大约是:
| 分岔 | 参数位置 |
|---|---|
| (r_1)(1周期→2周期) | 3.0000 |
| (r_2)(2周期→4周期) | 3.4495 |
| (r_3)(4周期→8周期) | 3.5441 |
| (r_4)(8周期→16周期) | 3.5644 |
由此计算出:
(\delta = \frac{3.5644 - 3.5441}{3.5441 - 3.4495} \approx \frac{0.0203}{0.0946} \approx 4.66)
与理论的4.6692相差不到1%,实验报告写到这个程度就已经相当有说服力了。
3. 实操过程、参数选型与常见坑
3.1 实验步骤的合理推进顺序
我建议你实际动手时不要直接一股脑跑分岔图,而是按下面的顺序一步步来,效率最高也最不容易卡壳。
第一步,先用最传统的迭代循环脚本跑出实验一的时间序列图,把收敛、2周期、4周期这几张图印在脑子里。第二步,调用蛛网图函数,针对同样的参数画几张蛛网图,进一步建立几何直觉。第三步,进阶到 Logistic 映射,先用固定 r 看序列,比如 r=2.8 应该收敛,r=3.2 应该2周期,r=3.5 应该4周期,r=3.9 应该混沌。第四步再画完整分岔图,第五步做分岔点测量。这样每走一步都能验证前面的认知,而不是等画出一张复杂的分岔图后一脸茫然。
3.2 参数选型的“为什么”
有几个关键参数我单独说明一下为什么这么设置。
迭代瞬态丢弃数 N_transient=200:因为从任意初值出发的迭代序列,前一部分数据是受初值影响较大的“暂态过程”,最终会趋于某种长期行为(也就是吸引子)。我们画分岔图关心的是“长期行为”,所以前面那些点必须丢掉。200这个数对于 Logistic 映射来说绝大多数情况下足够,但如果你想精确观察 r 非常接近分岔点的行为,建议加到500。
画图密度 r_list = 2.5:0.005:4.0:0.005是较为常用的分辨率,肉眼观察足够,但如果你想把图放大找精细结构,这个分辨率会显得稀疏。我建议,如果计算机内存允许,可以先用0.005得到全貌,然后对特定区间如[3.8, 3.9]再用0.001去补一个局部放大图,这样的两图组合在你的报告中非常加分。
3.3 实际操作中最常遇到的坑
第一个大坑是MATLAB里面忘记用点运算。很多同学定义了 (f(x)=x^2+c),然后写成 x.^2 + c 没问题,但在画蛛网图的时候如果传入的是一个向量 t,写了 t^2 + c,就会直接报错“矩阵维度不一致”。我建议自定义函数时一律养成用.^、.*、./的习惯,这样函数自动支持向量输入,画图调用会方便很多。
第二个坑是分岔图细节不够,出图后感觉一团糊。这种情况大多是 r_list 的步长太粗,或者画点太多导致图片文件过大。解决方法是分区间绘制,不要指望一张图展示从2.5到4.0的全部细节,局部放大永远是更好的展示方式。
第三个坑比较隐蔽,就是迭代序列“周期数判断”失效。很多人用 uniquetol 去重时发现 r=3.2 时最终长期迭代竟然出现了6个不同值,怎么看都不对。原因在于迭代次数不够,序列尚未完全收敛到2周期点;尤其是 r 非常接近分岔点的时候,收敛速度特别慢,迭代300次都不一定稳定。这里我建议用两种方法交叉验证:一是增加迭代次数到2000以上,二是同时计算相邻两点的差值是否小于 1e-8,如果小于则认为已到达数值意义上的周期。
3.4 牛顿迭代法在模块二里的出现
有些实训版本在模块二里还附带了一个牛顿迭代法的应用实验,用来求解非线性方程,比如求 (x^3 - x - 1 = 0) 的根。虽然牛顿法的式子是
(x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)})
但它本质上同样是一种函数迭代。这个实验如果出现,核心目的应该是强调迭代收敛的条件:在解附近必须有较好初值,否则可能会发散或收敛到另一个根。我当时在这个问题上吃过亏,用 MATLAB 求解时给了一个远离真根的初值,结果迭代20步后跑到一个无意义的极大值去了。后面我总结出一个技巧:先用 fplot 画一下函数曲线,从图上目测零点大概位置,再用这个位置作初值,百试百灵。
4. 结果分析与常见问题排查
4.1 如何对图形结果做分析
报告里最容易写空的地方就是“结果分析”。很多同学写“从图中可以看出,当c=-1时,迭代序列在两个点之间振荡”,然后就没有下文了。这只能算描述,不叫分析。
我更推荐用这种“描述现象→给出数值证据→对应数学概念”的三段式结构。以 c=-1 为例,你可以这样写:观察图1(c),迭代100步后序列在 x≈0 和 x≈-1 之间交替出现,进一步打印第50到60步数值可知 xn 近似为-1、0交替,这说明系统进入了一个稳定的2周期轨道;从几何上,蛛网图(图2)呈现一个闭合矩形回路,表明迭代在两个点之间循环往返;这一现象与非线性映射的倍周期分岔理论相符,也为后续在Logistic映射中观察周期倍增提供了基础。
这样写,老师一眼就能看出你既读懂了图,也读懂了背后的理论。
4.2 常见异常结果速查表
| 现象 | 可能原因 | 排查或修正方法 |
|---|---|---|
| 迭代序列出现NaN或Inf | 参数选取过大导致迭代值溢出 | 限制参数范围,初值尽量取在[0,1]内,程序里加边界判断 |
| 蛛网图只画出一个闪烁的点 | 迭代次数n设置太小 | 增加迭代次数到50以上,曲线会逐步缠绕成型 |
| 分岔图在混沌区有间歇性空白 | 瞬态点数不够 | 将N_transient提升到500或1000 |
| 分岔图左侧只有一条线,看不到分岔 | r范围起点太小或初值x0选为0 | Logistic映射中 x0=0 会永远停在0,换用0.2~0.4的初值 |
| 测量Feigenbaum常数偏差很大 | 分岔点定位不准确 | 用fzero做局部精扫,不要直接看像素读数 |
| 牛顿迭代不收敛 | 初值太远或函数导数接近0 | 先用画图法预估根的位置,或改用阻尼牛顿法 |
4.3 画图细节:让出图更专业
出了正确的结果,图好不好看、清不清晰,也会影响分数。我当时总结了一套画图规范,你可以直接抄作业。
分岔图里,点的大小不要超过2,颜色尽量用深色系;多条曲线在同一张图里时,要用不同颜色加图例,但注意不要花里胡哨。线程图中,收敛序列建议用带圆点的实线,发散序列用虚线,能让人一眼分辨。蛛网图的曲线和迭代折线颜色要有明显区分,比如蓝色曲线加红色折线。所有图必须加 xlabel、ylabel,有多个子图时需要统一标题格式。
还有一个小细节,输出图的时候不要截屏再贴到Word里,不然清晰度很差。我建议用MATLAB的 exportgraphics 或者 saveas,把图存成 300dpi 的 PNG 再插入文档,字迹和线条都清楚很多。
4.4 几个值得做的进阶验证
如果你想在报告中体现一些超出课件的思考,我强烈建议试一下以下两个验证实验,代码量很小但分析空间很大。
第一个是验证“初值敏感性”。取 r=3.9,从 x0=0.300 和 x0=0.301 分别迭代50步,在第1到10步两条序列几乎重合,第20步左右开始明显分离,到第50步已经完全不相关。这个对比图几乎就是“蝴蝶效应”的数值证明,写进报告会让你的结果分析立刻高一个层次。
第二个是画出“周期窗口”的局部放大图。将 r 范围取[3.82, 3.86],同样画分岔图,你会看到经典的3周期窗口,并且在窗口内部又出现3周期态向6周期态的倍周期分岔。这说明混沌区内嵌套着精细的自相似结构,也直接呼应了Feigenbaum常数的普适性。
5. 思考题与报告加分项
5.1 为什么收敛条件要看 ( |f'(x^*)| < 1 )
很多思考题会问“为什么不动点的稳定性取决于导数绝对值是否小于1”。这个问题用一阶泰勒展开就能解释。
假设 (x^) 是不动点,小扰动 (e_n = x_n - x^),那么:
(e_{n+1} = x_{n+1} - x^* = f(x_n) - f(x^) \approx f'(x^) (x_n - x^) = f'(x^) e_n)
所以每迭代一次,扰动大约被乘以 (f'(x^*))。这个值绝对值小于1,扰动会逐次缩小,迭代收敛;大于1,扰动被不断放大,即使初始值离不动点非常近,也会被推离。我在报告中建议不只是写公式,还要画一个小图,用 (f(x)=x^2+c) 在不动点处的切线斜率来说明,对比 c=-0.5 和 c=2 两种情况的斜率大小,直观得多。
5.2 周期点与混沌之间是什么关系
这个问题在思考题中经常换着法子出现。本质上是说:周期点是迭代后经过有限步回到自身的点,而混沌轨道的显著特征之一是“永不精确重复但长期有界”。
从实验数据上,你可以这样区分:打印 c=-1 下序列的后20项,值只会是0或-1;打印 c=-2 下序列的后20项,每个数都不同,但都在[-2,2]里。再进一步,可以计算相邻项差分的绝对值,周期情况下差分是固定模式重复,混沌情况下差分的变化毫无规律。如果能结合自相关或者功率谱分析,那报告的专业程度就能拔高很多,不过这些通常超出课程要求,作为选做即可。
5.3 报告里适合额外写的“实验感悟”
如果你想在结尾部分写一点感悟,我建议千万不要写“通过本次实验,我掌握了……”这种套话。更好的方式是写一个非常具体的小观察,比如:“在测量Feigenbaum常数时,我原本以为需要极高的计算精度才能得到接近4.669的结果,但实际在分岔点测量误差达到0.001左右时,常数就已经稳定在4.66附近。这让我体会到倍周期分岔路径的内在普适性并非一个只能通过复杂计算验证的抽象结论,而是一个可以被简单数值实验直观感受的性质。”
这种具体化的感性描述,通常比空泛的总结更受老师认可。
6. 关于工具选型与跨平台实现
6.1 MATLAB / Octave / Python 怎么选
南邮这门课的正规环境是 MATLAB,但有的同学电脑装不上或序列号过期,用 GNU Octave 也能跑绝大多数脚本,因为基本语法兼容。你只需要注意Octave的图形窗口交互不如MATLAB顺畅,分岔图点多了之后拖动会卡,但这不影响出图。
如果你本来就更熟悉 Python,用 NumPy 和 Matplotlib 完全可以把所有实验复现一遍。核心逻辑没有任何区别,Python 的 matplotlib 甚至在做局部放大图时更容易控制坐标轴范围。我个人的建议是:如果目标是课程拿高分,先老老实实按MATLAB的作业要求来,Python 可以作为验证手段;如果以后打算走数据分析或算法方向,多用 Python 熟悉一下也不是坏事。
6.2 一段可复用的Python参考代码
这里给一段对应Logistic分岔图的Python代码,方便想做交叉验证的同学参考:
import numpy as np import matplotlib.pyplot as plt r_values = np.linspace(2.5, 4.0, 2000) x0 = 0.3 N_transient = 500 N_plot = 200 plt.figure(figsize=(10, 6)) for r in r_values: x = x0 for _ in range(N_transient): x = r * x * (1 - x) for _ in range(N_plot): x = r * x * (1 - x) plt.plot(r, x, ',', color='navy', markersize=0.5) plt.xlabel('r') plt.ylabel('x') plt.title('Bifurcation diagram of Logistic map') plt.xlim(2.5, 4.0) plt.show()这段代码跑出来的图和MATLAB版本几乎没有差别。要注意的是 r_values 的密度,Python 这边我用2000个点,画图阶段每点只画最后一次迭代值,速度非常快,但如果要画每个r对应200个点,总点数就到了40万以上,matplotlib可能有点慢。这种情况下可以把图形后端换成 Agg,或者用散点图一次性传数组,而不是循环里一个个 plot。
6.3 运行环境与报错处理
MATLAB新版和旧版在 plot 颜色简写和 hold on 行为上略有差异。如果你用的是R2019b之前的版本,某些地方可能不会自动加 hold 开关,建议显式在每个子图循环开头写上 hold on。Octave 对匿名函数 @(x) x.^2+c 支持良好,但 older 版本中 fzero 求解行为可能不同,如果报错,改用 fsolve 或自己写二分法也能达到同样目的。
Python 端最容易遇到的问题是中文字体显示成方框。在绘图前加上:
plt.rcParams['font.sans-serif'] = ['SimHei'] plt.rcParams['axes.unicode_minus'] = False否则标题 “分岔图” 会变成一堆乱码。这个不算什么高深技术,但每年都有人在这里卡半天。
7. 从做实验到真正理解迭代
做完这套实验我最大的感触就是:数学里很多概念,光靠黑板推导很难在脑子里“活”起来,但当你看到那条 Logistic 映射分岔图一点点从一根线长成两支、四支、八支,最后变成一片混沌区域的时候,你才真正理解为什么物理学家费根鲍姆发现那个常数时会那么兴奋——原来这背后有一种跨越具体系统的普适规律,而你的电脑屏幕上就能复现它。
如果你只是照着参考答案把图跑出来交了作业,说实话,你亏了。我建议你在交完报告之后,自己再花一个小时,把参数 r 慢慢从2.5推到4.0,一次只改一点点,盯着那个序列值的变化,那种从有序走向混沌的渐变感,比任何文字描述都有说服力。这也是这门实验藏在作业背后的真正目的。