MATLAB SINDy稀疏辨识：慢时间尺度与logistic混沌代码解析

简介：面向混沌与复杂动力系统数据驱动建模的学习者与科研人员，这套 MATLAB 代码围绕 SINDy 架构实现慢时间尺度动态发现，可用于复现 Bramburger、Dylewsky 与 Kutz 在 Physical Review E（2020）中的数据和图形，并借助滑动窗口 DMD 提取快速周期成分。内容覆盖玩具模型数值积分、连续时间 SINDy 模型拟合、粗粒化演化的离散时间映射发现及慢速预测流程，也提供 Jupiter、Saturn、Logistic 与 Chaos 等示例场景，便于理解从高维数据中分离快慢尺度、构建可解释模型的完整思路。代码将快速周期识别与慢尺度建模衔接起来，形成从数据生成、稀疏回归到预测验证的闭环。资源包共 10 个文件，含 7 个 .m 脚本、2 个 .mat 数据文件和 1 个 .md 说明文档，压缩包约 24.17MB，结构清晰，适合按 README 指引逐步运行与对照实验。已有 380 人学习下载，可作为非线性动力学、模型发现与论文复现训练的实践材料。

1. 慢时间尺度上的稀疏辨识：这份 logistic 混沌 MATLAB 代码包在做什么

包名挂着 logistic 混沌 MATLAB 代码，解压后翻两页就会发现真正的主角是 SINDy：不从方程出发，把轨迹数据丢进候选函数库，用稀疏回归把右端项一项一项挑出来，挑不中的系数压成零。Slow-Discovery 这一版把重心压在慢时间尺度上，先滤掉快振荡，再在慢流形上做辨识。

代码量不大，几个 .m 文件加两组示例数据，依赖基本只有 MATLAB 自带函数和信号处理工具箱。适合手上有仿真轨迹、想让模型自己说话的人；如果你的数据只有几百个点还全是噪声，前面两章的导数估计会先给你上一课。

2. 结构与调用链：五个函数撑起一条从轨迹到方程的流水线

这套代码的价值不在示例本身，而在它把「辨识」拆成了边界清楚的几段。你换系统、换数据，也只需要替换其中一两段，剩下的骨架不动。先看目录，再看一条完整调用是怎么走完的。

2.1 目录布局：lib 里每个文件负责一段

这类示例包的目录通常很扁平，我拆下来的版本大致长这样。不同提交之间文件名会改，但分层逻辑没动过：数据一份、算法一组、示例脚本若干，画图单独拎出来。

Slow-Discovery/ ├── startup.m % 把 lib/ 和 data/ 加进 path，跑示例前先执行一次 ├── data/ │ ├── logistic_r4.mat % r=4 的混沌序列，5000 步 │ └── vdp_mu50.mat % 刚性快慢系统积分结果，等距采样 ├── lib/ │ ├── build_library.m % 拼候选函数库 Θ(X) │ ├── exponent_matrix.m % 多指标枚举，给 build_library 打底 │ ├── derivative_sg.m % Savitzky-Golay 平滑求导 │ ├── stlsq.m % 序贯阈值最小二乘 │ └── rollout_error.m % 用辨识出的模型积分，和真值比对 ├── examples/ │ ├── ex1_logistic_map.m % 离散映射，验证链路对不对 │ ├── ex2_slow_fast_vdp.m % 快慢分离，慢流形辨识 │ └── ex3_noisy_observer.m % 加噪声版本 └── plots/ └── plot_coefficients.m % 系数条形图，一眼看出哪些项活下来了

几个文件的分工值得单独说清楚，因为它们对应辨识流程里最容易出问题的几个环节：

examples这一层是可以整体替换的，它只负责组织数据、调用 lib、打印结果。真正的算法都在 lib 里，改 lib 会影响全部示例，改 examples 不会互相牵连。这种切法在学术代码里算克制了，很多同规模的包把所有东西塞进一个脚本，改一个参数要滚三百行。

2.2 主链路：五步走完一次辨识

不管换成什么系统，这条链路是不变的。SINDy 的核心假设就一句话：真实的右端项是候选函数库里少数几项的线性组合，所以回归完再阈值化，剩下的就是方程。

先用离散的 logistic 映射跑一遍，因为它有解析解可以对账。映射x_{k+1} = r·x_k·(1-x_k)没有导数，用一步增量代替，库只需要[1, x, x²]三列：

