企业官网建设流程全解析

简介：这份代码包面向从事非线性动力学、混沌系统与数据驱动建模的研究生和科研人员，用于复现慢时间尺度动态发现流程，并配合SINDy方法完成玩具模型及天体、Logistic、混沌等示例的数值实验。压缩包共10个文件，以7个MATLAB脚本（.m）为核心，辅以2个.mat观测数据和1份README说明，整体约24.17MB，脚本覆盖数值积分生成数据、连续时间SINDy模型拟合、粗粒度演化离散映射发现与慢速预测等环节。计算基于公开SINDy架构，快速周期部分可结合滑动窗口DMD方法处理，便于读者对照论文复现数据与图表，理解多尺度分解、模型发现和预测验证的完整链路。已有381人学习，适合具备MATLAB基础、希望将稀疏回归与动力系统发现落到代码实践中的读者。

1. 从一段抖动的轨迹里抠出慢方程：Slow-Discovery 与它的 logistic 混沌测试床

手头有一组多尺度仿真轨迹：快变量在毫秒级来回抖动，慢变量在秒级缓慢漂移。直接拿全部状态做回归，拟合出来的方程几乎全是噪声项；如果只盯着慢变量看，又不知道它到底受哪几项支配。Slow-Discovery 这份 MATLAB 代码处理的就是这一步：用 SINDy（稀疏识别非线性动力学）从数据里把慢时间尺度上的控制方程抠出来，再配一段 logistic 混沌映射代码当测试床，验证整条管线在强非线性、对初值极度敏感的数据上还能不能收敛。它适合两类人：手上有轨迹、想反推低维慢模型的研究型开发者；以及想拿 SINDy 当入门脚手架、先把流程跑通再换自己数据的人。纯 MATLAB 实现，不依赖深度学习框架，这份代码是系统开源的结构，求解器、示例和数据接口分开放，改起来不用动底层。

2. 项目定位与技术选型：Slow-SINDy 为什么用稀疏回归而不是神经网络

2.1 多尺度系统里，慢变量其实是被"平均"出来的

一个典型的多尺度系统长这样：快变量 $y$ 以很高的频率振荡，慢变量 $x$ 在长时间上缓慢演化，两者耦合在一起。真正决定系统长期行为的是 $x$ 的方程，但 $x$ 的方程里往往含着一堆快变量的平均效应，写不出解析形式。传统做法是解析平均或者中心流形约化，前提是你得知道快方程的解析结构；工程上更常见的情况是快方程根本写不出来，只有仿真数据。

Slow-Discovery 走的是另一条路：不碰解析推导，直接把轨迹喂给稀疏回归，让数据告诉你慢方程里有哪几项。它成立的前提有两个，而且都挺硬：时间尺度分离要足够大（常见经验是快慢频率比在 10 以上，差得越多越好），采样率要能覆盖快变量，否则混叠之后快变量的高频信息会折叠进低频，慢方程怎么拟合都是错的。这两条不满足，后面所有调参都是玄学。

代码里对这一点的处理是分层的：数据接口层负责重采样和切片，中间层负责构造候选函数库，最上层才是稀疏求解。这样你换一个尺度分离比不同的系统时，只需要动前面两层。

2.2 SINDy 的内核只有一行线性代数

SINDy 的整个思想可以压成一个式子：$\dot{X} = \Theta(X)\Xi$。$\Theta(X)$ 是候选函数库，每一列是一个候选项，比如常数、$x_1$、$x_1^2$、$x_1x_2$、$\sin(x_1)$；$\Xi$ 是待求的系数矩阵，每一列对应一个状态量的方程。求解就是带稀疏约束的最小二乘，理想目标是让 $\Xi$ 里大部分元素为零，剩下的非零项就是"真实"的动力学项。

选择稀疏回归而不是神经网络，理由很实际。第一，系数可以直接读成方程，你拿到的是 $\dot{x} = -0.98x + 1.02x^2y$ 这种东西，而不是一堆权重；第二，小样本友好，几百个采样点就能出结果，神经网络在这个数据量上基本学不出物理；第三，可验证，你可以把学到的系数和已知的真值逐项对比，误差是确定的数字。代价是候选库必须包含真实项，真实项如果不在库里（比如某个特殊函数），SINDy 永远找不出来——这是它最大的局限，改代码之前必须先想清楚。

