简介：一份面向Matlab建模学习者与数学优化实践者的变分法入门教程，系统梳理泛函、泛函变分、极值条件以及Euler-Lagrange方程等核心概念，并围绕最速降线、悬链线两个经典问题说明变分法如何将现实优化需求转化为极值求解模型。资料包内含1个doc格式文档，压缩包共617KB，文件虽小但知识密度较高，兼具理论推导与场景引入，适合数学建模竞赛备赛、物理工程优化入门及相关课程自学使用。此文档已有138人学习，内容从约翰·伯努利时代的历史溯源切入，逐步过渡到古典变分法的数学表达，并提示可借助Matlab的数值计算与可视化功能模拟验证最优曲线形态。通过阅读，读者能够理解泛函极值的基本思想，掌握Euler-Lagrange方程的推导思路，并形成将能量最优、时间最优等问题转化为标准变分形式的能力。整体安排简洁清晰，可作为使用Matlab开展变分法建模实验的前置参考。

1. 变分法不是高数余料：它是 Matlab 建模里「求最优函数」的正门

做过数学建模竞赛或者物理仿真的读者，多半会遇到一类「反直觉」的问题：要求的东西不是一个数，而是一条曲线、一个函数，目标是让某个积分最小或最大。最速降线、最短光程、悬链线、机械臂的能耗最优轨迹，全都可以装进同一个数学框架——变分法。这份「Matlab 建模教程」把这些内容做成了一条能落地的主线：先搞懂变分法的核心是什么，再在 Matlab 里把符号推导、数值求解、参数调试全部跑通。对刚装好 Matlab、还在摸索建模流程的人，这是一条不依赖额外工具包的可复现路线；对已经写过不少仿真代码的人，它也能把「最优函数」这类问题从玄学变成可以调参、可以验证的工程活。

2. 从变分到建模：把欧拉-拉格朗日方程的符号推导做成可复用脚本

2.1 泛函极值：从求一个数到求一条曲线

普通优化里，目标函数f(x)的自变量是一个数，极值条件是一阶导数为零；变分法里，目标函数J[y]的自变量是一个函数y(x)，被积函数里通常同时含x、y(x)和y'(x)，写成积分形式：