% examples/ex1_logistic_map.m rng(0); % 混沌系统对初值极敏感，随机数不固定就没法比对结果 r = 4.0; % 落在混沌区，轨迹能遍历 (0,1) 整段，回归的条件数才好 N = 5000; x = zeros(N,1); x(1) = 0.2; for k = 1:N-1 x(k+1) = r * x(k) * (1 - x(k)); end X = x(1:end-1); % 当前状态 dXdt = x(2:end) - x(1:end-1); % 映射没有导数，用一步增量顶上 [Theta, names] = build_library(X, 2, false, true); % 三列：1, x, x^2 [Xi, hist] = stlsq(Theta, dXdt, 1e-6, 10); for k = 1:numel(names) fprintf('%-8s %+9.5f\n', names{k}, Xi(k)); end fprintf('迭代过程中的非零项数：%s\n', mat2str(hist'));

跑完屏幕上应该是1 → 0.00000、x1 → +3.00000、x1^2 → -4.00000。这个结果不是巧合：把x_{k+1} - x_k展开就是(r-1)x - r x²，代进r = 4正好是3x - 4x²。系数能精确对上，说明库、导数、回归三段都没写错。

拿到这个对照之后，再去看ex2才知道偏差是数据本身的问题还是代码的问题。这也是我建议的跑法：先用有解析解的算例把链路验一遍，别一上来就拿真实数据调参。

2.3 跑起来：三条命令和两个检查点

在 MATLAB 命令窗口里切到包目录，先跑startup把路径加进去，再直接调示例脚本名。脚本是函数化的，不会在 workspace 里留一堆同名变量。

cd Slow-Discovery startup % 把 lib/ 和 data/ 加进搜索路径 ex1_logistic_map % 第一个算例：离散映射，用来对账 ex2_slow_fast_vdp % 第二个算例：刚性快慢系统，慢流形辨识

第一个检查点是系数表。ex1里的三个系数如果差得离谱，别再往下跑，先回头查导数估计和库的列顺序。第二个检查点是hist里的项数序列。正常情况它应该在前几轮快速下降然后停住，比如[20 9 4 4]。如果一直不变，说明阈值根本没起作用；如果一路掉到 1，说明阈值太狠，把真项也杀了。

提示：sgolayfilt和sgolay属于信号处理工具箱。只有基础 MATLAB 的话，把derivative_sg.m换成中心差分也能跑，只是噪声放大得厉害，ex3那个带噪算例会直接失效。

3. 核心模块拆解：库、导数、稀疏回归三段怎么写

SINDy 看着公式简单，落地时三个地方最容易翻车：库怎么枚举、导数怎么估、阈值怎么给。这三段决定了你最后拿到的是方程还是一堆数值噪声。

3.1 候选函数库：多指标枚举与列名回溯

库的本质是一个矩阵，每一列是一个候选函数在采样点上的取值。1 个变量 3 阶就是 4 列，3 个变量 3 阶就变成 20 列，10 个变量 3 阶直接跳到 286 列。枚举写错一位，后面全废。

function E = exponent_matrix(n, d) % 枚举 n 个变量、总次数不超过 d 的所有指数向量，一行一个 % 返回 E 的第 k 行是第 k 个单项式各变量的幂次 E = zeros(1, n); % 常数项，全零 for k = 1:d prev = E(sum(E,2) == k-1, :); % 上一轮里总次数恰为 k-1 的那些 add = []; for j = 1:n tmp = prev; tmp(:, j) = tmp(:, j) + 1; % 第 j 个变量次数加一 add = [add; tmp]; end E = unique([E; add], 'rows'); % 去重，避免 x1*x2 和 x2*x1 算两次 end end

build_library.m拿这个指数矩阵逐行展开成列：外层遍历指数向量，内层遍历变量，把X(:,j).^E(k,j)累乘起来。列名也顺手拼出来，后面打表、画条形图都靠它对齐，否则你看到Xi(7) = -0.33根本不知道是哪一项。

function [Theta, names] = build_library(X, polyOrder, useTrig, useConst) [N, n] = size(X); E = exponent_matrix(n, polyOrder); cols = {}; names = {}; if useConst cols{end+1} = ones(N,1); names{end+1} = '1'; end for k = 1:size(E,1) if all(E(k,:) == 0), continue; end % 常数项已经在上面加过 col = ones(N,1); for j = 1:n col = col .* X(:,j).^E(k,j); % 逐变量次幂累乘 end cols{end+1} = col; names{end+1} = mat2str(E(k,:)); end if useTrig for j = 1:n cols{end+1} = sin(X(:,j)); names{end+1} = sprintf('sin(x%d)', j); cols{end+1} = cos(X(:,j)); names{end+1} = sprintf('cos(x%d)', j); end end Theta = [cols{:}]; end

polyOrder是最需要克制的参数。阶数越高，能表达的非线性越丰富，但列数按组合数涨，而且高阶列的量级往往比常数项大好几个数量级，不归一化的话最小二乘会把权重全压到高阶列上。我一般从 2 阶开始试，能复现再往上加。useTrig只在数据里有明显周期成分时才打开，混沌系统和刚性振子系统基本用不上，加了反而给回归添乱。

3.2 导数估计：Savitzky-Golay 比直接差分强在哪

SINDy 最脆弱的一环就是导数。状态是测出来的，导数是算出来的，一步差分等于把噪声放大1/Δt倍。ex3那个算例就是专门演示这件事的：同一份数据，直接差分求出来的系数能差出两倍。

function dX = derivative_sg(X, dt, win, polyOrder) % 在滑动窗口内做局部多项式拟合，同时拿到平滑值和一阶导 % win 必须为奇数，polyOrder 要小于 win if mod(win,2) == 0 win = win + 1; end dX = zeros(size(X)); [~, g] = sgolay(polyOrder, win); % g 的列对应 0 阶、1 阶、2 阶导的卷积核 for j = 1:size(X,2) dX(:,j) = conv(X(:,j), -g(:,2), 'same') / dt; end edge = (win-1)/2; % 两端卷积边界不可靠，裁掉 dX(1:edge, :) = []; end

窗口长度win的选法是最玄学的一步。窗口太窄，噪声没压住，导数抖得跟毛刺一样，回归会硬生生多拟合出几项；窗口太宽，真实的快速变化被抹平，辨识出来的方程对不上。判断标准是：窗口要能装下一个完整的高频周期，但明显短于你关心的那个时间尺度。快慢系统里这个约束很紧，因为两个尺度差了几十倍，经常得先降采样再求导。-g(:,2)前面那个负号是符号约定，第一次用最好拿x = sin(t)试一下，导数核输出对了再往真数据上套。

3.3 稀疏回归：列归一化加硬阈值的两处关键处理

标准 SINDy 用的是序贯阈值最小二乘，思路很朴素：先普通最小二乘求一版系数，把绝对值小于阈值的项清零，再对剩下的项重做一次最小二乘，反复几轮直到项数不再变。

function [Xi, hist] = stlsq(Theta, dXdt, lambda, maxIter) % Theta : N×p 候选函数库；dXdt : N×n 导数 % lambda: 硬阈值，归一化之后无量纲，一般 1e-3 ~ 1e-1 [N, p] = size(Theta); scale = vecnorm(Theta, 2, 1); % 每列的 2-范数 scale(scale < eps) = 1; Tn = Theta ./ scale; % 关键一步：列归一化 Xi = Tn \ dXdt; % 首轮退化为普通最小二乘 hist = nan(maxIter,1); for it = 1:maxIter Xi(abs(Xi) < lambda) = 0; % 阈值化 for j = 1:size(dXdt,2) active = Xi(:,j) ~= 0; if any(active) Xi(active,j) = Tn(:,active) \ dXdt(:,j); % 只对活着的项重解 end end hist(it) = nnz(Xi); if it > 1 && hist(it) == hist(it-1), break; end % 结构不再变就停 end Xi = Xi ./ scale'; % 系数换回原量纲 hist = hist(1:it); end

列归一化这一步容易被跳过，但跳过之后lambda就失去意义了。x1的取值在 1 附近，x1^3能到几十，两列的范数差一个量级，同一个阈值对前者是「很小」，对后者可能是「很大」。归一化之后所有列都在同一尺度上，lambda才变成一个可以跨系统迁移的参数——我一般直接设在1e-2左右起步，再上下扫一个数量级。

hist这个输出比系数本身更有诊断价值。正常收敛是单调下降到某个数就停住；来回跳说明阈值卡在两项系数的量级之间，得微调；一直不降说明lambda太小或者导数噪声太大，回归在硬拟合噪声。

注意：当两列近似共线时（比如数据只覆盖很窄的一段状态空间），Tn(:,active) \ dXdt(:,j)会给出很大的系数，不同初值跑出来的方程完全不一样。这种情况不是调参能救的，得多采几段不同初值的数据拼起来。

4. 慢时间尺度发现：从快振荡里把慢流形剥出来

这一章是这个包真正的重点。前面那套流水线直接套在快慢系统上也能跑，跑完误差还挺小，但那个模型没有任何意义。慢时间尺度发现要解决的是：怎么让回归只看见慢的那一层。

4.1 为什么直接拟合会学歪

以经典的刚性快慢振子为例，μ取 50：

ẋ₁ = x₂ ẋ₂ = μ(1 - x₁²)x₂ - x₁

x₁是慢变量，振幅在 2 附近慢慢变；x₂大部分时间贴在慢流形x₂ ≈ x₁³/3 - x₁上，只在转折点附近瞬时跳变。你把原始高频采样的(x₁, x₂)直接丢进二阶库，回归确实能拟合出东西，但学到的是一条把跳变也硬描下来的曲线，项数多、系数大，往外推几步就崩。因为真实的慢动态只有ẋ₁ = x₁³/3 - x₁这一条，剩下的全是快层瞬态。

辨识的时间尺度必须和数据的时间尺度对齐。这一点是所有多尺度数据驱动方法的共同前提，SINDy 尤其敏感，因为它连项数都要求稀疏。

4.2 低通滤波加慢网格重采样

处理办法分两步。先用 Savitzky-Golay 平滑把快跳变抹掉，窗长要远大于快层的特征时间；再在慢网格上重采样，让相邻样本之间的变化体现的是慢流形上的演化。

% examples/ex2_slow_fast_vdp.m mu = 50; % 越大快慢分离越明显 f = @(t,x) [x(2); mu*(1 - x(1)^2)*x(2) - x(1)]; tspan = linspace(0, 2000, 400001); % 强制等距，不让求解器自己挑点 opts = odeset('RelTol',1e-9,'AbsTol',1e-11); [t, X] = ode15s(f, tspan, [2; 0], opts); % 刚性系统用 ode15s，ode45 会跑到天黑 dt = t(2) - t(1); xs = sgolayfilt(X(:,1), 2, 2001); % 窗长远大于快跳变时间 ts = t(1):200*dt:t(end); % 慢网格，一步跨过一个慢相位 Xs = interp1(t, xs, ts)'; dXs = derivative_sg(Xs, ts(2)-ts(1), 21, 3); Theta = build_library(Xs, 3, false, true); % 慢流形上三次项刚好够 Xi = stlsq(Theta, dXs, 5e-2, 10);

tspan传等距向量的写法是刻意的。ode15s自适应步长，在你完全不指定的情况下会给出疏密不均的时间点，两个相邻点可能差 1e-6 也可能差 1e-1，直接拿来差分，等于把步长抖动当成了动态。强制等距输出是最省事的做法，实在不行就interp1二次重采样。

三次库[1, x, x², x³]辨识出来应该是[0, -1, 0, 1/3]，也就是ẋ₁ = -x₁ + x₁³/3。实测会有点偏差，μ = 50时系数能对上九成左右，μ越大越准。原因很清楚：低通滤波抹掉的是瞬态，留下的慢段占时间比例是μ量级，权重上去了，跳变段的污染就相对小了。想要更高精度就得继续加大μ，或者改用分段拟合只取慢相位那一段。

4.3 评估：rollout 误差比系数误差更值得看

系数对上了当然好，但对不上时你还需要知道模型能不能用。判断标准是拿辨识出来的方程积分一遍，和真值轨迹比。rollout_error.m干的就是这件事：

% 一维慢模型：Xi 按 [1, x, x^2, x^3] 排列 rhs = @(t, x) polyval(flipud(Xi), x); [~, Xp] = ode45(rhs, ts(1:200:end), Xs(1)); err = norm(Xp - Xtrue, 'fro') / norm(Xtrue, 'fro');

用相对 Frobenius 范数是为了让不同量级的系统能横向比。我一般看三个数：训练段误差、外推段误差、以及外推能撑多少个慢周期还不发散。

系数误差和 rollout 误差是两回事，经常出现系数差 15% 但轨线几乎重合的情况，因为小项对定性行为影响很小。反过来，如果系数完美而 rollout 三步就发散，八成是库里有共线列，回归求出来的解在参数空间里离真值很远，只是在采样点上碰巧对上了。

5.

本文还有配套的精品资源，点击获取