实际求解很少直接上 L1 正则，因为那要调惩罚系数，而且系数会有压缩偏差。工程上更常见的是 STLSQ（序列阈值最小二乘）：先做一次普通最小二乘，把绝对值小于阈值的系数硬置零，再对剩下的项重做最小二乘，反复迭代到系数不再变化。这份代码用的就是 STLSQ，好处是阈值有明确的物理量纲含义，比 L1 的 $\lambda$ 好调。

2.3 目录结构与 MATLAB 环境

按这类项目的常见组织方式，目录大致是下面这个形状（具体文件名以你下载到的包为准，结构上不会差太多）：

环境上，MATLAB R2019b 之后的版本都能跑，核心求解只用mldivide和pinv，不需要额外的工具箱。只有在你开启 Savitzky-Golay 平滑求导时会用到sgolay，这个函数在信号处理工具箱里，没有的话代码里一般也留了基于polyfit的滑动窗口实现。跑之前先确认工作目录切到了项目根目录，否则addpath找不到src，报的错会是"未定义函数"，新手很容易在这里卡住。

3. 核心链路拆解：一条轨迹从求导到稀疏系数的四步

3.1 入口脚本与配置项

入口脚本不干活，只负责把参数摆清楚、把四步串起来。这种写法看着啰嗦，但改参数的时候不用在几个文件里来回翻，血泪经验是：凡是把配置散落在各处的项目，三个月后自己都看不懂当时设了什么。

%% main.m —— 整条发现链路的入口 clc; clear; close all; addpath(genpath('src')); % 把求解器目录加进搜索路径 cfg = struct(); cfg.dt = 0.01; % 采样间隔，必须和轨迹文件里的真实步长一致 cfg.polyOrder = 3; % 候选库多项式最高阶，一般 2~3 足够 cfg.threshold = 0.05; % STLSQ 稀疏阈值，相对归一化后的系数 cfg.maxIter = 20; % 阈值迭代次数上限 cfg.smoothWin = 11; % Savitzky-Golay 窗口点数，奇数 cfg.polyFit = 3; % 平滑用的多项式阶数 X = load_trajectory('data/multiscale_demo.mat'); % N x d，每列一个状态量 dX = derivative(X, cfg.dt, cfg.smoothWin, cfg.polyFit); % 平滑微分 [Theta, names, scale] = build_library(X, cfg.polyOrder); % 候选库 + 列名 + 归一化尺度 Xi = stlsq(Theta, dX, cfg.threshold, cfg.maxIter); % 稀疏系数 Xi = Xi ./ scale; % 反归一化回原始量纲 print_equations(Xi, names); % 打印成可读方程

cfg.dt是最容易设错的参数：很多人从仿真里导出数据时顺手写了 1，结果学出来的系数整体差了一个数量级，还以为是算法有问题。cfg.threshold是归一化之后的阈值，所以它和状态量的物理单位无关，换了系统也不用重调，这一点比 L1 正则友好得多。

3.2 求导：整条链路最容易翻车的一步

SINDy 拟合的是 $\dot{X}$，但采样数据里只有 $X$。用一阶差分求导会把噪声放大约 $1/\Delta t$ 倍，噪声稍大一点，候选库拟合出来的就全是噪声项。常见做法是先用 Savitzky-Golay 滤波做局部多项式拟合，再对拟合多项式求解析导数，这样平滑和求导一步完成。