J[y] = ∫ F(x, y, y') dx

建模任务往往是给出一族候选函数，找出让J取极值的那一条。变分法的核心结论是：如果J在某个y上取极值，那么这个y必须满足欧拉-拉格朗日方程：

∂F/∂y − d/dx(∂F/∂y') = 0

这个方程在建模里太常用了。数学建模竞赛里，只要题目里出现「使旅行时间最短」「使能耗最小」「使路径长度最短」这类措辞，十有八九最终会落到这个式子上。问题在于，手推欧拉-拉格朗日方程并不难，难的是 F 一旦复杂起来，求导、链式法则、整理合并很容易出错，尤其当 F 里还混着根号、分式和非线性项。所以我的固定做法是：让 Matlab 的符号计算去推，人来检查方程对不对，而不是人充当求导引擎。

2.2 欧拉-拉格朗日方程的符号推导：一段可以直接抄的 matlab 代码

Matlab 的 Symbolic Math Toolbox 提供了functionalDerivative，可以一眼看出 F 对 y 的完整变分结果。以经典的「两点之间最短路径」为例，被积函数是弧长微元：

% 变分法符号推导：一行求出欧拉-拉格朗日方程 syms x y(x) % 声明自变量 x 和函数 y(x)，注意 y 必须是函数 F = sqrt(1 + diff(y, x)^2); % 泛函的被积函数：弧长元，求最短路径 EL = functionalDerivative(F, y) == 0; % 变分核心：对函数 y 求变分 EL = simplify(EL); % 化简，把方程整理成最简二阶 ODE disp(EL) % 输出结果：2*diff(y(x), x, x) == 0

逻辑说明：syms x y(x)和常见写法syms x y的区别在最后一行输出里体现得很明显——前者把y声明为依赖x的函数，diff(y, x)才能被识别为y'，functionalDerivative才能返回真正意义上的欧拉-拉格朗日表达式。functionalDerivative(F, y)返回的不是一个数字，而是一个表达式，它代表∂F/∂y − d/dx(∂F/∂y')；把它等于 0 就是完整的 EL 方程。

参数说明：这一步要求 Symbolic Math Toolbox 可用，老版本 Matlab 偶尔提示Undefined function时，去命令行跑ver symbolic确认工具箱存在即可。算完simplify很关键，因为原始输出里通常是一堆diff(diff(y(x),x),x)嵌套项，不化简根本没法看出y''=0这样的物理意义。上面这段代码的结果会化简成2*y'' == 0，也就是「最短路径是直线」，符合直觉。

换一个被积函数，只要把F改掉就行。比如把 F 换成y * sqrt(1 + diff(y,x)^2)，求解结果会直接指向悬链线方程。整个脚本的可复用性就在这里：以后做建模，只需要维护一个F的表达式，符号推导是一次性的。

2.3 悬链线建模：把符号解落到具体参数

悬链线这个例子很适合验证符号推导的正确性：给定两个悬挂点，求一条处处受力平衡的曲线，泛函写成势能积分，F 取y*sqrt(1+y'^2)，y 方向取向上为正。求出 EL 方程后，可以手工把候选解代回去验证是不是y = a*cosh((x-b)/a)。

% 悬链线：用符号工具箱把解析形式 y=a*cosh((x-b)/a) 代回 EL 验证 syms x y(x) a b F = y * sqrt(1 + diff(y, x)^2); % 悬链线泛函 el = functionalDerivative(F, y); % 只取方程左侧，不着急等于 0 cand = a * cosh((x - b) / a); % 手工写出的候选解 el_cand = subs(el, {y, diff(y, x)}, {cand, diff(cand, x)}); simplify(el_cand) % 期望输出 0，证明候选解确实满足 EL

逻辑说明：这里的技巧是用subs把符号函数y和它的导数y'同时替换成候选解。如果只替换y不替换diff(y,x)，化简结果里有残留的y'项，看起来永远对不上。两行subs把候选解拆成两个替换目标，是这类验证的固定写法。simplify输出 0 等价于 EL 恒成立。

参数说明：悬链线的两个参数a和b有明确的物理含义——a与缆索张力、线密度之比有关，b是最低点的横坐标。实际建模时给的是两个端点坐标，比如端点分别位于(0,1)和(2,3)，要让a*cosh((x-b)/a)穿过这两个点：

syms a b real eq1 = a * cosh((0 - b) / a) == 1; eq2 = a * cosh((2 - b) / a) == 3; ab = vpasolve([eq1, eq2], [a, b], [1; 1]); % 初值猜 a=1, b=1

这里容易踩的坑是vpasolve的初值非常影响结果，换一组端点就要重新试初值。我的经验是：先观察两个端点连线的中点位置，把b初值设成中点横坐标，a初值设成一个与「下垂程度」同量级的数，大概率一次收敛。悬链线的符号验证和参数求解做完，变分法的第一印象就建立起来了：先由 EL 拿到方程，再从方程里解出参数，这是变分建模的标准节奏。

3. 用 Matlab 数值求解变分问题：以最速降线为例走通 bvp4c

3.1 什么时候必须放弃解析解转向数值

不是所有变分问题都有漂亮解析解。最速降线的 F 里带着根号和除法，EL 方程是非线性二阶 ODE，虽然能通过参数化技巧求到解析解，但那个过程的推导量已经不适合建模竞赛里现场手推。更别提 F 里一旦混入摩擦、控制输入、分段约束，解析这条路基本就走死了。

这时可以把 EL 方程当成一个边值问题来处理：二阶 ODE 加上两个端点条件，数值求解。Matlab 内置的bvp4c就是干这个的。它的内部原理是配置法加牛顿迭代，对大多数光滑边值问题表现稳定，不需要自己写打靶循环。但它的输入参数比ode45多一层，很多人第一次用会卡在「初始猜测怎么给」上。下面用一个完整例子把整个调用链走一遍。

3.2 打靶法与 bvp4c 求解边值问题的完整脚本

最速降线问题的泛函是：

J[y] = ∫ sqrt((1+y'^2) / (2g y)) dx

其中 y 向下为正，g 是重力加速度。常数2g不影响极值点，可以合并进 F 的常数因子。由 EL 方程可以化简出一个比较紧凑的二阶 ODE：

y'' = −(1 + y'^2) / (2y)

这个形式可以直接交给bvp4c。注意边界条件：起点y(0)=0会让 ODE 右侧除以 0，所以在数值上必须把起点抬到1e-3量级，否则第一轮迭代就 NaN。

% 最速降线：最小化 J = ∫ sqrt((1+y'^2)/y) dx，y 向下为正 function brachistochrone_demo() x0 = 0; x1 = 1; % 初始猜测：20 个网格点起步，跑通后再加密 solinit = bvpinit(linspace(x0, x1, 20), [0.01; 1]); opts = bvpset('RelTol', 1e-4, 'AbsTol', 1e-6, 'NMax', 5000); sol = bvp4c(@odefun, @bcfun, solinit, opts); % 绘制数值解 plot(sol.x, sol.y(1, :), 'LineWidth', 1.5) xlabel('x'); ylabel('y'); title('最速降线数值解') end function dydx = odefun(x, y) % y(1) 是高度，y(2) 是斜率 y' dydx = [y(2); -(1 + y(2)^2) / (2 * y(1))]; end function res = bcfun(ya, yb) % 边界条件：起点高度为 1e-3，终点高度为 1 res = [ya(1) - 1e-3; yb(1) - 1]; end

逻辑说明：bvp4c需要三个输入，ODE 右端函数、边界条件函数、初始猜测结构体solinit。solinit由bvpinit生成，第二个参数[0.01; 1]是给两个状态变量[y; y']各猜一个常数值。边界条件函数bcfun传入的是起点状态ya和终点状态yb，返回值是一个列向量，每个元素代表一个边界约束。这里只约束了位置，没有约束终点斜率，因此解出来的曲线真的是「自由末端斜率」的极小值曲线。

参数说明：如果终点高度或跨度变了，只需要改x1和bcfun里的终值；如果物理模型换了，改odefun里的 ODE 表达式，其他结构可以完全复用。这里已经埋了一个典型翻车点——y(0)直接取 0 会让 ODE 里出现除零，后面第 5 章会专门展开。

3.3 bvp4c 的参数怎么设：网格、容差与初值猜测

bvp4c不像ode45那样「给个初值沿着时间推」，它是在整个区间上反复迭代修正解，所以那几个控制参数对成败的影响很大。下面这张表是我常用的设置。

失败时最常见的警告是Unable to meet the tolerance without reducing the step size，翻译过来就是迭代收敛不了。对策按优先级来：先把NMax拉上去，再放宽RelTol到 1e-3，最后还是不行就回去改初始猜测。初始猜测不需要精确，但方向和量级必须对。比如最速降线的y(1)高度在 0 到 1 之间，猜测值猜成 10 就会让 ODE 里那个分母先飙到离谱的值，迭代自然散掉。

判断是否收敛，不要只看bvp4c的返回值不报错。要画图看曲线平滑度，再把网格加密一倍、两次解之间的最大差压到 1e-4 以下。这部分在第 6 章会再展开。记住一个原则：bvp4c是一个需要你理解它行为的求解器，不是喂进去就出结果的魔法盒；一旦不理解参数的含义，报错时就只能对着黑匣子发呆。

4. 变分法在建模实战中的高频形态：最优控制与轨迹规划

4.1 从欧拉-拉格朗日到哈密顿：控制输入怎么进泛函

很多建模题看起来不是「求一条曲线」，而是「求一个控制输入随时间变化的规律」。比如小车从 A 到 B，控制量是驱动力，要最小化总能耗；机械臂从初始位形运动到目标位形，控制量是关节力矩，要最小化冲击。这类问题仍然属于变分法的范畴，只是多了一层约束：状态变量 x(t) 的变化不是自由选择的，必须满足状态方程 x' = f(x, u, t)。

把约束用拉格朗日乘子 λ(t) 放进目标泛函，就得到了哈密顿函数：

H(x, u, λ, t) = L(x, u, t) + λ' * f(x, u, t)

极值条件变成三组方程：状态方程x' = ∂H/∂λ，协态方程λ' = −∂H/∂x，以及最优性条件∂H/∂u = 0。这就是庞特里亚金极小值原理的骨架。因为篇幅，这里不展开推导，只需记住：变分法的思想换了一层包装，从「求 y(x)」变成了「求 u(t)」，解法的落地路径完全不同——解析求解基本不可能，工程上几乎全是离散化加数值优化。

4.2 不想碰解析变分？用 fmincon 做离散化轨迹规划

当控制输入受到不等式约束（比如力矩上限 5 N·m）时，解析变分和 bvp4c 都很难处理。我的常用做法是：把时间轴离散成 N 段，把每个时刻的状态和输入都当作优化变量，用 fmincon 去求。这是一个标准的直接配点法，Matlab 里写起来非常直接。

% 最优控制离散化：最小化控制能量，状态方程 x' = x + u function trajectory_planning_demo() N = 40; T = 1; dt = T / N; % 设计变量 v = [x(1..N); u(1..N)]，前 N 个是各时刻状态，后 N 个是控制 v0 = [linspace(0, 1, N)'; zeros(N, 1)]; % 初值：直线轨迹 + 零控制 lb = [-inf(N, 1); -5 * ones(N, 1)]; % 控制下限：u >= -5 ub = [ inf(N, 1); 5 * ones(N, 1)]; % 控制上限：u <= 5 % 目标：控制能量的矩形法离散化 cost = @(v) dt * sum(v(N+1:end).^2); opts = optimoptions('fmincon', 'Display', 'iter', ... 'MaxIterations', 500, 'StepTolerance', 1e-8); vopt = fmincon(cost, v0, [], [], [], [], lb, ub, ... @(v) mycon(v, N, dt), opts); % 画图：上边是状态轨迹，下边是最优控制 t = (0:N-1) * dt; figure subplot(2,1,1); plot(t, vopt(1:N)); grid on; title('状态轨迹 x(t)'); subplot(2,1,2); plot(t, vopt(N+1:end)); grid on; title('最优控制 u(t)'); end function [c, ceq] = mycon(v, N, dt) c = []; % 不等式约束已放到 lb/ub 里，这里不用再写 ceq = zeros(N, 1); % 离散状态方程：x(k+1) = x(k) + (x(k)+u(k))*dt for k = 1:N-1 ceq(k) = v(k+1) - (v(k) + (v(k) + v(N+k)) * dt); end ceq(N) = v(N) - 1; % 末端状态约束：x(T)=1 end

逻辑说明：状态变量和控制变量被拼成一个大的设计向量v。目标函数只惩罚控制能量，状态本身的轨迹完全由约束方程决定。fmincon 的ceq返回一组等式约束，既有状态方程，又有末端条件，求解器会在每次迭代里同时满足它们。lb/ub直接对应执行器限幅，这是 bvp4c 处理不了的部分。

参数说明：N 取 40 是一个平衡点，网格太粗时矩形法离散误差会让结果偏离真解，太细时 fmincon 收敛变慢且更依赖初值。MaxIterations给到 500 一般够，遇到复杂模型可以提高；StepTolerance设 1e-8 是防止迭代提前停在精度不足的平台上。fmincon 默认算法 interior-point，对这类中等规模问题表现稳定。

4.3 选型判断：解析变分、数值边值、参数优化怎么选

拿到一个实际问题，我的选择顺序是这样的：先看约束是不是纯等式、有没有不等式限幅。如果只是两端固定、求最优曲线，优先走第 2 章的符号推导加第 3 章的 bvp4c；一旦出现控制输入限幅、路径避障这类不等式约束，直接跳到 fmincon 离散化。如果系统状态维数高，比如四轴飞行器十几维状态，bvp4c 的初值猜测基本猜不动，这时也认命用配置法。

这里给个后悔药：fmincon 解完，如果想确认离散误差，可以把 N 翻倍再算一次，两个最优目标值差在百分之一以内就可以接受。差太多就说明原来的网格不够密，不是什么模型出了问题。

5. 变分法上手的 5 个翻车现场：现象、原因与解决

5.1 functionalDerivative 求导结果里永远带着未知的 diff(y(x),x)

现象：functionalDerivative不报错，但返回的表达式没法化简，始终含有一堆diff(y(x), x)散件，看不出方程的样子。

原因：最常见的是把syms x y(x)写成了syms x y，后者把 y 声明成了普通符号变量而不是函数，导致diff(y,x)只能象征性地挂在表达式里。

解决：声明函数型变量必须是syms x y(x)。如果声明没问题，检查 F 里有没有把diff(y,x)写成别的符号别名，比如dy，这样符号引擎无法把它识别为 y 的导数。再一个容易被忽略的坑是，F 写成.^2还是^2都行，但 F 里一旦混入数值数组的.*就会在后续化简时报类型错。保持 F 是纯符号表达式，不要混入任何数值向量。

5.2 bvp4c 第一次迭代就 NaN，或者解直接飞到了 1e10

现象：最速降线的 ODE 很简单，但 bvp4c 就是跑不出来，返回 NaN，或者画出来的曲线离谱。

原因：ODE 右侧是-(1+y'^2)/(2y)，起点边界条件直接写成y(0)=0，分母为 0，第一次求导就爆掉。这属于模型的解析延拓问题，不是 bvp4c 的锅。

解决：数值上把起点抬到1e-3，边界条件同步写成ya(1)-1e-3，解和真实曲线的差别肉眼不可见。更稳妥的做法是初始猜测不要给常数，给一个接近解析解的曲线，比如最速降线在起点附近满足y ~ x^(2/3)：

function yinit = guess(x) yinit = [(x + 1e-3)^(2/3); (2/3)*(x+1e-3)^(-1/3)]; end solinit = bvpinit(linspace(0,1,50), @guess);

这个改进同时治好了「起点处斜率无穷」的数值病态问题。以后遇到 ODE 里带1/y或1/(y-a)这类项，都先想想边界点是不是碰上了奇异位置。

5.3 悬链线符号解 dsolve 给了一长串看不懂的表达式

现象：dsolve解悬链线 EL 方程，输出一个复杂到无法辨认的表达式，既不是 cosh 也不是 sinh，代值也验证不了。

原因：dsolve求的是通解，默认把所有数学分支都保留了，不化简成物理上有意义的形式。悬链线方程本身有多余的参数自由度，不给全边界条件时解的形式可以长得非常难看。

解决：先接受「符号通解验证」，用subs把候选的 cosh 形式代回 EL 验证，而不是指望dsolve直接吐出一个漂亮的 cosh。真正建模时直接用vpasolve对着端点条件求a、b的数值，符号形式反而容易把人误导进公式堆里。

5.4 fmincon 的优化结果对初值极其敏感，换个 v0 换个轨迹

现象：同样的模型，把v0从「零控制」改成「常数控制 0.1」，解出来的最优轨迹形状完全不同。

原因：优化问题非凸，fmincon 是局部优化器，不同初值落在不同的局部极值上。这在配点法里非常普遍，尤其是状态方程非线性时。

解决：先跑一个粗网格（N=10）拿到一个理性的初值，再把这个解插值到细网格上作为新的v0，这叫延续法。或者用多起点，固定随机种子生成 20 组v0，取目标值最小的结果。GlobalSearch函数也能干这事，但矩阵规模大了很费时间，建模阶段手动多起点更可控。

5.5 目标函数里 N=40，约束函数里 N=30，结果怎么都对不上

现象：fmincon 迭代正常结束，但画出来的状态轨迹在末端突然跳变，末端约束检查也不为零。

原因：把 N 和 dt 直接写死在了目标函数和约束函数里，改了一处漏了一处。更隐蔽的是，矩形法离散在 N 小时误差偏大，末端约束表面满足，但插值到连续时间时对不上。

解决：N 和 dt 一律通过函数参数传入，不要在工作区和子函数里重复写死。我的习惯是封装成一个config结构体，字段名带单位，比如config.N = 40; config.dt = config.T / config.N;，所有子函数只读这个结构体。离散误差问题则用第 6 章的网格加密检查来兜底，两层防线叠在一起基本不会翻车。

6. 让数值解经得起追问：变分解的验证与参数收敛检查

bvp4c 给出一条曲线，fmincon 给出一组数，这些都只是中间产物。评审、答辩或者论文审稿人会问一个问题：这个解对不对？我的验证套路固定分三步，一步都省不了。

第一步，把解代回原泛函算目标值。最速降线解出来以后，用trapz(sol.x, sqrt((1+sol.y(2,:).^2)./sol.y(1,:)))重新算一遍积分，看结果是否符合物理量级。如果画出来的曲线看着正常但积分值是负的，多半是 y 的正负方向定义错了。

第二步，网格加密做收敛性检查。写一个循环，N 从 20、50、100 到 200 各跑一遍，把每次的目标值画在semilogx图上。理想情况下目标值会随着 N 增大而单调逼近一个水平平台。如果在加密后仍剧烈摆动，说明问题出在边界条件或离散格式上，而不是网格密度。

Ns = [20 50 100 200]; Jvals = zeros(size(Ns)); for i = 1:numel(Ns) sol = run_brachistochrone(Ns(i)); % 封装好的求解函数 F = sqrt((1 + sol.y(2,:).^2) ./ sol.y(1,:)); Jvals(i) = trapz(sol.x, F); end semilogx(Ns, Jvals, '-o')

第三步，找特解对拍。悬链线对拍 cosh 解析解，最短路径对拍直线。特解对拍能一票否决符号推导里的漏项，比任何数值误差分析都直观。我的血泪经验是：第一次跑最速降线时，J值在网格加密到 100 后仍然缓慢漂移，最后发现是起点边界抬升1e-3带来的误差被AbsTol放大了。从那以后我把收敛曲线作为默认输出，不再只看单条轨迹。

这套验证流程同样适用于 fmincon 的配点法。优化目标值随网格收敛、与低阶特解对拍、约束残差检查，三件事做完，变分法在 Matlab 里的流程才算真正闭环。希望帮到你。

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