Matlab变分法实战：从欧拉-拉格朗日方程到最速降线数值求解

简介一份围绕变分法展开的Matlab建模教程文档聚焦最速降线、悬链线等经典问题面向正在学习数学建模、最优化理论或物理工程优化的读者适合用于掌握泛函极值的基本思想与Euler-Lagrange方程的求法。文中从变分法的历史背景讲起逐步深入到泛函、变分、极值条件等核心概念并以曲线长度最短、最速降线为例展示Euler-Lagrange方程的建立与求解过程为后续使用Matlab进行数值求解与可视化打下基础。资源包共包含1个doc格式文档整体大小617KB内容结构清晰包含完整的定义、推导与示例叙述前后衔接自然。目前已有138人学习下载适合需要在Matlab中使用变分法解决路径规划、能量最小化等实际问题的中高级学习者参考。1. 变分法和 Matlab 建模从“求数”到“求函数”大多数人用 Matlab 是在找数解方程求根、优化找极值、拟合出系数输入输出都是向量或矩阵。变分法恰恰反过来它求的是一个函数整条曲线、整条轨迹、整个场的分布。同样是“极值”两个字普通优化站在固定点上算方向变分法站在函数空间里找形状。这正是数学建模里一类高频问题的解法已知某个物理量依赖函数形状下落时间、弯曲能量、系统泛函让这个量取最小的那个形状是什么。这篇 Matlab 建模教程从欧拉-拉格朗日方程讲起用符号工具箱、bvp4c 求解器和优化工具箱把三条落地路线走通最后用最速降线案例完整验证一遍适合刚接触变分法但想直接出数值结果的建模参赛者也适合轨迹规划、图像分割里需要泛函极值的老手。2. 变分法核心泛函极值与欧拉-拉格朗日方程的推导路径2.1 泛函极值的变分形式为什么不是普通求导普通函数取极值看的是自变量变化时函数值怎么变。变分问题里自变量本身是一个函数 $y(x)$目标量写成一个积分$$ J[y] \int_a^b F(x, y, y),dx $$$J$ 不再是一个数到数的映射而是一个“函数到数”的映射这样的映射叫泛函。要问 $J$ 在哪个 $y(x)$ 上取极小不能直接求导而是要让 $y(x)$ 整体动一下$y(x) \to y(x) \varepsilon \eta(x)$。其中 $\eta(x)$ 是区间内任意一个光滑小扰动且在两个端点满足 $\eta(a)\eta(b)0$。这个“整体挪动”就是变分取极值的必要条件是对任意 $\eta$ 都有 $\left.\frac{d}{d\varepsilon}J[y\varepsilon\eta]\right|_{\varepsilon0}0$。这个条件看起来抽象但拆开就是一次分部积分的事。把 $F$ 展开到一阶$$ \frac{d}{d\varepsilon}J \int_a^b \left( F_y \eta F_{y}\eta \right) dx $$对第二项分部积分并把端点条件带入就能把 $\eta$ 换成 $\eta$$$ \int_a^b \eta \left( F_y - \frac{d}{dx}F_{y} \right)dx 0 $$因为 $\eta$ 可以任意取括号里的表达式必须在每个点上等于 0这就得到了核心方程。2.2 欧拉-拉格朗日方程的手推与必背结构从上一步直接得到$$ \frac{\partial F}{\partial y} - \frac{d}{dx}\frac{\partial F}{\partial y} 0 $$这就是欧拉-拉格朗日方程。注意这里的 $\frac{d}{dx}$ 是沿着解曲线的全导数展开是$$ \frac{d}{dx}F_{y} F_{yy}y F_{yy}y $$所以它本质上是一个二阶常微分方程。建模时拿到一个变分问题第一步就是把这个方程写出来写出来以后问题就变成了“解 ODE 边值问题”Matlab 里那一整套求解器就可以上场了。有一类常见情况值得单独背下来当被积函数 $F$ 不显含 $x$ 时欧拉-拉格朗日方程可以积分一次得到 Beltrami 恒等式$$ F - yF_{y} C $$这个式子把二阶问题降成一阶很多经典变分问题最速降线、最小旋转曲面都是靠它直接写出解析解的。判断 $F$ 是否显含 $x$在建模论文里是写清楚推导过程的关键一步数值求解时也能省掉一个未知函数。2.3 把变分问题翻译成 Matlab 能解的两种形式拿到一个具体建模题下一步要想清楚交给 Matlab 哪种形式。我把常见做法归纳成两类第一类把欧拉-拉格朗日方程整理成二阶 ODE再写成两个一阶 ODE 的标准形式 $y_1y_2,\ y_2g(x,y_1,y_2)$交给bvp4c或bvp5c求解。这是数值上最稳定、最通用的方式适合边界条件明确、端点固定的问题比如悬链线、测地线、弹性梁的挠曲。第二类干脆不解微分方程直接把泛函离散化。把 $y(x)$ 表示成一组基函数或一组网格节点值把积分 $J$ 写成离散求和然后调用优化工具箱里的fminunc或fmincon把“求函数极值”变成“求参数极值”。这种做法在边界复杂、或者除了变分之外还带不等式约束时特别好用比如带最大曲率约束的路径规划。下表是几个典型变分模型的 $F$ 和对应的欧拉-拉格朗日结果建模时可以直接对照选用。问题被积函数 $F$欧拉-拉格朗日方程最短路径$\sqrt{1y^2}$$y0$悬链线$y\sqrt{1y^2}$$yy1y^2$最速降线$\sqrt{(1y^2)/y}$$2yy1y^20$最小旋转曲面$y\sqrt{1y^2}$$yy1y^2$同悬链线注意悬链线和最小旋转曲面的 $F$ 长得一样但边界条件和目标物理含义不同解虽然都是悬链线参数选择完全不同。这就是为什么建模时不能只背公式要把边界条件写清楚。3. Matlab 实现变分法三条可落地的数值路线3.1 路线一符号工具箱导出欧拉-拉格朗日方程手推欧拉-拉格朗日方程容易在 $\frac{d}{dx}F_{y}$ 的链式法则上出错尤其是 $F$ 同时依赖 $y$、$y$ 的时候。我一般的做法是先用符号计算把方程导出来确认结构再写数值求解器。以最速降线为例设 $y$ 是向下为正的下落深度$F\sqrt{(1y^2)/y}$syms x y dy d2y F sqrt((1 dy^2) / y); % 最速降线的被积函数去掉了常数 1/2g % 欧拉-拉格朗日F_y - d/dx(F_dy) 0 Fy diff(F, y); Fdy diff(F, dy); % d/dx(F_dy) F_dyy * y F_dydy * y dFdy_dx diff(Fdy, y) * dy diff(Fdy, dy) * d2y; EL simplify(Fy - dFdy_dx)这段代码的关键在于dFdy_dx的展开。diff(Fdy, y)和diff(Fdy, dy)只对符号变量做偏导要还原全导数必须自己乘上 $y$ 和 $y$这是符号推导最容易漏的一步。运行后EL化简为 $- (1 dy^2 2yd2y)/(2*y^2)$令分子为零就得到 $2yy 1 y^2 0$和查表结果一致。这个路线适合建模初期把方程形式定下来同时也能顺手检查自己有没有把物理量放错位置。符号推导针对的是“求出方程”后续真正跑数值还靠下面两条路线。3.2 路线二bvp4c 数值边值流程与悬链线示例拿到二阶 ODE 后第一步永远是把二阶降成一阶。以悬链线为例方程 $yy1y^2$令 $y_1y$, $y_2y$得到function dydx catenaryODE(x, y) dydx [y(2); (1 y(2)^2) / y(1)]; end边界条件取悬链线的经典场景绳子两端固定在 $(-1,\cosh 1)$ 和 $(1,\cosh 1)$function res catenaryBC(ya, yb) L 1; res [ya(1) - cosh(L); yb(1) - cosh(L)]; end solinit bvpinit(linspace(-1, 1, 50), catenaryGuess); opts bvpset(RelTol, 1e-6, AbsTol, 1e-8); sol bvp4c(catenaryODE, catenaryBC, solinit, opts); xq linspace(-1, 1, 200); yq deval(sol, xq, 1); % deval 的第三个参数 1 表示取第 1 个分量 plot(xq, yq, r-, LineWidth, 1.5); hold on; plot(xq, cosh(xq), k--); legend(bvp4c 数值解, 解析解 cosh(x), Location, north);其中catenaryGuess提供 $y_1\cosh x$、$y_2\sinh x$ 的初值猜测。bvpset里的RelTol和AbsTol控制残差精度对于边界值问题默认值有时不够建模比赛里我会至少压到 1e-6。deval负责在任意加密网格上插值画图前一定要用它把解取出来不能直接用sol里的稀疏网格画。3.3 路线三直接离散泛函用优化工具箱求解不解微分方程直接对泛函下手。把区间 $[-1,1]$ 分成 $N$ 段节点值 $u_i$ 就是未知量积分用中点公式离散N 60; h 2 / N; xq -1:h:1; uIni cosh(xq); % 初始猜测用解析解 % 离散泛函: sum h * u * sqrt(1 (du/dx)^2) costF (u) sum(h * sqrt(1 (diff(u) ./ h).^2) .* ... 0.5 .* (u(1:end-1) u(2:end))); % 固定端点: u(-1)cosh(1), u(1)cosh(1) Aeq zeros(2, N1); Aeq(1,1) 1; Aeq(2,N1) 1; beq [cosh(1); cosh(1)]; opts optimoptions(fmincon, Display, final, ... Algorithm, interior-point, MaxIterations, 2000); uOpt fmincon(costF, uIni, [], [], Aeq, beq, [], [], [], opts); plot(xq, uOpt, b-o, MarkerSize, 3); hold on; plot(xq, cosh(xq), k--);这里Aeq是等式约束矩阵第一行锁定左端点第二行锁定右端点。costF里的0.5*(u(1:end-1)u(2:end))是中点处的 $y$ 值配合diff(u)得到的中点导数正好是梯形法则的精度。跑出来数值解和解析解在 $N60$ 时已经画不出区别。这条路线的好处是边界条件里如果还有不等式约束比如曲率不能超过某值直接在fmincon的A、b参数里加线性约束或者在nonlcon里加非线性约束比改造 bvp4c 省事得多。代价是网格加密后变量数上涨收敛速度不如专门的边值求解器。4. 完整建模实战最速降线的求解、初值与奇异点排错4.1 最速降线的目标泛函和边界条件设定最速降线问题的标准设定是质点在重力作用下从原点沿无摩擦曲线滑到点 $(1,1)$这里 $y$ 向下为正问曲线形状使滑行时间最短。由能量守恒 $v\sqrt{2gy}$时间泛函为$$ T[y] \int_0^1 \frac{\sqrt{1y^2}}{\sqrt{2gy}},dx $$常数 $\sqrt{2g}$ 不影响极值点建模时直接去掉。注意 $y(0)0$ 处被积函数分母为 0泛函在起点奇异这会在数值求解里埋一个坑下面专门处理。4.2 解析解先行用 Beltrami 恒等式避开奇异点$F$ 不显含 $x$直接用 Beltrami 恒等式 $F - yF_{y}C$。代入 $F\sqrt{(1y^2)/y}$ 后整理得到一阶方程$$ y(1y^2) C_0 $$这个式子本身就能在论文里完成最核心的建模推导。再做参数化令 $y a(1-\cos\theta)$代入可得 $x a(\theta-\sin\theta)$这就是摆线。Matlab 里用fsolve把 $a$ 和 $\theta$ 定下来% 边界: x(终点)1, y(终点)1 % x a(theta - sin(theta)), y a(1 - cos(theta)) fun (a) a .* (acos(1 - 1/a) - sin(acos(1 - 1/a))) - 1; a fsolve(fun, 0.5); thetaEnd acos(1 - 1/a); theta linspace(0, thetaEnd, 200); xSol a .* (theta - sin(theta)); ySol a .* (1 - cos(theta)); plot(xSol, ySol, b-, LineWidth, 1.5);fsolve的初值选 0.5 比较安全因为终点 $(1,1)$ 对应的摆线参数通常在 0.3 到 1.5 之间。如果初始值给成负数或零acos会直接出错。4.3 数值解验证遇到的 singular Jacobian 该怎么处理拿到解析解后再用 bvp4c 做独立验证。前面已经推导出二阶方程 $2yy1y^20$写成标准形式function dydx brachODE(x, y) dydx [y(2); -(1 y(2)^2) / (2 * y(1))]; end function res brachBC(ya, yb) res [ya(1) - 1e-3^(2/3); yb(1) - 1]; end x0 1e-3; solinit bvpinit(linspace(x0, 1, 100), (x) [x.^(2/3); (2/3)*x.^(-1/3)]); opts bvpset(RelTol, 1e-4, AbsTol, 1e-6, Stats, on); try sol bvp4c(brachODE, brachBC, solinit, opts); xq linspace(x0, 1, 300); yNum deval(sol, xq, 1); plot(xq, yNum, r-, LineWidth, 1.5); hold on; plot(xSol, ySol, k--); catch ME disp(bvp4c 失败通常是奇异 Jacobian); rethrow(ME); end这里最关键的排错点在于起点。原问题 $y(0)0$ 直接作为边界会给求解器带来奇异 Jacobian因为方程右端除 $y(1)$ 项在零点爆掉。常见做法是把起点挪到 $x_010^{-3}$左端点边界取 $y(x_0)x_0^{2/3}$这是从摆线近端点的渐近行为 $y\sim x^{2/3}$ 推出来的。bvpinit的初值猜测函数也要匹配这个奇异性。给 $y_1 x^{2/3}$、$y_2 (2/3)x^{-1/3}$比起始猜测给成线性函数要稳得多。如果照抄 $y_1 x$ 这类常规初值bvp4c 大概率在第一次迭代就报 singular Jacobian这是最典型的失败路径。4.4 网格、容差和初值对结果的影响对比数值解和解析解的误差来源主要有三个起点偏移、网格密度、容差设置。我用三组参数跑对比起点 x0网格数RelTol终点处误差现象1e-2501e-32.1e-2明显偏离摆线1e-31001e-44.2e-4基本重合1e-42001e-63.1e-5与解析解无法区分如果发现曲线在终点附近出现抖动优先把RelTol调到 1e-5 以下其次是加密网格。起点偏移的影响是系统性的$x_0$ 越大数值解越往终点方向偏。实际建模时我会先按 $x_010^{-3}$ 跑一遍再缩到 $10^{-4}$ 对比一次两次结果一致才认为数值解可信。5. 收敛性验证与最优控制里的变分法延伸5.1 用解析解当标尺做收敛性验证最速降线最大的优势是有解析解这给数值实现提供了一个免费标尺。我通常会固定容差只加密网格记录数值解与解析解的最大偏差err max(abs(yNum - interp1(xSol, ySol, xq, pchip)));这个误差应该随网格数 $N$ 增大而下降如果反而上升或者振荡多半是初值猜测和真实解相差太远bvp4c 收敛到了错误的解上。对于没有解析解的变分问题一个替代做法是用 $N200$ 和 $N400$ 两套网格各解一次看最大偏差是否在容差量级这就是数学建模论文里常说的“网格无关性验证”。5.2 从变分法到最优控制LQR 视角变分法和最优控制是同一个框架的两副面孔。状态方程 $\dot{x} Ax Bu$性能指标$$ J \int_0^T (x^TQx u^TRu),dt $$的极值条件就是欧拉-拉格朗日方程在向量情形下的推广配合协态变量得到的就是最优控制里的 Pontryagin 最小值原理。建模竞赛里出现轨迹规划题时最省力的路径是把控制变量消掉写出 Hamilton 方程用 bvp4c 解状态和协态的两点边值问题。这与本文第 3、4 章的做法完全一致只是 $y$ 从标量变成向量边界条件从两个点变成状态终点约束。5.3 一个实用技巧用参数扫描的上一组解做初值bvp4c 对初值敏感尤其是带非线性项或奇异项的问题。一个很少写进教程但非常有效的技巧是先把某个参数设成容易收敛的值解出来再把解作为下一步附近参数问题的bvpinit。相同的思路用来处理最速降线的端点移动很直观——把终点横坐标从 0.5 逐步扫到 1每一步都用上一步的解做初值基本不会翻车。% 参数续贯法示意 solPrev sol; for farX 0.6:0.1:1.0 % 更新边界条件中的 farX solinit bvpinit(solPrev.x, solPrev.y); solPrev bvp4c(brachODE, (ya, yb) [ya(1)-1e-3^(2/3); yb(1)farX], solinit); endbvpinit(solPrev.x, solPrev.y)直接把上一组解的整体形状传进去比重新猜一个解析函数快得多。对于新手来说与其纠结初值函数怎么写不如先把问题在“好解”的参数下跑通再一点点逼近目标参数大多数奇异 Jacobian 问题都能这样绕过去。本文还有配套的精品资源点击获取