说实话南邮《数学实验》这门课前几周大家还能对着PPT敲敲代码应付过去到了模块二“函数的迭代”基本就开始两极分化了。有的同学照抄了老师给的示例跑出图来却看不懂在干什么有的同学代码报错连迭代是收敛还是发散都判断不了还有的卡在期末报告那道“分析迭代收敛性与参数关系”的题上。借这篇文章我把当时整理的一份完整参考思路和配套代码全部放出来从原理、代码到测量结果一条线讲透。1. 这个模块到底在做什么1.1 先看懂“迭代”这个动作“函数的迭代”听起来是个很高大上的概念其实说白了就是你手里有一个函数 (f(x))任意给定一个初始值 (x_0)然后反复执行同一个操作(x_1 f(x_0))(x_2 f(x_1))一直到 (x_{n1} f(x_n))。这样产生的序列(x_0, x_1, x_2, \dots, x_n, \dots)就是我们说的迭代序列。很多同学一开始觉得这有什么好实验的不就是套公式循环算吗真正到了实验课你才会发现同一个函数、不同参数、不同初值跑出来的序列可能完全不同。有的序列会稳定在一个数附近有的会在几个数之间来回跳有的看起来完全没有规律。这门实验的核心任务就是通过数值实验观察、分类、解释这些“奇怪”的现象。1.2 这个模块在整门课里的定位如果你把《数学实验》这门课看作一次“从公式到现象”的思维训练那模块二就是第一个真正让你脱离解析推导的限制靠计算机去“看”数学行为的模块。前面的线性代数实验大多还能手算验证到了函数迭代这里很多现象你根本没法用笔算出来比如混沌状态下的迭代序列解析上只能证明它有界但具体怎么跑必须上机器。我在做这个实验时最大的感受是这门课的重点并非“编程”而是“用程序做数学观察”。老师最终想考察的是你能否通过迭代观察到不动点、周期点、混沌等动力学行为并且把观察结果用图形和数值清晰地表达出来。所以参考答案给代码只是第一步更重要的是理解每张图、每个数值到底对应哪个数学概念。1.3 最终要交出什么样的报告南邮这个实验模块的作业通常要求包含实验目的、迭代原理说明、程序代码、运行结果图、结果分析、思考题回答。其中“结果分析”是最拉分的部分。同一个班很多人的代码是从学长那边拷贝的图都差不多但你能不能在几幅图里准确说出“这是一个2周期点”“分岔点大约在 r3.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(n1) 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 实验二蛛网图可视化迭代过程第二个实验通常是画蛛网图也叫科布韦布图。它的画法非常直观在坐标平面里画出 (yf(x)) 曲线和 (yx) 对角线然后从 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_{n1}); grid on; end调用方式比如f (x) x.^2 - 1; cobweb_plot(f, 0.2, 50, [-1.5, 1.5]);蛛网图的真正价值不是好看而是把“序列收敛/周期/混沌”变成了几何信息。收敛时蛛网线绕成一个“螺旋”粘在曲线与对角线的交点上周期时蛛网会形成一个闭合的矩形回路混沌时蛛网在某一区域内杂乱缠绕把整个区域越画越满。2.3 实验三Logistic 映射与倍周期分岔图如果你只做前两个实验期末肯定不够。模块二的重头戏是研究经典的Logistic映射(x_{n1} 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 看序列比如 r2.8 应该收敛r3.2 应该2周期r3.5 应该4周期r3.9 应该混沌。第四步再画完整分岔图第五步做分岔点测量。这样每走一步都能验证前面的认知而不是等画出一张复杂的分岔图后一脸茫然。3.2 参数选型的“为什么”有几个关键参数我单独说明一下为什么这么设置。迭代瞬态丢弃数 N_transient200因为从任意初值出发的迭代序列前一部分数据是受初值影响较大的“暂态过程”最终会趋于某种长期行为也就是吸引子。我们画分岔图关心的是“长期行为”所以前面那些点必须丢掉。200这个数对于 Logistic 映射来说绝大多数情况下足够但如果你想精确观察 r 非常接近分岔点的行为建议加到500。画图密度 r_list 2.5:0.005:4.00.005是较为常用的分辨率肉眼观察足够但如果你想把图放大找精细结构这个分辨率会显得稀疏。我建议如果计算机内存允许可以先用0.005得到全貌然后对特定区间如[3.8, 3.9]再用0.001去补一个局部放大图这样的两图组合在你的报告中非常加分。3.3 实际操作中最常遇到的坑第一个大坑是MATLAB里面忘记用点运算。很多同学定义了 (f(x)x^2c)然后写成 x.^2 c 没问题但在画蛛网图的时候如果传入的是一个向量 t写了 t^2 c就会直接报错“矩阵维度不一致”。我建议自定义函数时一律养成用.^、.*、./的习惯这样函数自动支持向量输入画图调用会方便很多。第二个坑是分岔图细节不够出图后感觉一团糊。这种情况大多是 r_list 的步长太粗或者画点太多导致图片文件过大。解决方法是分区间绘制不要指望一张图展示从2.5到4.0的全部细节局部放大永远是更好的展示方式。第三个坑比较隐蔽就是迭代序列“周期数判断”失效。很多人用 uniquetol 去重时发现 r3.2 时最终长期迭代竟然出现了6个不同值怎么看都不对。原因在于迭代次数不够序列尚未完全收敛到2周期点尤其是 r 非常接近分岔点的时候收敛速度特别慢迭代300次都不一定稳定。这里我建议用两种方法交叉验证一是增加迭代次数到2000以上二是同时计算相邻两点的差值是否小于 1e-8如果小于则认为已到达数值意义上的周期。3.4 牛顿迭代法在模块二里的出现有些实训版本在模块二里还附带了一个牛顿迭代法的应用实验用来求解非线性方程比如求 (x^3 - x - 1 0) 的根。虽然牛顿法的式子是(x_{n1} 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选为0Logistic映射中 x00 会永远停在0换用0.2~0.4的初值测量Feigenbaum常数偏差很大分岔点定位不准确用fzero做局部精扫不要直接看像素读数牛顿迭代不收敛初值太远或函数导数接近0先用画图法预估根的位置或改用阻尼牛顿法4.3 画图细节让出图更专业出了正确的结果图好不好看、清不清晰也会影响分数。我当时总结了一套画图规范你可以直接抄作业。分岔图里点的大小不要超过2颜色尽量用深色系多条曲线在同一张图里时要用不同颜色加图例但注意不要花里胡哨。线程图中收敛序列建议用带圆点的实线发散序列用虚线能让人一眼分辨。蛛网图的曲线和迭代折线颜色要有明显区分比如蓝色曲线加红色折线。所有图必须加 xlabel、ylabel有多个子图时需要统一标题格式。还有一个小细节输出图的时候不要截屏再贴到Word里不然清晰度很差。我建议用MATLAB的 exportgraphics 或者 saveas把图存成 300dpi 的 PNG 再插入文档字迹和线条都清楚很多。4.4 几个值得做的进阶验证如果你想在报告中体现一些超出课件的思考我强烈建议试一下以下两个验证实验代码量很小但分析空间很大。第一个是验证“初值敏感性”。取 r3.9从 x00.300 和 x00.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_{n1} x_{n1} - x^* f(x_n) - f(x^) \approx f(x^) (x_n - x^) f(x^) e_n)所以每迭代一次扰动大约被乘以 (f(x^*))。这个值绝对值小于1扰动会逐次缩小迭代收敛大于1扰动被不断放大即使初始值离不动点非常近也会被推离。我在报告中建议不只是写公式还要画一个小图用 (f(x)x^2c) 在不动点处的切线斜率来说明对比 c-0.5 和 c2 两种情况的斜率大小直观得多。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, ,, colornavy, markersize0.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.^2c 支持良好但 older 版本中 fzero 求解行为可能不同如果报错改用 fsolve 或自己写二分法也能达到同样目的。Python 端最容易遇到的问题是中文字体显示成方框。在绘图前加上plt.rcParams[font.sans-serif] [SimHei] plt.rcParams[axes.unicode_minus] False否则标题 “分岔图” 会变成一堆乱码。这个不算什么高深技术但每年都有人在这里卡半天。7. 从做实验到真正理解迭代做完这套实验我最大的感触就是数学里很多概念光靠黑板推导很难在脑子里“活”起来但当你看到那条 Logistic 映射分岔图一点点从一根线长成两支、四支、八支最后变成一片混沌区域的时候你才真正理解为什么物理学家费根鲍姆发现那个常数时会那么兴奋——原来这背后有一种跨越具体系统的普适规律而你的电脑屏幕上就能复现它。如果你只是照着参考答案把图跑出来交了作业说实话你亏了。我建议你在交完报告之后自己再花一个小时把参数 r 慢慢从2.5推到4.0一次只改一点点盯着那个序列值的变化那种从有序走向混沌的渐变感比任何文字描述都有说服力。这也是这门实验藏在作业背后的真正目的。