简介一份面向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)才能被识别为yfunctionalDerivative才能返回真正意义上的欧拉-拉格朗日表达式。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)嵌套项不化简根本没法看出y0这样的物理意义。上面这段代码的结果会化简成2*y 0也就是「最短路径是直线」符合直觉。换一个被积函数只要把F改掉就行。比如把 F 换成y * sqrt(1 diff(y,x)^2)求解结果会直接指向悬链线方程。整个脚本的可复用性就在这里以后做建模只需要维护一个F的表达式符号推导是一次性的。2.3 悬链线建模把符号解落到具体参数悬链线这个例子很适合验证符号推导的正确性给定两个悬挂点求一条处处受力平衡的曲线泛函写成势能积分F 取y*sqrt(1y^2)y 方向取向上为正。求出 EL 方程后可以手工把候选解代回去验证是不是y a*cosh((x-b)/a)。% 悬链线用符号工具箱把解析形式 ya*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]); % 初值猜 a1, b1这里容易踩的坑是vpasolve的初值非常影响结果换一组端点就要重新试初值。我的经验是先观察两个端点连线的中点位置把b初值设成中点横坐标a初值设成一个与「下垂程度」同量级的数大概率一次收敛。悬链线的符号验证和参数求解做完变分法的第一印象就建立起来了先由 EL 拿到方程再从方程里解出参数这是变分建模的标准节奏。3. 用 Matlab 数值求解变分问题以最速降线为例走通 bvp4c3.1 什么时候必须放弃解析解转向数值不是所有变分问题都有漂亮解析解。最速降线的 F 里带着根号和除法EL 方程是非线性二阶 ODE虽然能通过参数化技巧求到解析解但那个过程的推导量已经不适合建模竞赛里现场手推。更别提 F 里一旦混入摩擦、控制输入、分段约束解析这条路基本就走死了。这时可以把 EL 方程当成一个边值问题来处理二阶 ODE 加上两个端点条件数值求解。Matlab 内置的bvp4c就是干这个的。它的内部原理是配置法加牛顿迭代对大多数光滑边值问题表现稳定不需要自己写打靶循环。但它的输入参数比ode45多一层很多人第一次用会卡在「初始猜测怎么给」上。下面用一个完整例子把整个调用链走一遍。3.2 打靶法与 bvp4c 求解边值问题的完整脚本最速降线问题的泛函是J[y] ∫ sqrt((1y^2) / (2g y)) dx其中 y 向下为正g 是重力加速度。常数2g不影响极值点可以合并进 F 的常数因子。由 EL 方程可以化简出一个比较紧凑的二阶 ODEy −(1 y^2) / (2y)这个形式可以直接交给bvp4c。注意边界条件起点y(0)0会让 ODE 右侧除以 0所以在数值上必须把起点抬到1e-3量级否则第一轮迭代就 NaN。% 最速降线最小化 J ∫ sqrt((1y^2)/y) dxy 向下为正 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那样「给个初值沿着时间推」它是在整个区间上反复迭代修正解所以那几个控制参数对成败的影响很大。下面这张表是我常用的设置。参数默认值我常用的值作用RelTol1e-31e-4 ~ 1e-6控制解的相对误差建模竞赛里 1e-4 够用AbsTol1e-61e-6 ~ 1e-8控制接近零的状态分量的绝对精度NMax默认较小3000 ~ 10000允许的最大网格点数防止迭代中途提前终止bvpinit网格点数无20 起步100 收尾影响解的分辨率与收敛速度失败时最常见的警告是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(N1: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(N1:end)); grid on; title(最优控制 u(t)); end function [c, ceq] mycon(v, N, dt) c []; % 不等式约束已放到 lb/ub 里这里不用再写 ceq zeros(N, 1); % 离散状态方程x(k1) x(k) (x(k)u(k))*dt for k 1:N-1 ceq(k) v(k1) - (v(k) (v(k) v(Nk)) * 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 的初值猜测基本猜不动这时也认命用配置法。问题特征推荐工具原因纯泛函极值边界干净符号推导 bvp4c精度高不引入离散误差泛函复杂但无强约束bvp4c避免手推 ODE有控制限幅/避障约束fmincon 离散化不等式约束天然支持高维状态10fmincon 多起点打靶法初值没法猜这里给个后悔药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 右侧是-(1y^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)*(x1e-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 是局部优化器不同初值落在不同的局部极值上。这在配点法里非常普遍尤其是状态方程非线性时。解决先跑一个粗网格N10拿到一个理性的初值再把这个解插值到细网格上作为新的v0这叫延续法。或者用多起点固定随机种子生成 20 组v0取目标值最小的结果。GlobalSearch函数也能干这事但矩阵规模大了很费时间建模阶段手动多起点更可控。5.5 目标函数里 N40约束函数里 N30结果怎么都对不上现象fmincon 迭代正常结束但画出来的状态轨迹在末端突然跳变末端约束检查也不为零。原因把 N 和 dt 直接写死在了目标函数和约束函数里改了一处漏了一处。更隐蔽的是矩形法离散在 N 小时误差偏大末端约束表面满足但插值到连续时间时对不上。解决N 和 dt 一律通过函数参数传入不要在工作区和子函数里重复写死。我的习惯是封装成一个config结构体字段名带单位比如config.N 40; config.dt config.T / config.N;所有子函数只读这个结构体。离散误差问题则用第 6 章的网格加密检查来兜底两层防线叠在一起基本不会翻车。6. 让数值解经得起追问变分解的验证与参数收敛检查bvp4c 给出一条曲线fmincon 给出一组数这些都只是中间产物。评审、答辩或者论文审稿人会问一个问题这个解对不对我的验证套路固定分三步一步都省不了。第一步把解代回原泛函算目标值。最速降线解出来以后用trapz(sol.x, sqrt((1sol.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 里的流程才算真正闭环。希望帮到你。本文还有配套的精品资源点击获取 
长河君 · 2026/10/29 12:11:55 · 阅读 1.2k Matlab变分法建模实战：从欧拉-拉格朗日方程到最优控制
简介一份面向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)才能被识别为yfunctionalDerivative才能返回真正意义上的欧拉-拉格朗日表达式。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)嵌套项不化简根本没法看出y0这样的物理意义。上面这段代码的结果会化简成2*y 0也就是「最短路径是直线」符合直觉。换一个被积函数只要把F改掉就行。比如把 F 换成y * sqrt(1 diff(y,x)^2)求解结果会直接指向悬链线方程。整个脚本的可复用性就在这里以后做建模只需要维护一个F的表达式符号推导是一次性的。2.3 悬链线建模把符号解落到具体参数悬链线这个例子很适合验证符号推导的正确性给定两个悬挂点求一条处处受力平衡的曲线泛函写成势能积分F 取y*sqrt(1y^2)y 方向取向上为正。求出 EL 方程后可以手工把候选解代回去验证是不是y a*cosh((x-b)/a)。% 悬链线用符号工具箱把解析形式 ya*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]); % 初值猜 a1, b1这里容易踩的坑是vpasolve的初值非常影响结果换一组端点就要重新试初值。我的经验是先观察两个端点连线的中点位置把b初值设成中点横坐标a初值设成一个与「下垂程度」同量级的数大概率一次收敛。悬链线的符号验证和参数求解做完变分法的第一印象就建立起来了先由 EL 拿到方程再从方程里解出参数这是变分建模的标准节奏。3. 用 Matlab 数值求解变分问题以最速降线为例走通 bvp4c3.1 什么时候必须放弃解析解转向数值不是所有变分问题都有漂亮解析解。最速降线的 F 里带着根号和除法EL 方程是非线性二阶 ODE虽然能通过参数化技巧求到解析解但那个过程的推导量已经不适合建模竞赛里现场手推。更别提 F 里一旦混入摩擦、控制输入、分段约束解析这条路基本就走死了。这时可以把 EL 方程当成一个边值问题来处理二阶 ODE 加上两个端点条件数值求解。Matlab 内置的bvp4c就是干这个的。它的内部原理是配置法加牛顿迭代对大多数光滑边值问题表现稳定不需要自己写打靶循环。但它的输入参数比ode45多一层很多人第一次用会卡在「初始猜测怎么给」上。下面用一个完整例子把整个调用链走一遍。3.2 打靶法与 bvp4c 求解边值问题的完整脚本最速降线问题的泛函是J[y] ∫ sqrt((1y^2) / (2g y)) dx其中 y 向下为正g 是重力加速度。常数2g不影响极值点可以合并进 F 的常数因子。由 EL 方程可以化简出一个比较紧凑的二阶 ODEy −(1 y^2) / (2y)这个形式可以直接交给bvp4c。注意边界条件起点y(0)0会让 ODE 右侧除以 0所以在数值上必须把起点抬到1e-3量级否则第一轮迭代就 NaN。% 最速降线最小化 J ∫ sqrt((1y^2)/y) dxy 向下为正 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那样「给个初值沿着时间推」它是在整个区间上反复迭代修正解所以那几个控制参数对成败的影响很大。下面这张表是我常用的设置。参数默认值我常用的值作用RelTol1e-31e-4 ~ 1e-6控制解的相对误差建模竞赛里 1e-4 够用AbsTol1e-61e-6 ~ 1e-8控制接近零的状态分量的绝对精度NMax默认较小3000 ~ 10000允许的最大网格点数防止迭代中途提前终止bvpinit网格点数无20 起步100 收尾影响解的分辨率与收敛速度失败时最常见的警告是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(N1: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(N1:end)); grid on; title(最优控制 u(t)); end function [c, ceq] mycon(v, N, dt) c []; % 不等式约束已放到 lb/ub 里这里不用再写 ceq zeros(N, 1); % 离散状态方程x(k1) x(k) (x(k)u(k))*dt for k 1:N-1 ceq(k) v(k1) - (v(k) (v(k) v(Nk)) * 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 的初值猜测基本猜不动这时也认命用配置法。问题特征推荐工具原因纯泛函极值边界干净符号推导 bvp4c精度高不引入离散误差泛函复杂但无强约束bvp4c避免手推 ODE有控制限幅/避障约束fmincon 离散化不等式约束天然支持高维状态10fmincon 多起点打靶法初值没法猜这里给个后悔药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 右侧是-(1y^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)*(x1e-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 是局部优化器不同初值落在不同的局部极值上。这在配点法里非常普遍尤其是状态方程非线性时。解决先跑一个粗网格N10拿到一个理性的初值再把这个解插值到细网格上作为新的v0这叫延续法。或者用多起点固定随机种子生成 20 组v0取目标值最小的结果。GlobalSearch函数也能干这事但矩阵规模大了很费时间建模阶段手动多起点更可控。5.5 目标函数里 N40约束函数里 N30结果怎么都对不上现象fmincon 迭代正常结束但画出来的状态轨迹在末端突然跳变末端约束检查也不为零。原因把 N 和 dt 直接写死在了目标函数和约束函数里改了一处漏了一处。更隐蔽的是矩形法离散在 N 小时误差偏大末端约束表面满足但插值到连续时间时对不上。解决N 和 dt 一律通过函数参数传入不要在工作区和子函数里重复写死。我的习惯是封装成一个config结构体字段名带单位比如config.N 40; config.dt config.T / config.N;所有子函数只读这个结构体。离散误差问题则用第 6 章的网格加密检查来兜底两层防线叠在一起基本不会翻车。6. 让数值解经得起追问变分解的验证与参数收敛检查bvp4c 给出一条曲线fmincon 给出一组数这些都只是中间产物。评审、答辩或者论文审稿人会问一个问题这个解对不对我的验证套路固定分三步一步都省不了。第一步把解代回原泛函算目标值。最速降线解出来以后用trapz(sol.x, sqrt((1sol.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 里的流程才算真正闭环。希望帮到你。本文还有配套的精品资源点击获取