function dX = derivative(X, dt, win, order) % 输入：X 为 N x d 轨迹，dt 采样间隔，win 窗口点数（奇数），order 局部多项式阶数 % 输出：dX 为同尺寸的导数估计 if mod(win, 2) == 0 win = win + 1; % 窗口必须是奇数，否则无法居中 end [~, g] = sgolay(order, win); % g 的每一列对应 0~order 阶导数的卷积核 dX = conv2(X, g(:, 2)', 'same') / dt; % 第 2 列是一阶导核，再除以 dt end

窗口win和阶数order是一对矛盾：窗口越大平滑越狠，但会把慢变量的真实曲率一起抹掉；阶数越高保真度越好，但对噪声越敏感。经验值是order = 3、win取轨迹点数的百分之一到五十分之一，同时保证窗口内状态变化不超过整个轨迹幅度的十分之一。另一个坑在边界：conv2的'same'模式在两端是补零的，首尾各win/2个点的导数不可信，做验证的时候要把这段切掉，不然误差统计会被边界拖垮。

3.3 候选函数库：列归一化比阶数更重要

候选库决定了 SINDy 能搜到哪一类方程，构造逻辑就是枚举多项式项、可选地加三角项和常数项。但真正影响结果的是列归一化：如果某个状态量的量级是 $10^{-3}$，它的平方项就是 $10^{-6}$，在最小二乘里几乎不起作用，阈值一卡就全被置零了，看起来像是"系统里没有二次项"。

function [Theta, names, scale] = build_library(X, order) % 枚举到 order 阶的多项式项，外加常数项，返回归一化后的库和每列的尺度 [N, d] = size(X); cols = {}; names = {}; for total = 0:order % 按总阶数枚举，避免组合爆炸 idx = nchoosek_terms(d, total); % 返回所有和为 total 的指数组合 for k = 1:size(idx, 1) e = idx(k, :); cols{end+1} = prod(X .^ e, 2); % 逐列相乘得到单项式 names{end+1} = make_name(e); % 生成 x1^2*x2 这种可读名字 end end Theta = [ones(N, 1), cell2mat(cols)]; % 第 1 列固定为常数项 scale = max(abs(Theta), [], 1); % 每列的最大绝对值作为尺度 scale(scale == 0) = 1; % 防止常数列或全零列触发除零 Theta = Theta ./ scale; % 归一化，让各列在同一量级上竞争 end

常数项提前算进scale里是有意为之，它保证了阈值对所有列一视同仁。另一个细节是nchoosek_terms按总阶数枚举而不是按嵌套循环，状态量维度一高，嵌套循环的写法会直接写死，而按总阶数枚举只要改一个参数。归一化尺度必须原样返回给主流程，最后反乘回系数，这一步忘了，打印出来的方程系数就是错的。

3.4 STLSQ：阈值迭代与系数回收

求解器本身不长，逻辑就是"最小二乘 → 卡阈值 → 再最小二乘"。关键在于每一列独立迭代，因为不同状态量的方程稀疏程度不一样，共用一套掩码会把某一列的有效项误删。

function Xi = stlsq(Theta, dX, lambda, maxIter) % Theta: N x p 候选库；dX: N x d 导数；lambda: 稀疏阈值 Xi = pinv(Theta) * dX; % 初始解：所有项都保留的最小二乘 for iter = 1:maxIter mask = abs(Xi) < lambda; % 小于阈值的项标记为待删除 Xi(mask) = 0; % 硬阈值，直接置零 changed = false; for j = 1:size(Xi, 2) % 逐列回收，各列独立 keep = ~mask(:, j); if any(keep) sol = Theta(:, keep) \ dX(:, j); % 只对保留下来的项重解 if any(abs(sol - Xi(keep, j)) > 1e-10) changed = true; end Xi(keep, j) = sol; end end if ~changed break; % 系数稳定，提前退出 end end end

pinv而不是\做初始解，是因为候选库里常有高度相关的列（$x$ 和 $x^2$ 在窄区间上就近似线性相关），\会给出警告甚至病态解，pinv的伪逆至少能给出最小范数解，方便后续阈值筛选。lambda的调法是从大到小扫：先给一个大值让所有系数归零，再逐步减小，记录每档的非零项个数，画出"非零项数—残差"曲线，拐点就是合适的阈值。一次性猜一个阈值然后反复重跑，是效率最低的做法。

4. logistic 混沌映射测试床：用已知答案验证整条管线

4.1 生成序列与参数区间

logistic 映射 $x_{k+1} = r x_k (1 - x_k)$ 是验证回归管线的好材料：形式极简、真值已知、参数一调行为就从收敛切到混沌，非常适合做对照实验。

function x = logistic_series(r, x0, N, burn) % r: 增长参数；x0: 初值（建议落在 0.1~0.9）；N: 需要的点数；burn: 丢弃的暂态步数 x = zeros(N + burn, 1); x(1) = x0; for k = 2:(N + burn) x(k) = r * x(k-1) * (1 - x(k-1)); % 二次映射，注意迭代顺序 end x = x(burn + 1:end); % 丢掉暂态，只留稳态段 end

参数区间值得记一下：$r < 3$ 时序列收敛到不动点；$3 < r < 1+\sqrt{6} \approx 3.449$ 出现二周期；$r$ 继续增大进入倍周期分岔；到 $r \approx 3.57$ 之后大范围混沌。中间还夹着周期窗口，比如 $r \approx 3.83$ 附近是明显的三周期。burn这一步不能省，从初值到吸引子的暂态段不属于动力系统的稳态行为，混进去之后回归出来的系数会带系统偏差。

4.2 离散映射回归：绕开求导这一步

logistic 是离散映射，不是连续系统，所以不该套用"求导 + 连续 SINDy"的流程，正确做法是把 $x_{k+1}$ 当作 $x_k$ 的函数直接回归：$x_{k+1} = \Theta(x_k)\Xi$。这样完全避开了数值求导，噪声不被放大。

r = 3.9; x0 = 0.31; x = logistic_series(r, x0, 2000, 500); % 构造回归对：用前一时刻预测后一时刻 Xnow = x(1:end-1); Xnext = x(2:end); % 只保留常数项、一次项、二次项，因为 logistic 正好是二次的 Theta = [ones(numel(Xnow),1), Xnow, Xnow.^2]; Xi = stlsq(Theta, Xnext, 1e-3, 20); r_hat = -Xi(3); % x^2 项系数是 -r c_hat = Xi(1); % 常数项应为 0 fprintf('估计 r = %.6f，真值 r = %.6f，常数项 = %.2e\n', r_hat, r, c_hat);

这里有一个很爽的自检点：一次项系数和二次项系数应该严格互为相反数，都等于 $\pm r$，常数项应该被稀疏化掉。跑出来的结果如果一次项系数是 3.9、二次项系数是 -3.9、常数项小于 $10^{-12}$，说明求导之外的整条链路（库构造、归一化、阈值迭代、系数回收）都是对的。反过来，如果常数项跑出一个非零值，八成是burn没做够，或者数据里混进了别的段。

4.3 噪声、采样与初值敏感性对照

真正拉开差距的是加噪和降采样这两组实验。给序列加不同水平的高斯噪声，看估计出来的 $r$ 偏多少：

趋势很清楚：噪声在 $10^{-2}$ 量级以下，二次结构还稳得住；到 $5\times10^{-2}$ 时，$r$ 的估计已经明显偏低，而且开始冒出多余项——这是稀疏回归的典型行为，噪声被候选库里的高阶项"接住"，变成了假动力学项。这时候该做的是加大平滑（连续系统）或者增加采样点做多次平均（离散映射），而不是一味调小阈值，调小只会让假项更多。

降采样是另一个隐蔽的坑。把序列每隔 5 个点取一个，虽然仍然是 $x_{k+1}$ 对 $x_k$ 的回归，但实际回归关系变成了 $x_{k+5} = f^{(5)}(x_k)$，那是五次复合映射，用二次库去拟合当然对不上。表现为估计的 $r$ 完全离谱、残差下不去，而人会误以为是代码 bug。处理办法很简单：回归对之间的步长必须和真实采样间隔一致，降采样之后要么重建库、要么把复合映射的真值算出来做对照。

初值敏感性这条要单独说。logistic 在混沌区对初值是指数敏感的，两条相差 $10^{-8}$ 的初值，几十步之后就完全分开。所以验证 SINDy 结果时，绝对不能用"长时积分后轨迹是否重合"作为唯一指标，否则连真值方程都通不过。可行的指标只有两个：短期（几步到十几步）预测误差，以及系数本身的稳定性——把轨迹切成若干段分别回归，看同一项的系数在各段之间抖不抖。系数稳、短期误差小，就是可信的结果。

5. 避坑与排查：方程学不出来时的 5 个具体原因

5.1 拿到结果先看残差，别看系数

很多人跑完第一件事是盯着系数矩阵找规律，但系数是否可信，得先看拟合残差 $| \dot{X} - \Theta(X)\Xi |$ 的分布。残差如果集中在轨迹的某一段（比如快变量幅值最大的地方），说明这一段动力学没被候选库覆盖，或者求导在那里失真；残差如果整体均匀但数值偏大，通常是库的阶数不够或阈值过大。我一般会先画残差随时间的曲线，再去读系数，顺序反了容易把假项当发现。

5.2 系数全是零

现象：跑完stlsq，Xi整个矩阵全零，或者只剩常数项。原因通常是阈值设得过大，或者归一化没生效——某一列的尺度是 $10^{-6}$，归一化漏了，它在最小二乘里的贡献本来就接近零，阈值一卡必然被删。解决：先把lambda调到 $10^{-6}$ 量级确认能出非零解，再从大到小扫一遍；同时打印scale，检查有没有列的最大值远小于其他列。如果某列全零，那说明这一列在数据上恒为零，候选库里本来就不该存在。

5.3 系数看着对、数值差一个数量级

现象：方程的项选对了，但系数整体偏大或偏小十倍、百倍。原因基本是两个：一是cfg.dt和真实采样间隔不一致，求导时除错了；二是归一化之后忘了反乘scale。前者在从仿真软件导出数据时极其常见，导出脚本里写dt = 1，实际步长是 0.01，系数就整体差 100 倍。后者是纯粹的程序疏漏，检查一下主流程里有没有Xi = Xi ./ scale这一行。

5.4 一步预测很准，长时积分直接发散

现象：用学到的方程做单步预测误差极小，但积分几十步之后轨迹跑飞。原因多半不在代码，而在系统本身是混沌的——初值的微小误差被指数放大，任何模型都撑不住长时间积分。解决思路是换评价体系：改成多起点、短窗口的预测误差统计，或者直接比较系数与真值的相对误差。如果你的目标系统不是混沌的，那才需要怀疑模型漏项，这时候回到残差分析，看漏掉的是哪一段动力学。

5.5 高阶项系数巨大且不稳定

现象：候选库里本来只想要二次项，结果三次、四次项的系数跑出几百上千，重跑一次符号还变。原因几乎都是共线性：数据在相空间里只覆盖了很窄的一块，$x$ 和 $x^2$ 在这段区间上高度线性相关，最小二乘解不唯一，系数就会乱跳。解决：扩大激励范围或者换一段数据，让状态量覆盖更宽的区域；再不行就在库里去掉明显冗余的项，或者先对候选库做一次 QR 分解剔掉相关列。这一点上，数据质量比算法调参重要得多。

6. 二次开发：把 Slow-Discovery 换到自己的数据上

这份代码的真正价值不在 demo 跑出来的那条方程，而在它的三层结构：数据接口、候选库、求解器。换系统时，改动点基本集中在前两层。

第一个替换点是数据接口。全流程约定X是一个 $N \times d$ 的矩阵，每列一个状态量，dt是标量采样间隔。你要做的就是把load_trajectory换掉，保证输出满足这个约定，时间列如果是非均匀的，先重采样到均匀网格。这一步不动求解器，改完直接用demo_logistic的回归测试验证接口是否正确。

第二个替换点是候选库。build_library里的多项式枚举很容易扩展，加三角项、加绝对值项、加某个已知的物理项都可以，只要保证返回的names和列一一对应、scale同步更新。扩展的时候注意别把库撑得太大：$d=5$、阶数 3 的时候列数已经上百，再往上加组合项，样本量不够会直接过拟合。

第三个替换点是求解器。stlsq的接口是(Theta, dX, lambda, maxIter)，想换成 L1 正则、ADMM 或者加权最小二乘，只要保持这个签名不变，main.m一行都不用动。

改完怎么验证？我的习惯是先用已知系数的合成方程做一遍回归测试：自己造一个 $\dot{x} = -x + 0.5x^2$ 的系统，用ode45生成轨迹，加白噪声，再跑整条管线，看能不能把 $-1.0$ 和 $0.5$ 原样找回来。这一步过了，再上真实数据。真实数据永远有你想不到的问题，拿它当第一轮测试，出了问题根本分不清是代码错还是数据脏。从那以后我每次动求解器或者候选库，都强制走一遍这个合成测试，宁可多花十分钟，也不想再经历一次把假项写成结论、两周后回头推翻的翻车。希望这个习惯也能帮到你。

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