5.11 曲线曲面几何与变分法初步

前一节处理积分的理论基础，本节回到几何与优化之间的交界处。前面的内容处理的都是「给定函数求积分或导数」，变分法处理的问题形式不同：在所有函数中找使某个泛函取极值的那一个。最短路径、最小表面积、最小能量、最速下降曲线都属于这一类，「求一个函数」比「求一个数」高一个层次，因此需要新的工具。

两类工具构成本节的主线。第一类是微分几何：曲线与曲面的局部形状由曲率、度量张量（第一基本形式）与曲率张量（第二基本形式）描述，这些量决定了「长度」「面积」「弯曲程度」这些几何量如何计算。第二类是变分原理：泛函的一阶变分为零给出 Euler-Lagrange 方程，它把无穷维的极值问题化为一个微分方程；反过来，物理与几何中的大量微分方程都可以写成变分形式（最小作用量原理、极小曲面方程、测地线方程）。两条线索的交点是测地线：曲面上长度最短的曲线，它的方程既可以用曲率刻画（测地曲率为零），也可以用变分法导出（度量给出的 Euler-Lagrange 方程）。

解析求解 Euler-Lagrange 方程通常只能处理少数情形（线性方程、可分离的方程、有守恒量的方程）。实际问题的求解依靠直接方法：在有限维子空间中寻找泛函的极小值，把问题化为线性代数或非线性规划。这条路线在现代计算中最成功的产物是有限元方法，它的理论基础是弱形式与 Sobolev 空间，本节会给出弱形式的来源、$H^1$ 空间的定义与一维有限元的完整实现（实验三的数据给出 $O(h^2)$ 的收敛率与节点值的超收敛现象）。

5.11.1 曲线、弧长与曲率

参数化与弧长

平面曲线通常用参数表示：$\gamma (t)=(x(t),y(t))$，$t\in [a,b]$。弧长是速度模长的积分：

$$L={\int }_a^b\Vert {\gamma }^{\prime }(t)\Vert dt={\int }_a^b\sqrt{x^{\prime }(t{)}^2+y^{\prime }(t{)}^2}dt$$

弧长是曲线的内在量，与参数化方式无关（换参数时速度与 $dt$ 的变化互相抵消）。弧长参数化用 $s$ 代替 $t$，使 $\parallel {\gamma }^{\prime }(s)\parallel \equiv 1$，这时的速度是单位切向量 $\tau (s)={\gamma }^{\prime }(s)$，几何量的公式最简洁。

曲率

曲线的曲率 $\kappa$ 度量它偏离直线的程度：单位切向量的转动速率，$\kappa =\parallel {\tau }^{\prime }(s)\parallel$。在一般参数下

$$\kappa =\frac{|x^{\prime }{y}^{″}-y^{\prime }{x}^{″}|}{{(x^{\prime 2}+y^{\prime 2})}^{3/2}}$$

分子是「速度与加速度的外积」（二维时的叉积），分母是速度模长的三次方，量纲上是长度分之一次方。$\kappa$ 的倒数 $R=\frac{1}{\kappa }$ 称曲率半径，它是该点处最佳逼近圆（密切圆）的半径。

三条基本结论。

直线的曲率恒为零（$x^{\prime }{y}^{″}-y^{\prime }{x}^{″}=0$），这是「曲率度量偏离直线的程度」的直接体现。

圆的曲率恒为半径的倒数 $\frac{1}{R}$。实验四的数据验证这一点：单位圆的数值曲率在 $[1.0000,1.0000]$。

曲率与参数化无关：换参数后分子分母同时变化，比值不变（这一点可以把上式写成 $\kappa =\parallel {\gamma }^{\prime }\times {\gamma }^{″}\parallel /\parallel {\gamma }^{\prime }{\parallel }^3$ 看清）。

Frenet 标架把曲线的局部几何完整描述下来。在弧长参数下定义切向量 $\tau$、法向量 $n$（$\tau$ 逆时针旋转 $\frac{\pi }{2}$），则

$$\frac{d\tau }{ds}=\kappa n,\frac{dn}{ds}=-\kappa \tau$$

这两条Frenet 公式说明曲率是「标架转动」的唯一速率。三维情形多一个挠率 ${\tau }_{\text{torsion}}$ 描述曲线偏离平面的程度，公式的形态类似。

例子：椭圆 ${\left(\frac{x}{a}\right)}^2+{\left(\frac{y}{b}\right)}^2=1$ 的曲率在长轴端点取最小值 $\frac{b}{a^2}$，在短轴端点取最大值 $\frac{a}{b^2}$。实验四的数据：$a=1.5$、$b=1$ 的椭圆曲率范围 $[0.4444,1.5000]$（理论值 $\frac{1}{2.25}=0.4444$ 与 $1.5$），$a=2$、$b=1$ 的范围 $[0.2500,2.0000]$（理论值 $0.25$ 与 $2$）。曲率的动态范围随偏心率增大而扩大，这条性质在图像处理与形状分析中用作形状特征。

5.11.2 曲面、基本形式与测地线

参数曲面与第一基本形式

曲面用两个参数表示：$r(u,v)=(x(u,v),y(u,v),z(u,v))$。两个偏导数 $r_u$、$r_v$ 张成切平面，其叉积给出法向量。第一基本形式是切平面上的内积：

$$I=Ed{u}^2+2Fdudv+Gd{v}^2,E=r_u\cdot r_u,F=r_u\cdot r_v,G=r_v\cdot r_v$$

它的系数矩阵 $\left(\begin{array}{cc}E& F\\ F& G\end{array}\right)$ 是切平面上的度量张量，决定曲面上长度的计算：曲线 $u(t),v(t)$ 的长度是 $\int \sqrt{E{\dot{u}}^2+2F\dot{u}\dot{v}+G{\dot{v}}^2}dt$。面积元是

$$dS=\sqrt{EG-F^2}dudv$$

这与 5.9 节的 $\parallel r_u\times r_v\parallel$ 是同一个量（由 Lagrange 恒等式 $\parallel r_u\times r_v{\parallel }^2=EG-F^2$）。

第二基本形式与曲率

第二基本形式度量曲面如何偏离切平面：

$$II=Ld{u}^2+2Mdudv+Nd{v}^2,L=r_{uu}\cdot n,M=r_{uv}\cdot n,N=r_{vv}\cdot n$$

其中 $n$ 是单位法向量。两个基本形式给出两类的曲率量。

法曲率：曲面沿切方向 $\dot{u}\dot{v}$ 的弯曲程度，${\kappa }_n=\frac{II}{I}$。

主曲率 ${\kappa }_1$、${\kappa }_2$ 是法曲率的极值，由广义特征值问题 $det(II-\kappa I)=0$ 给出，即

$$(EG-F^2){\kappa }^2-(EN-2FM+GL)\kappa +(LN-M^2)=0$$

高斯曲率 $K={\kappa }_1{\kappa }_2=\frac{LN-M^2}{EG-F^2}$，平均曲率 $H=\frac{{\kappa }_1+{\kappa }_2}{2}$。高斯曲率是曲面内在的（Gauss 绝妙定理）：它只依赖第一基本形式，可以用曲面上的长度与角度测量出来，不需要把曲面放在三维空间中观察。这条定理是微分几何最重要的结论之一，也是相对论中「时空弯曲的内禀性」这一说法的来源。

三类曲面的分类：$K>0$（如球面，椭球面）、$K=0$（如平面、圆柱面、可展曲面）、$K<0$（如马鞍面、双曲面）。半径 $R$ 的球面 $K=\frac{1}{R^2}$、$H=\frac{1}{R}$，是 $K$ 为正常数曲面中 $H$ 最小的（等周性质）。

测地线

曲面上两点之间的最短曲线称测地线。它有两条等价的刻画。

变分刻画：在连接两点的所有曲线中使长度最小，即长度泛函的变分取零。球面上的测地线是大圆（过两点的平面与球面的交线）。

几何刻画：测地线的「测地曲率」为零，即它的加速度向量没有切向分量（加速度只沿法向）:

$${\gamma }^{″}(s)\perp T_{p}S\text{（在弧长参数下）}$$

第二种刻画把测地线方程写成 ${\gamma }^{″}={\kappa }_{n}n$（$n$ 是曲面法向量），这是曲面上的「直线运动」的准确含义：一个沿曲面自由运动的质点，其加速度只沿法向。

平面上的测地线是直线（曲率为零），球面上的测地线是大圆（曲率等于球面半径的倒数），圆柱面（$K=0$）上的测地线展开后是直线（螺旋线）。这些例子说明「最短」与「最直」两条刻画给出同一批曲线，它们在弯曲空间中的角色与欧氏空间中的直线相同。

5.11.3 变分法与 Euler-Lagrange 方程

泛函与变分

泛函是把函数映到实数的映射，例如

$$J[y]={\int }_a^{b}F(x,y(x),y^{\prime }(x))dx$$

其中 $F$ 是给定的三元函数（称 Lagrange 函数或密度）。带边界条件 $y(a)=y_a$、$y(b)=y_b$ 的变分问题指：在上述边界条件的函数类中找使 $J$ 取极值的 $y$。

变分是把「函数的变化」形式化：取扰动 $\eta$（满足 $\eta (a)=\eta (b)=0$，保证扰动后仍满足边界条件），构造函数族 $y_{\epsilon }=y+\epsilon \eta$，把

$$\phi (\epsilon )=J[y+\epsilon \eta ]$$

看成 $\epsilon$ 的普通函数。极值点处 ${\phi }^{\prime }(0)=0$ 对任意 $\eta$ 成立——这一条就是「一阶变分为零」的定义，它与多元微积分中「梯度为零」的地位相同。

Euler-Lagrange 方程

计算 ${\phi }^{\prime }(0)$：

$$\delta J\text{defeq}{\phi }^{\prime }(0)={\int }_a^b[\frac{\partial F}{\partial y}\eta +\frac{\partial F}{\partial y^{\prime }}{\eta }^{\prime }]dx$$

第二项含 ${\eta }^{\prime }$，用分部积分把它转成 $\eta$ 的形式：

$${\int }_a^b\frac{\partial F}{\partial y^{\prime }}{\eta }^{\prime }dx={\left[\frac{\partial F}{\partial y^{\prime }}\eta \right]}_a^b-{\int }_a^b\frac{d}{dx}\left(\frac{\partial F}{\partial y^{\prime }}\right)\eta dx$$

边界项因 $\eta (a)=\eta (b)=0$ 而消失，因此

$$\delta J={\int }_a^b[\frac{\partial F}{\partial y}-\frac{d}{dx}\left(\frac{\partial F}{\partial y^{\prime }}\right)]\eta dx$$

由 $\eta$ 的任意性（变分法基本引理：若连续函数与所有这样的 $\eta$ 的积分都为零，则函数恒为零），得 Euler-Lagrange 方程

$$\frac{\partial F}{\partial y}-\frac{d}{dx}\left(\frac{\partial F}{\partial y^{\prime }}\right)=0$$

三个要点。第一，变分问题（无穷维）被化成一个微分方程（含二阶导数），这是变分法最核心的一步。第二，边界项消失用的是边界条件给出的约束（$\eta$ 在端点为零），若边界自由则给出自然边界条件 $\frac{\partial F}{\partial y^{\prime }}{|}_{a,b}=0$。第三，$F$ 不显含 $x$ 时存在守恒量（Beltrami 恒等式）：

$$F-y^{\prime }\frac{\partial F}{\partial y^{\prime }}=\text{常数}$$

这条恒等式在求最速降线、极小旋转曲面等问题的显式解时是关键工具。

三个经典例子

最短路径：$F=\sqrt{1+y^{\prime 2}}$，EL 方程给出 $\frac{d}{dx}\frac{y^{\prime }}{\sqrt{1+y^{\prime 2}}}=0$，即 $y' = $ 常数，解为直线。这条结果与欧氏几何的直觉一致，但它是从变分原理推出来的，说明「最短」这个几何性质可以由分析方法导出。

最速降线：求一条曲线使质点从 $(0,0)$ 无初速下滑到 $(x_1,y_1)$ 的时间最短。由能量守恒，速度 $v=\sqrt{2gy}$，因此

$$T[y]={\int }_0^{x_1}\frac{\sqrt{1+y^{\prime 2}}}{\sqrt{2gy}}dx$$

$F$ 不显含 $x$，用 Beltrami 恒等式得 $y(1+y^{\prime 2})=2a$（$a$ 为常数），解出参数形式的摆线

$$x=a(\theta -\sin \theta ),y=a(1-\cos \theta )$$

下滑时间 $T=\sqrt{\frac{a}{g}}{\theta }_1$。实验二的数据：取 $a=0.5$、${\theta }_1=\pi$，终点 $(1.5708,1.0)$，解析值 $T=0.7096134476$ 秒，数值积分给出 $0.7096134480$ 秒（相对偏差 $6\times {10}^{-10}$）；直线的下滑时间 $0.8412091759$ 秒（比摆线多 $18.5\mathrm{\%}$），起点斜率取直线斜率的 $0.7$ 倍与 $0.4$ 倍的抛物线分别是 $0.9192494657$ 秒与 $1.0628828407$ 秒。摆线在起点处切线竖直（斜率无穷），这一条是它优于其他曲线的原因之一：竖直的起始段让质点迅速获得速度。起点斜率趋于零的曲线（例如顶点在起点的抛物线）下滑时间对数发散，因此不属于可比的候选。

等周问题：在给定周长的闭曲线中找面积最大的一条。用弧长参数化并写成变分问题，EL 方程给出曲率恒为常数，即圆。实验四的数据用等周不等式的形式核对：

$$L^2\ge 4\pi A$$

单位圆取等号（比值 $1.000000$），椭圆 $a=1.5,b=1$ 的比值 $1.062657$、$a=2,b=1$ 的比值 $1.188827$，超椭圆 $|x{|}^4+|y{|}^4=1$ 的比值 $1.056872$，正方形 $1.273240$。数据说明两件事：不等式的等号只在圆处成立；面积的「效率」随形状偏离圆形而下降，椭圆的偏心率越大越不经济。

5.11.4 变分问题的直接方法与弱解

从经典解到弱解

Euler-Lagrange 方程要求解有二阶导数（方程中出现 $y^{″}$），而变分问题本身只要求 $J[y]$ 有意义（$y$ 与 $y^{\prime }$ 可积即可）。这中间的落差引出弱形式：把 EL 方程两边乘以试探函数 $\eta$ 并分部积分，得到

$${\int }_a^b(\frac{\partial F}{\partial y}\eta +\frac{\partial F}{\partial y^{\prime }}{\eta }^{\prime })dx=0,\mathrm{\forall }\eta \in H_0^1(a,b)$$

这个等式只含 $y$ 与 $\eta$ 的一阶导数，因此对更「粗糙」的函数也有意义。满足弱形式的 $y$ 称弱解，弱解在足够光滑时就是经典解。

线性情形：若 $F=\frac{1}{2}p(x)(y^{\prime }{)}^2-f(x)y$（$p>0$），弱形式成为

$${\int }_a^b(p{y}^{\prime }{\eta }^{\prime }-f\eta )dx=0,\mathrm{\forall }\eta \in H_0^1$$

对应的经典方程是 ${\left(p{y}^{\prime }\right)}^{\prime }+f=0$（一维的椭圆型方程）。取 $p=1$、$f=\sin \pi x$ 得 $-y^{″}=\sin \pi x$，解 $y=\frac{\sin \pi x}{{\pi }^2}$，这正是实验一与实验三处理的问题。

Sobolev 空间

弱形式需要一个新的函数空间：$H^1(a,b)$ 由「本身与弱导数都平方可积」的函数组成：

$$H^1(a,b)=\{u\in L^2(a,b):u^{\prime }\in L^2(a,b)\},\parallel u{\parallel }_{H^1}^2=\parallel u{\parallel }_{L^2}^2+\parallel u^{\prime }{\parallel }_{L^2}^2$$

其中 $u^{\prime }$ 是弱导数：若存在 $v\in L^2$ 使得对一切光滑且紧支撑的 $\phi$ 有 $\int u{\phi }^{\prime }=-\int v\phi$，则 $v$ 是 $u$ 的弱导数。弱导数的定义来自分部积分，它允许函数有拐角（分段线性函数的弱导数就是分段常数）。

$H^1$ 是完备的内积空间（与 5.10 节的 $L^p$ 完备性同源）；$H_0^1$ 是 $H^1$ 中在边界取零值的函数的闭包，它正是弱形式中试探函数所在的空间（对应边界条件 $\eta (a)=\eta (b)=0$）。

一维的嵌入定理：$H^1(a,b)\subset C[a,b]$，即 $H^1$ 函数可以取到连续的代表元。这条性质使一维问题的弱解自动连续，也让有限元的节点值有明确含义。高维情形没有这么好的性质（$H^1$ 在 ${\mathbb{R}}^2$ 上嵌入 $L^q$ 对任意有限 $q$，但仍可能有奇点），需要更细致的嵌入定理。

Ritz 方法与 Galerkin 方法

直接方法的核心思想是把无穷维的极小化问题限制在有限维子空间：取 $V_n=\mathrm{span}\{{\phi }_1,\dots ,{\phi }_n\}\subset H_0^1$（基函数），在 $V_n$ 中最小化 $J$。线性情形下极小化条件给出线性方程组

$$Ku=b,K_{ij}={\int }_a^{b}p{\phi }_j^{\prime }{\phi }_i^{\prime }dx,b_i={\int }_a^{b}f{\phi }_{i}dx$$

$K$ 称刚度矩阵，是对称正定的（$p>0$ 时）。这条推导就是 Ritz 方法；若直接从弱形式出发要求对每个 ${\phi }_i$ 成立（不经过极小化），得到的是 Galerkin 方法，在线性自伴问题中两者给出同一方程组。

Galerkin 正交性给出误差估计的基础：真解 $u$ 与近似解 $u_n$ 之差在能量内积下与整个 $V_n$ 正交，因此 $u_n$ 是 $u$ 在能量内积下的正交投影，$\parallel u-u_n{\parallel }_{\text{energy}}\le \parallel u-v{\parallel }_{\text{energy}}$ 对任意 $v\in V_n$ 成立。这条结论把误差估计化归为逼近误差，于是收敛率由 $V_n$ 的逼近能力决定。

有限元方法

有限元方法把基函数取为分片多项式：在网格上取分片线性（一维）或分片线性/双线性（高维）的「帽函数」，每个基函数在一个节点附近支集很小。三条优点：刚度矩阵稀疏（带状或分块带状），局部基函数使装配可以逐单元进行，且基函数的逼近性质与函数的二阶导数有关，收敛率可以证明。

一维分片线性的收敛率：设 $u\in H^2$，则能量范数误差 $O(h)$、$L^2$ 误差 $O(h^2)$、节点值误差 $O(h^2)$（在均匀网格上）。实验三的数据给出完整的验证：$L^2$ 误差从 $N=4$ 时的 $3.980\times {10}^{-3}$ 降到 $N=128$ 时的 $3.940\times {10}^{-6}$，每次网格加倍误差降到约四分之一（比值 $3.96$、$3.99$、$4.00$、$4.00$、$4.00$），对应 $O(h^2)$。

节点值的超收敛是本节数据的另一条信息：节点最大误差按每次网格加倍降到约 $\frac{1}{65}$ 的速度下降（$6.075\times {10}^{-8}\to 9.284\times {10}^{-10}\to 1.443\times {10}^{-11}\to 2.253\times {10}^{-13}$），远快于二阶。原因是一维问题的一个特殊性质：用分片线性基与精确积分装配刚度矩阵时，有限元解的节点值在所有节点上与精确解完全相等（这是「节点精确性」定理），这里观测到的残余误差来自载荷积分的数值求积（用的是三点 Gauss，误差 $O(h^6)$，与实测的比值 $2^6=64$ 相符）。$N=64$ 之后误差落到机器精度（$4.455\times {10}^{-15}$），比值失去意义。这条数据说明：误差的收敛阶受「最弱的一环」控制，有限元方法与求积公式的精度需要匹配，否则精度提升会在求积误差处停住。

5.11.5 变分方法的其他用途

变分表述的价值不限于求解经典问题，它提供了一种「把目标写成泛函再求极值」的建模方式。

极小曲面：面积泛函 $A=\iint \sqrt{1+z_x^2+z_y^2}dxdy$ 的 EL 方程是 $\mathrm{\nabla }\cdot \frac{\mathrm{\nabla }z}{\sqrt{1+|\mathrm{\nabla }z{|}^2}}=0$，解称极小曲面（肥皂膜的形状）。平均曲率为零是它的几何刻画（$H=0$）。

变分推断：把后验分布 $p(\theta \mid D)$ 用简单的分布族 $q$ 逼近，极小化 $\mathrm{KL}(q\parallel p)$；由 5.10 节的 Radon-Nikodym 定理，$q$ 有密度时 KL 是积分泛函，极小化问题就落在「在所有密度中找最优的一个」这个变分框架内。第四章 4.2 节讨论的 KL 散度与这里的变分表述是同一对象的两种视角。

正则化：把目标写成「数据项 + 正则项」的泛函，例如图像去噪中的总变差模型 $J[u]=\frac{1}{2}\parallel u-g{\parallel }^2+\lambda \int |\mathrm{\nabla }u|$。EL 方程是含 $\mathrm{\nabla }u/|\mathrm{\nabla }u|$ 的非线性方程，其弱形式与 Sobolev 空间的框架相同。

最优控制：状态方程约束下的目标极小化，Euler-Lagrange 方程推广为Pontryagin 极大值原理（引入协态变量，与 Lagrange 乘子同源）。这条路线在机器人轨迹规划、强化学习的最优控制表述中都出现。

守恒律：Noether 定理把泛函的对称性翻译为守恒量（时间平移对称给出能量守恒、空间平移给出动量守恒、旋转对称给出角动量守恒）。这条定理把「找守恒量」这件在物理中很困难的事化归为「找对称性」，是变分表述在理论物理中最有力的应用。

5.11.6 数值实验

实验一：Euler-Lagrange 方程与直接离散化

求解变分问题 $minJ[y]={\int }_0^1(y^{\prime 2}-2fy)dx$，$y(0)=y(1)=0$，$f(x)=\sin \pi x$。EL 方程是 $y^{″}=-f$，解为 $y=\frac{\sin \pi x}{{\pi }^2}$。直接离散化把 $J$ 用梯形公式近似，对节点值求导得线性方程组。

import math
import numpy as np

def u_exact(x):
    return np.sin(math.pi * x) / math.pi ** 2

print(f"{'N':>8} {'节点最大误差':>16} {'误差比值':>10} {'真 L2 误差':>16} {'L2 比值':>10}")
prev_i = prev_l = None
gp = np.array([-0.906179845938664, -0.538469310105683, 0.0,
               0.538469310105683, 0.906179845938664])
gw = np.array([0.236926885056189, 0.478628670499366, 0.568888888888889,
               0.478628670499366, 0.236926885056189])
for N in [4, 8, 16, 32, 64, 128]:
    h = 1.0 / N
    nodes = np.linspace(0, 1, N + 1)
    fv = np.sin(math.pi * nodes)
    A = np.zeros((N - 1, N - 1))
    for i in range(1, N):
        A[i - 1, i - 1] = 2.0 / h
        if i - 1 >= 1:
            A[i - 1, i - 2] = -1.0 / h
        if i + 1 <= N - 1:
            A[i - 1, i] = -1.0 / h
    b = h * fv[1:-1]      # 梯形法：端点处 y = 0，权重 1/2 的项为零
    u = np.zeros(N + 1)
    u[1:-1] = np.linalg.solve(A, b)
    err_i = float(np.max(np.abs(u - u_exact(nodes))))
    e2 = 0.0
    for e in range(N):                      # 真 L2 误差：逐单元 5 点 Gauss
        a0, b0 = nodes[e], nodes[e + 1]
        for q in range(5):
            xi = 0.5 * (a0 + b0) + 0.5 * h * gp[q]
            uh = u[e] + (u[e + 1] - u[e]) * (xi - a0) / h
            e2 += 0.5 * h * gw[q] * (uh - float(u_exact(xi))) ** 2
    err_l = math.sqrt(e2)
    ri = "-" if prev_i is None else f"{prev_i / err_i:.2f}"
    rl = "-" if prev_l is None else f"{prev_l / err_l:.2f}"
    print(f"{N:>8} {err_i:>16.3e} {ri:>10} {err_l:>16.3e} {rl:>10}")
    prev_i, prev_l = err_i, err_l

输出：

       N           节点最大误差       误差比值          真 L2 误差      L2 比值
       4        5.373e-03          -        1.770e-03          -
       8        1.312e-03       4.09        4.193e-04       4.22
      16        3.261e-04       4.02        1.034e-04       4.06
      32        8.142e-05       4.01        2.576e-05       4.01
      64        2.035e-05       4.00        6.436e-06       4.00
     128        5.086e-06       4.00        1.609e-06       4.00

两条误差列的比值都稳定在 $4.00$，即 $N$ 加倍时误差降到四分之一，对应 $O(h^2)$。这个阶数与理论一致：离散泛函用梯形公式逼近积分，误差 $O(h^2)$；离散化的 EL 方程（$-(y_{j-1}-2y_j+y_{j+1})/h=h{f}_j$）是二阶精度格式。节点误差与 $L^2$ 误差的常数因子不同（$N=4$ 时 $5.37\times {10}^{-3}$ 对 $1.77\times {10}^{-3}$，$L^2$ 误差是整体的平均意义，小于最大误差），但收敛阶相同。这条实验给出「解析路线」与「直接离散化路线」的一致性：两条路线给出的解在同一网格上收敛到同一个函数。

实验二：最速降线与摆线

比较摆线（EL 解）与其他曲线，并用直接离散化 + 优化求数值极小。

import math
import numpy as np
from scipy.optimize import minimize

g = 9.8
a = 0.5
x1 = a * (math.pi - math.sin(math.pi))
y1 = a * (1 - math.cos(math.pi))
T_cyc = math.sqrt(a / g) * math.pi
print(f"终点：x1 = {x1:.10f}, y1 = {y1:.10f}")
print(f"摆线（解析）：T = sqrt(a/g)*theta1 = {T_cyc:.10f} s")

def T_of(x, y, M=400001):
    """改用 y 作积分变量并令 y = u^2，消除 x = 0 处的奇性"""
    x = np.asarray(x, dtype=float)
    y = np.asarray(y, dtype=float)
    yg = np.linspace(0.0, y[-1], M)
    xg = np.interp(yg, y, x)
    dxdy = np.gradient(xg, yg)
    u = np.sqrt(yg)
    integrand = 2.0 * np.sqrt(1.0 + dxdy ** 2) / math.sqrt(2.0 * g)
    return float(np.trapezoid(integrand, u))

thn = np.linspace(0, math.pi, 400001)
print(f"摆线（数值积分）：T = {T_of(a * (thn - np.sin(thn)), a * (1 - np.cos(thn))):.10f} s")
xs = np.linspace(0, x1, 400001)
m = y1 / x1
for m0 in [m, 0.7 * m, 0.4 * m]:
    yv = m0 * xs + (y1 - m0 * x1) * (xs / x1) ** 2
    T = T_of(xs, yv)
    print(f"起始斜率 {m0 / m:.1f} 倍：T = {T:.10f} s，比摆线多 {(T / T_cyc - 1) * 100:.4f}%")

def T_element(yv, h):
    """分段线性曲线的下滑时间（单元内解析积分）"""
    yv = np.asarray(yv)
    d = np.diff(yv) / h
    K = (1.0 + d ** 2) / (2.0 * g)
    sa = np.sqrt(np.maximum(yv[:-1], 1e-300))
    sb = np.sqrt(np.maximum(yv[1:], 1e-300))
    ad = np.abs(d)
    seg = np.where(ad < 1e-12, h / sa, 2.0 * np.abs(sb - sa) / np.maximum(ad, 1e-300))
    return float(np.sum(np.sqrt(K) * seg))

def cycloid_y(xq):
    th = np.linspace(0, math.pi, 400001)
    return np.interp(xq, a * (th - np.sin(th)), a * (1 - np.cos(th)))

for N in [20, 40, 80, 160, 320]:
    h = x1 / N
    xg = np.linspace(0, x1, N + 1)
    y0 = 0.7 * xg + (y1 - 0.7 * x1) * (xg / x1) ** 2
    res = minimize(lambda z: T_element(np.concatenate([[0.0], z, [y1]]), h),
                   y0[1:-1], method="L-BFGS-B", bounds=[(1e-9, 3.0)] * (N - 1),
                   options={"maxiter": 200000, "maxfun": 200000, "ftol": 1e-20, "gtol": 1e-16})
    yfull = np.concatenate([[0.0], res.x, [y1]])
    print(f"N = {N:>4}: T = {T_element(yfull, h):.10f}，"
          f"与摆线的最大偏差 = {float(np.max(np.abs(yfull - cycloid_y(xg)))):.6f}")

输出：

终点：x1 = 1.5707963268, y1 = 1.0000000000
摆线（解析）：T = sqrt(a/g)*theta1 = 0.7096134476 s
摆线（数值积分）：T = 0.7096134480 s
起始斜率 1.0 倍：T = 0.8412091759 s，比摆线多 18.5447%
起始斜率 0.7 倍：T = 0.9192494657 s，比摆线多 29.5423%
起始斜率 0.4 倍：T = 1.0628828407 s，比摆线多 49.7834%
N =   20: T = 0.7126498323，与摆线的最大偏差 = 0.025965
N =   40: T = 0.7110892891，与摆线的最大偏差 = 0.016090
N =   80: T = 0.7103384053，与摆线的最大偏差 = 0.010024
N =  160: T = 0.7099719115，与摆线的最大偏差 = 0.006267
N =  320: T = 0.7097965189，与摆线的最大偏差 = 0.003466

第一段数据有三条信息。摆线的解析值与数值积分值相差 $4\times {10}^{-10}$（相对偏差 $6\times {10}^{-10}$），说明数值积分（把 $y$ 作为积分变量并令 $y=u^2$）准确地算出了 $\sqrt{a/g}\cdot {\theta }_1$ 这个解析结果；这一步是关键，因为直接用 $x$ 作积分变量时被积函数在起点处按 $x^{-1/2}$ 发散，梯形公式会给出系统性偏差。直线比摆线慢 $18.5447\mathrm{\%}$（$0.8412091759$ 对 $0.7096134476$ 秒），这个数可以作为「最速降线问题的最优性」的定量说明。起始斜率减半（$0.7$ 倍与 $0.4$ 倍）时下滑时间分别增加 $29.5\mathrm{\%}$ 与 $49.8\mathrm{\%}$，趋势单调：起始越平缓，质点获得速度越慢，总时间越长。极限情形（起始斜率为零）下滑时间对数发散。

第二段数据给出直接离散化的收敛行为。离散最小值从 $0.7126498323$（$N=20$）单调降到 $0.7097965189$（$N=320$），从上方逼近连续问题的极小值 $0.7096134476$，这是变分原理的普遍性质：限制在有限维子空间中的极小值不小于全局极小值。与摆线的最大节点偏差从 $0.025965$ 降到 $0.003466$，相邻比值 $1.61$、$1.61$、$1.60$、$1.81$ 对应 $2^{2/3}=1.587$，即偏差按 $O(h^{2/3})$ 下降。这个阶数远低于光滑问题的 $O(h^2)$，原因是摆线在起点处的形状：$x\approx \frac{a{\theta }^3}{6}$、$y\approx \frac{a{\theta }^2}{2}$ 给出 $y\sim C{x}^{2/3}$，切线竖直；分段线性函数逼近 $x^{2/3}$ 的误差只能达到 $O(h^{2/3})$。这条数据说明直接离散化的效率取决于解的光滑程度，遇到奇性时应当改用非均匀网格（在起点附近加密）或更高阶的单元。

实验三：一维有限元方法

用分片线性有限元求解 $-u^{″}=\sin \pi x$、$u(0)=u(1)=0$，观察节点误差与 $L^2$ 误差的收敛率。

import math
import numpy as np

def u_exact(x):
    return np.sin(math.pi * x) / math.pi ** 2

gp3 = np.array([-math.sqrt(3.0 / 5), 0.0, math.sqrt(3.0 / 5)])
gw3 = np.array([5.0 / 9, 8.0 / 9, 5.0 / 9])
gp5 = np.array([-0.906179845938664, -0.538469310105683, 0.0,
                0.538469310105683, 0.906179845938664])
gw5 = np.array([0.236926885056189, 0.478628670499366, 0.568888888888889,
                0.478628670499366, 0.236926885056189])
for N in [4, 8, 16, 32, 64, 128]:
    h = 1.0 / N
    nodes = np.linspace(0, 1, N + 1)
    A = np.zeros((N + 1, N + 1))
    b = np.zeros(N + 1)
    for e in range(N):                       # 逐单元装配刚度矩阵与载荷向量
        a0, b0 = nodes[e], nodes[e + 1]
        ke = (1.0 / h) * np.array([[1.0, -1.0], [-1.0, 1.0]])
        fe = np.zeros(2)
        for q in range(3):
            xi = 0.5 * (a0 + b0) + 0.5 * h * gp3[q]
            w = 0.5 * h * gw3[q]
            phi = np.array([(b0 - xi) / h, (xi - a0) / h])
            fe += w * math.sin(math.pi * xi) * phi
        for i in range(2):
            for j in range(2):
                A[e + i, e + j] += ke[i, j]
            b[e + i] += fe[i]
    free = np.arange(1, N)                   # 齐次 Dirichlet 边界条件
    u = np.zeros(N + 1)
    u[free] = np.linalg.solve(A[np.ix_(free, free)], b[free])
    err_i = float(np.max(np.abs(u - u_exact(nodes))))
    e2 = 0.0
    for e in range(N):
        a0, b0 = nodes[e], nodes[e + 1]
        for q in range(5):
            xi = 0.5 * (a0 + b0) + 0.5 * h * gp5[q]
            uh = u[e] + (u[e + 1] - u[e]) * (xi - a0) / h
            e2 += 0.5 * h * gw5[q] * (uh - float(u_exact(xi))) ** 2
    print(f"N = {N:>4}: 节点最大误差 = {err_i:.3e}, L2 误差 = {math.sqrt(e2):.3e}")

输出：

N =    4: 节点最大误差 = 6.075e-08, L2 误差 = 3.980e-03
N =    8: 节点最大误差 = 9.284e-10, L2 误差 = 1.005e-03
N =   16: 节点最大误差 = 1.443e-11, L2 误差 = 2.519e-04
N =   32: 节点最大误差 = 2.253e-13, L2 误差 = 6.302e-05
N =   64: 节点最大误差 = 4.455e-15, L2 误差 = 1.576e-05
N =  128: 节点最大误差 = 9.437e-16, L2 误差 = 3.940e-06

两列误差的收敛速度截然不同。$L^2$ 误差的比值是 $3.96$、$3.99$、$4.00$、$4.00$、$4.00$，稳定在 $4$，对应 $O(h^2)$，与有限元的理论收敛率一致。节点最大误差的比值是 $65.44$、$64.35$、$64.02$、$50.59$、$4.72$，前四行接近 $64=2^6$，即 $O(h^6)$，与 $L^2$ 误差相差三个数量级；最后两行的比值失去意义，因为误差已经落到机器精度（$4.455\times {10}^{-15}$ 与 $9.437\times {10}^{-16}$ 都在 ${10}^{-15}$ 量级，是双精度的舍入下限）。

节点值的超收敛可以用一句话解释：在这个一维问题中，若刚度矩阵与载荷向量都用精确积分装配，有限元解的节点值与精确解的节点值完全相等（节点精确性），因此观测到的节点误差只来自载荷积分的数值求积误差，而三点 Gauss 公式的误差是 $O(h^6)$。这条性质的实用含义是精度由最弱的一环决定：想利用节点超收敛，必须让求积误差低于目标误差；若求积只用中点或梯形（误差 $O(h^2)$），节点误差也会退化到 $O(h^2)$，超收敛消失。$L^2$ 误差不受这一环的影响，因为它是整体量，由逼近误差主导。

实验四：曲率与等周不等式

用周期中心差计算闭曲线的周长、面积、曲率范围与等周比值，并核对摆线的曲率公式。

import math
import numpy as np

def curve_stats(xt, yt, t):
    """闭曲线的周期中心差分（t 上首尾相接）"""
    h = t[1] - t[0]
    ux = (np.roll(xt, -1) - np.roll(xt, 1)) / (2 * h)
    uy = (np.roll(yt, -1) - np.roll(yt, 1)) / (2 * h)
    uxx = (np.roll(xt, -1) - 2 * xt + np.roll(xt, 1)) / h ** 2
    uyy = (np.roll(yt, -1) - 2 * yt + np.roll(yt, 1)) / h ** 2
    speed = np.sqrt(ux ** 2 + uy ** 2)
    kappa = np.abs(ux * uyy - uy * uxx) / np.maximum(speed ** 3, 1e-300)
    L = float(np.sum(speed) * h)
    A = 0.5 * float(np.sum(xt * uy - yt * ux) * h)
    return L, abs(A), kappa

t = np.linspace(0, 2 * math.pi, 200000, endpoint=False)
tests = [
    ("圆（半径 1）", np.cos(t), np.sin(t)),
    ("椭圆 a=1.5, b=1", 1.5 * np.cos(t), np.sin(t)),
    ("椭圆 a=2, b=1", 2.0 * np.cos(t), np.sin(t)),
    ("超椭圆 |x|^4+|y|^4=1",
     np.sign(np.cos(t)) * np.abs(np.cos(t)) ** 0.5,
     np.sign(np.sin(t)) * np.abs(np.sin(t)) ** 0.5),
]
for name, xt, yt in tests:
    L, A, k = curve_stats(xt, yt, t)
    print(f"{name:>22}: L = {L:.6f}, A = {A:.6f}, "
          f"L^2/(4 pi A) = {L ** 2 / (4 * math.pi * A):.6f}, "
          f"曲率 [{np.min(k):.4f}, {np.max(k):.4f}]")

a = 0.5
th = np.linspace(1e-3, math.pi - 1e-3, 400001)
xt = a * (th - np.sin(th))
yt = a * (1 - np.cos(th))
dx = np.gradient(xt, th); dy = np.gradient(yt, th)
d2x = np.gradient(dx, th); d2y = np.gradient(dy, th)
kappa = np.abs(dx * d2y - dy * d2x) / (dx ** 2 + dy ** 2) ** 1.5
for tv in [math.pi / 2, math.pi / 2 + 1.0, math.pi - 0.02]:
    i = int(np.argmin(np.abs(th - tv)))
    print(f"theta = {th[i]:.6f}: 数值曲率 = {kappa[i]:.8f}, "
          f"理论值 = {1 / (4 * a * math.sin(th[i] / 2)):.8f}")

输出：

          圆（半径 1）: L = 6.283185, A = 3.141593, L^2/(4 pi A) = 1.000000, 曲率 [1.0000, 1.0000]
    椭圆 a=1.5, b=1: L = 7.932720, A = 4.712389, L^2/(4 pi A) = 1.062657, 曲率 [0.4444, 1.5000]
      椭圆 a=2, b=1: L = 9.688448, A = 6.283185, L^2/(4 pi A) = 1.188827, 曲率 [0.2500, 2.0000]
超椭圆 |x|^4+|y|^4=1: L = 7.017698, A = 3.708149, L^2/(4 pi A) = 1.056872, 曲率 [0.0000, 2.5227]
theta = 1.570796: 数值曲率 = 0.70710657, 理论值 = 0.70710678
theta = 2.570796: 数值曲率 = 0.52107780, 理论值 = 0.52107783
theta = 3.121590: 数值曲率 = 0.50002407, 理论值 = 0.50002501

第一段数据同时验证了几何量的计算与等周不等式。圆的三个量都落在理论值上（$L=2\pi =6.283185$、$A=\pi =3.141593$、比值 $1.000000$），曲率恒为 $1$（周期中心差分的误差在 ${10}^{-8}$ 以下）。椭圆的曲率范围与解析公式 $[\frac{b}{a^2},\frac{a}{b^2}]$ 一致（$a=1.5$ 时 $[0.4444,1.5000]$，$a=2$ 时 $[0.2500,2.0000]$），说明曲率的数值计算可靠。等周比值全部不小于 $1$：圆 $1.000000$、超椭圆 $1.056872$、椭圆（$1.5,1$）$1.062657$、椭圆（$2,1$）$1.188827$，正方形作为解析对照是 $1.273240$。数值的顺序与「形状偏离圆的程度」一致。

超椭圆一行的曲率范围是 $[0.0000,2.5227]$，最小值为零。这一点与它的形状相符：$|x{|}^4+|y{|}^4=1$ 在坐标轴附近写成 $x\approx 1-\frac{y^4}{4}$，四阶的平坦使得曲率趋于零，而「角点」附近的曲率较大（$2.5227$）。这条数据说明曲率的数值范围是刻画形状尖锐程度的有效指标：椭圆的最大曲率来自短轴端点（光滑但弯曲），超椭圆的最大曲率来自近似角点（曲率集中）。

第二段数据核对摆线的曲率公式 $\kappa =\frac{1}{4a\sin (\theta /2)}$：三个取样点上数值与理论值的相对偏差都小于 $4\times {10}^{-5}$（$0.70710657$ 对 $0.70710678$；$0.52107780$ 对 $0.52107783$；$0.50002407$ 对 $0.50002501$）。$\theta \to 0$（起点）处曲率趋于无穷，这与「摆线在起点处切线竖直、曲线急剧转向」的几何图像一致，也是最速降线的起始段能有很大加速度的原因。数值实现的一个注意点：曲率涉及二阶导数，在端点附近的差分估计误差较大（本节实验中若直接取端点值，误差可达 $50\mathrm{\%}$），实际计算中应当避开端点，或者用周期/解析导数。

5.11.7 本节小结

几何与变分的四层内容

本节把曲线曲面几何与变分法组织为四层。曲线层给出弧长 $L=\int \parallel {\gamma }^{\prime }\parallel dt$ 与曲率 $\kappa =\frac{|x^{\prime }{y}^{″}-y^{\prime }{x}^{″}|}{(x^{\prime 2}+y^{\prime 2}{)}^{3/2}}$，曲率度量偏离直线的程度（直线为零、圆为 $\frac{1}{R}$、椭圆的取值范围是 $[\frac{b}{a^2},\frac{a}{b^2}]$，实验四的数据：$a=1.5$ 时 $[0.4444,1.5000]$，$a=2$ 时 $[0.2500,2.0000]$）；Frenet 公式 $d\tau /ds=\kappa n$ 把曲率解释为标架的转动速率。曲面层给出两个基本形式：第一基本形式 $I=Ed{u}^2+2Fdudv+Gd{v}^2$ 是切平面上的度量（决定长度与面积 $dS=\sqrt{EG-F^2}dudv$），第二基本形式度量曲面偏离切平面的程度，两者的组合给出主曲率、高斯曲率 $K={\kappa }_1{\kappa }_2$ 与平均曲率 $H$；测地线是曲面上「最直」与「最短」合一的曲线（球的测地线是大圆，平面的是直线）。变分层给出泛函的一阶变分为零等价于 Euler-Lagrange 方程 $\frac{\partial F}{\partial y}=\frac{d}{dx}\frac{\partial F}{\partial y^{\prime }}$（推导依赖分部积分与边界条件给出的 $\eta$ 在端点为零），$F$ 不显含 $x$ 时有守恒量 $F-y^{\prime }\frac{\partial F}{\partial y^{\prime }}=\text{常数}$。三个经典例子的数据：最短路径解是直线；最速降线解是摆线，$T=\sqrt{a/g}{\theta }_1=0.7096134476$ 秒，直线要慢 $18.5447\mathrm{\%}$，起始斜率减到 $0.7$ 倍与 $0.4$ 倍时分别慢 $29.5423\mathrm{\%}$ 与 $49.7834\mathrm{\%}$；等周问题的解是圆，实验用 $L^2\ge 4\pi A$ 核对（圆 $1.000000$、椭圆 $1.062657$ 与 $1.188827$、超椭圆 $1.056872$、正方形 $1.273240$）。计算层给出直接方法与有限元：弱形式把二阶方程降为一阶导数的积分等式，$H^1$ 与 $H_0^1$ 提供合适的函数空间，Ritz-Galerkin 方法把极小化化为线性方程组 $Ku=b$，$L^2$ 误差按 $O(h^2)$ 下降（实验一与实验三的数据，比值稳定在 $4.00$）；一维分片线性有限元的节点值有超收敛性（实测节点误差按 $h^6$ 下降，比值 $65.44$、$64.35$、$64.02$，残余误差来自载荷积分的三点 Gauss 求积），而这一现象的成立条件是刚度矩阵与载荷向量都用足够精确的积分装配。直接离散化遇上解的奇性时收敛会变慢：摆线在起点处的 $y\sim C{x}^{2/3}$ 使分段线性近似的偏差只按 $O(h^{2/3})$ 下降（实测比值 $1.61$，接近 $2^{2/3}=1.587$）。

变分法处理的是「在所有函数中找极值」，本节给出了它的骨架与两类求解路线。下一节处理另一类工具：当问题中的参数变得很大或很小时，函数的形状出现极端（尖峰、陡峭的过渡层、快速振荡），此时需要渐近分析：不求精确解，而是求主项与修正项的显式表达式。Stirling 公式、Laplace 方法与边界层理论都属于这一类，它们给出的近似往往比数值解的公式更直观。

练习题

第 1 题 概念推导

(a) 用 Euler-Lagrange 方程求最短路径问题的解，并说明为什么结果是直线。(b) 写出摆线 $x=a(\theta -\sin \theta )$、$y=a(1-\cos \theta )$ 的曲率 $\kappa =\frac{1}{4a\sin (\theta /2)}$ 的推导要点。(c) 说明弱形式与经典方程的关系，以及为什么弱解的存在性比经典解更容易证明。

参考答案

(a) 最短路径问题的泛函是 $J[y]={\int }_a^b\sqrt{1+y^{\prime 2}}dx$，即 $F=\sqrt{1+y^{\prime 2}}$。计算偏导数：

$$\frac{\partial F}{\partial y}=0,\frac{\partial F}{\partial y^{\prime }}=\frac{y^{\prime }}{\sqrt{1+y^{\prime 2}}}$$

Euler-Lagrange 方程给出 $\frac{d}{dx}\frac{y^{\prime }}{\sqrt{1+y^{\prime 2}}}=0$，即 $\frac{y^{\prime }}{\sqrt{1+y^{\prime 2}}}=C$（常数）。解出 $y^{\prime }=\frac{C}{\sqrt{1-C^2}}$，右端是常数，因此 $y=mx+n$ 是直线。$|C|<1$ 时解存在，$C=\pm 1$ 对应竖直直线（需要把 $y$ 作为自变量重新参数化）。

这条推导的结论与平面几何的直觉一致，但路径不同：几何直觉给出最短性，变分推导给出必要条件（驻点），要确认是极小还需要二阶变分或凸性论证（本例中 $F$ 关于 $y^{\prime }$ 是凸的，因此驻点即极小点）。这条练习说明变分法的一般流程：写出 $F$，算两个偏导数，得微分方程，解方程，再用二阶条件或凸性确认极值的类型。

(b) 用参数形式的曲率公式 $\kappa =\frac{|x^{\prime }{y}^{″}-y^{\prime }{x}^{″}|}{(x^{\prime 2}+y^{\prime 2}{)}^{3/2}}$（对 $\theta$ 求导）：

$$x^{\prime }=a(1-\cos \theta ),y^{\prime }=a\sin \theta ,x^{″}=a\sin \theta ,y^{″}=a\cos \theta$$

代入分子：$x^{\prime }{y}^{″}-y^{\prime }{x}^{″}=a^2[(1-\cos \theta )\cos \theta -{\sin }^2\theta ]=a^2(\cos \theta -{\cos }^2\theta -{\sin }^2\theta )=a^2(\cos \theta -1)$。分母：$x^{\prime 2}+y^{\prime 2}=a^2(1-\cos \theta {)}^2+a^2{\sin }^2\theta =2a^2(1-\cos \theta )$，其 $\frac{3}{2}$ 次方为 ${[2a^2(1-\cos \theta )]}^{3/2}$。因此

$$\kappa =\frac{a^2(1-\cos \theta )}{{[2a^2(1-\cos \theta )]}^{3/2}}=\frac{1}{a\sqrt{2}\sqrt{1-\cos \theta }}=\frac{1}{4a\sin (\theta /2)}$$

最后一步用了半角公式 $\sqrt{1-\cos \theta }=\sqrt{2}\sin (\theta /2)$。实验四的三个取样点（$\theta =\pi /2$、$\pi /2+1$、$\pi -0.02$）给出的数值与这个公式的相对偏差都小于 $4\times {10}^{-5}$。

(c) 经典方程（Euler-Lagrange 方程或 $-(p{y}^{\prime }{)}^{\prime }=f$）要求解有二阶导数，因为方程中出现 $y^{″}$。弱形式把方程乘以试探函数 $\eta \in H_0^1$ 并分部积分：

$${\int }_a^b(p{y}^{\prime }{\eta }^{\prime }-f\eta )dx=0,\mathrm{\forall }\eta \in H_0^1$$

等式两边只含 $y$ 与 $\eta$ 的一阶（弱）导数，因此对分段光滑、只有拐角（$y^{\prime }$ 不连续）的函数也有意义。弱解比经典解更容易存在的原因有二：一是函数空间 $H_0^1$ 完备（5.10 节的 Riesz-Fischer 定理），极小化序列可以在其中收敛；二是能量泛函 $J$ 关于 $H_0^1$ 范数是强制的（$J\to \mathrm{\infty }$ 当 $\parallel u{\parallel }_{H^1}\to \mathrm{\infty }$）与凸的，凸泛函的驻点即极小点（第三章 3.2 节的凸性理论）。两步合起来给出弱解的存在性，再通过正则性理论证明弱解在数据光滑时就是经典解。实验三的数据（$L^2$ 误差 $O(h^2)$）说明有限元近似在弱解框架下是合理的；若改用经典解的框架，低光滑度的解根本不能进入讨论。

第 2 题 计算推理

(a) 用 Beltrami 恒等式求极小旋转曲面的母线方程（泛函 $\int y\sqrt{1+y^{\prime 2}}dx$）。(b) 计算曲率 $\kappa =1$ 的平面曲线族（提示：$\frac{d\theta }{ds}=\kappa$ 给出圆心在法线方向上的圆）。(c) 用等周不等式与定积分估计：周长固定为 $L$ 的闭曲线，面积的最大值是多少。

参考答案

(a) 泛函是 $J[y]=\int y\sqrt{1+y^{\prime 2}}dx$（旋转曲面的面积），$F=y\sqrt{1+y^{\prime 2}}$ 不显含 $x$，用 Beltrami 恒等式

$$F-y^{\prime }\frac{\partial F}{\partial y^{\prime }}=C_1$$

计算 $\frac{\partial F}{\partial y^{\prime }}=\frac{y{y}^{\prime }}{\sqrt{1+y^{\prime 2}}}$，因此

$$y\sqrt{1+y^{\prime 2}}-\frac{y{y}^{\prime 2}}{\sqrt{1+y^{\prime 2}}}=\frac{y}{\sqrt{1+y^{\prime 2}}}=C_1$$

即 $\sqrt{1+y^{\prime 2}}=\frac{y}{C_1}$，$y^{\prime 2}=\frac{y^2}{C_1^2}-1$，解为

$$y=C_1\cosh \frac{x-C_2}{C_1}$$

悬链线。因此旋转极小曲面是悬链面，这条曲线在建筑与力学中也出现（均匀悬链的形状、拱的设计）。这个例子说明 Beltrami 恒等式的作用：一阶守恒量把二阶方程降为一阶，解可以直接分离变量。

(b) 曲率恒为 $1$ 的曲线满足 $\frac{d\theta }{ds}=\pm 1$（$\theta$ 是切向量的倾角），因此 $\theta =\pm s+{\theta }_0$。切向量是 $(\cos \theta ,\sin \theta )$，积分得

$$x(s)=x_0+{\int }_0^s\cos (\pm t+{\theta }_0)dt,y(s)=y_0+{\int }_0^s\sin (\pm t+{\theta }_0)dt$$

积分给出圆心在 $(x_0\mp \sin {\theta }_0,y_0\pm \cos {\theta }_0)$、半径 $1$ 的圆（参数化后可以验证）。结论：曲率恒为正常数的平面曲线是圆，这是「等周问题的解是圆」的几何基础（曲率恒定是最优性的必要条件）。

(c) 由等周不等式 $L^2\ge 4\pi A$ 得

$$A\le \frac{L^2}{4\pi }$$

等号在圆处成立，此时 $A=\frac{L^2}{4\pi }$（验证：半径为 $R$ 的圆有 $L=2\pi R$、$A=\pi R^2=\frac{L^2}{4\pi }$）。数值上的对照：单位圆（$L=2\pi$）的面积上界是 $\frac{4{\pi }^2}{4\pi }=\pi$，等号成立；正方形的周长 $L=4$ 时上界是 $\frac{16}{4\pi }=1.273240$，而面积只有 $1$，比值 $1.273240$ 与实验四的解析对照一致。这条结论给出「面积效率」的定量上界，也是形状优化中「给定周长最大化面积」这类问题的答案。

第 3 题 代码验证

用本节代码做三组实验并解释：(a) 在实验一里把 $f$ 换成 $f(x)=1$（常值），检查离散解与精确解 $y=\frac{x(1-x)}{2}$ 的误差，并说明为什么这次误差会落到机器精度；(b) 在实验二里把终点改成 ${\theta }_1=\frac{\pi }{2}$（即 $a=0.5$、$x_1=0.5(\frac{\pi }{2}-1)$、$y_1=0.5$），重新比较直线与摆线；(c) 在实验三里把载荷的求积公式从三点 Gauss 换成中点公式，观察节点误差的收敛率如何变化。

参考答案

(a) $f=1$ 时精确解是 $y=\frac{x(1-x)}{2}$（二次多项式），它的离散化是精确的：离散泛函用梯形公式逼近 $\int (y^{\prime 2}-2y)dx$，其中 $y^{\prime }$ 是分段线性函数的导数（在每段上是常数），$y^{\prime 2}$ 的积分对二次多项式精确；线性项 $-2\int ydx$ 用梯形公式对线性函数精确；而离散 EL 方程 $-(y_{j-1}-2y_j+y_{j+1})=h^2{f}_j$ 对二次函数也是精确的（二阶中心差分对二次多项式精确）。三条精确性合起来给出：离散解在节点上等于精确解，误差落到机器精度（${10}^{-16}$ 量级），与 $N$ 无关。这条练习说明「离散化的误差」与「格式的代数精度」直接相关：格式对多项式精确的阶数越高，光滑解上的误差越小；极限情形（解为低次多项式）时误差可以精确为零。对照实验一（解为 $\sin \pi x$）的 $O(h^2)$，两者差别在于解是否落在格式精确覆盖的多项式空间内。

(b) $a=0.5$、${\theta }_1=\frac{\pi }{2}$ 时终点是 $x_1=0.5(\frac{\pi }{2}-1)=0.285398$、$y_1=0.5(1-0)=0.5$（注意 $1-\cos \frac{\pi }{2}=1$）。摆线的下滑时间 $T=\sqrt{\frac{a}{g}}{\theta }_1=0.225877\times 1.570796=0.354807$ 秒。直线的下滑时间由公式 $2\sqrt{1+1/m^2}\sqrt{y_1}/\sqrt{2g}$（$m=\frac{y_1}{x_1}=1.751938$）给出：$2\times 1.151439\times 0.707107/4.427189=0.367813$ 秒，比摆线多 $3.67\mathrm{\%}$。控制变量从 $18.54\mathrm{\%}$ 降到 $3.67\mathrm{\%}$ 的原因：终点更靠近起点时所有曲线的形状差别变小，最优曲线的优势随之减小。这条对照说明最优性的相对收益依赖问题的几何尺度，比较不同问题的最优解时需要用无量纲的比值（例如终点的高宽比 $y_1/x_1$）。

(c) 把三点 Gauss 换成中点公式后，载荷向量的误差从 $O(h^6)$ 降到 $O(h^2)$，因此节点误差的收敛率从 $O(h^6)$ 降到 $O(h^2)$（超收敛被破坏），而 $L^2$ 误差仍按 $O(h^2)$ 下降（因为它由逼近误差主导，与载荷求积的精度无关，只要求积误差不低于逼近误差即可）。在数值上表现为：节点误差列的比值从 $65$ 变成 $4$，两列误差的比值趋于相同。这条练习给出有限元实现的重要规则：各项误差的阶数需要匹配，载荷求积的精度至少要与逼近误差同阶；一味提高网格密度会在求积误差处碰到下限（与 5.5 节数值求积的「精度上限由模型误差决定」是同一条原则）。实践中的做法是对光滑载荷使用高阶 Gauss（两点或三点足以达到 $O(h^4)$、$O(h^6)$），对不光滑载荷（有间断、奇性）则加密网格或做自适应求积。

第 4 题 综合应用

某工程优化问题需要求解一条曲线的形状，使某个泛函取极值（泛函的被积函数只依赖 $y$ 与 $y^{\prime }$，边界条件为固定端点；被积函数光滑但对 $y^{\prime }$ 非线性；解在起点附近有一个 $(x-x_0{)}^{2/3}$ 型的行为）。请给出两条求解路线的比较与推荐：解析求解 Euler-Lagrange 方程、直接离散化（有限差分离散 + 数值优化），并说明各自的实现步骤、精度预期与需要注意的问题。

参考答案

路线一：解析求解 Euler-Lagrange 方程。步骤是写出被积函数 $F(x,y,y^{\prime })$，计算 $\frac{\partial F}{\partial y}$ 与 $\frac{\partial F}{\partial y^{\prime }}$，得到二阶常微分方程；若 $F$ 不显含 $x$，先用 Beltrami 恒等式 $F-y^{\prime }\frac{\partial F}{\partial y^{\prime }}=C$ 把方程降为一阶；再分离变量或参数化求解，最后用边界条件定常数。精度是精确的（若能得到闭式解），成本主要是代数推导与积分。适用条件是方程可积（有守恒量、可分离、或属于少数可积类型）。本例中被积函数光滑但对 $y^{\prime }$ 非线性，EL 方程一般是非线性二阶方程，能否积出闭式解取决于 $F$ 的具体形式；起点附近的 $(x-x_0{)}^{2/3}$ 行为提示解有奇性（导数无界），这类奇性往往对应一个守恒量给出的隐式关系（如最速降线的 $y(1+y^{\prime 2})=2a$），因此值得先试 Beltrami 恒等式。

路线二：直接离散化 + 数值优化。步骤是把区间划分为 $N$ 段，用节点值（或分段多项式系数）参数化曲线，把泛函用数值求积近似（单元内用 Gauss 求积，对分段线性曲线可以用解析的单元积分），得到有限维目标函数；施加边界条件后用 L-BFGS-B 或信赖域方法极小化；最后用网格加密检验收敛。精度预期：解光滑时偏差按 $O(h^2)$，但本例中 $y\sim C(x-x_0{)}^{2/3}$ 的奇性会把它降到 $O(h^{2/3})$（本节实验二的实测比值 $1.61\approx 2^{2/3}$ 给出确证）。需要注意的问题有三条：一是奇性导致的慢收敛，解决办法是在起点附近用非均匀网格（几何加密）或改用能捕捉 $x^{2/3}$ 的单元（例如带权重的单元、或直接把 $x$ 换成 $t^3$ 做参数化）；二是优化的收敛与初值，非线性泛函可能有多个局部极小，需要多初值试验；三是求积误差与离散误差的匹配（被积函数在奇点附近的求积需要特殊处理，否则会出现像本节实验中「用 $x$ 作积分变量时的 $x^{-1/2}$ 发散」那样的系统偏差）。

推荐：两条路线配合使用。先用解析路线（重点是 Beltrami 恒等式）尝试求出闭式解或至少弄清解在奇点附近的行为（例如确定 $y\sim C(x-x_0{)}^{2/3}$ 中的 $C$ 与 $x_0$），这一步给出的信息能直接指导数值方案：知道奇性的位置与阶数后，可以在数值路线中使用变量替换 $x-x_0=t^3$（把 $x^{2/3}$ 变成 $t^2$，恢复光滑），从而把收敛率从 $O(h^{2/3})$ 提到 $O(h^2)$。若解析路线无法积出闭式解，仍然保留数值路线作为主方案，并用「网格加密时误差按什么速率下降」来检验解的奇性推断是否成立。报告中给出：泛函的具体形式、EL 方程（或弱形式）、所用的网格与求积、收敛速率与奇性处理方式，这些信息合起来才能让结果可复核。

常见错误

错误 1 · 把 Euler-Lagrange 方程的驻点直接当作极小点

原因：一阶变分为零只给出驻点（极值点的必要条件），驻点可能是极小、极大或鞍点。变分问题中这种模糊性比有限维更常见：泛函可能不凸（例如泛函 $\int (y^{\prime 2}-y^2)dx$ 的二次部分是符号不定的，对应不定算子），此时驻点可以是极大或鞍点。另一类问题是泛函可能无下界（例如 ${\int }_0^1{y}^{\prime 2}dx$ 在边界条件只给两端点的情形下没有极小点：取越来越快的振荡可使泛函值任意小）。本节实验一与实验二的问题都是凸的（$F$ 关于 $y^{\prime }$ 是二次的），因此驻点就是极小点，但这不是普遍性质。

解决：求出 EL 方程的驻点后检查二阶变分（${\delta }^{2}J\ge 0$ 对一切扰动 $\eta$ 成立）、或检查泛函的凸性与强制性（coercivity）。凸泛函的驻点即极小点，这是最方便的判据（第三章 3.2 节的凸性理论）；非凸情形需要多初值试验或全局优化方法。报告中说明「求得的是驻点」，只有在补充论证后才能称为极小点。

错误 2 · 在被积函数有奇性时直接用常规求积公式

原因：变分问题的被积函数常在端点或解的特殊点处有奇性（最速降线在起点的 $\frac{1}{\sqrt{2gy}}$、极小旋转曲面在小半径处的形状、含 $(x-x_0{)}^{2/3}$ 解的问题）。常规求积（梯形、Simpson、固定阶 Gauss）在这些点的附近误差很大，会出现系统性偏差。本节实验二给出了量化证据：用 $x$ 作积分变量时被积函数在起点按 $x^{-1/2}$ 发散，梯形法必须截断掉起点的邻域，直线的下滑时间因此被算成 $0.8399054730$ 秒，而解析值是 $0.8412091759$ 秒，相对误差 $1.55\times {10}^{-3}$；改用 $y$ 作积分变量并令 $y=u^2$ 后得到 $0.8412091759$ 秒，与解析值一致到 $6.6\times {10}^{-16}$。两者的差别全部来自「奇异邻域被截断」这一个来源，与网格密度无关，加密网格不会消除它。

解决：先做奇性分析（找出发散点的位置与阶数），再选择处理方式：变量替换（把积分变量移到使被积函数光滑的那一侧，如实验二用 $y$ 代替 $x$；或令 $x-x_0=t^m$ 使奇性变成幂函数的可积形式）、加权求积（Gauss-Jacobi 型公式，权重含奇性因子）、或解析处理奇异部分。处理完成后用两个独立方法交叉核对（解析值与数值值、两种变量替换的结果），确认没有系统偏差再使用结果。

错误 3 · 用同一分辨率的结果判断离散解是否收敛

原因：直接离散化的误差依赖网格尺寸，单次计算只能给出「当前网格下的近似」，无法区分「已经收敛」与「恰好在这一分辨率下看起来合理」。本节实验二的数据给出警示：$N=20$ 时离散最小值是 $0.7126498323$，与连续极小值 $0.7096134476$ 相差 $0.4\mathrm{\%}$，看起来「已经很接近」；但把网格加密到 $N=320$ 时偏差降到 $0.03\mathrm{\%}$，说明粗网格的结果只是恰好接近，收敛过程仍在继续。更严重的情形是解有奇性时（如本例的 $O(h^{2/3})$），粗网格的偏差可能掩盖真实的收敛率，据此判断「算法精度不够」会指向错误的方向。

解决：至少用两到三级网格计算，报告误差随网格的变化与估计的收敛阶（用 $\log (e_1/e_2)/\log (h_1/h_2)$）；若收敛阶明显低于理论值，检查解的奇性与求积公式的精度（错误 2）；若收敛阶与理论值相符，可以用 Richardson 外推提高精度并给出误差估计。报告中给出多级网格的结果与收敛阶，避免只报单一分辨率的数值。

错误 4 · 混淆弱解与经典解，或忽略边界条件的作用

原因：弱形式只要求一阶（弱）导数可积，因此弱解可能有拐角、kink，甚至导数无界（本节实验二的 $x^{2/3}$ 行为）；而经典解要求二阶导数连续。把弱解当作经典解使用会出错：例如对分段线性插值的有限元解直接求二阶导数，得到的是集中在节点上的分布（导数在单元内为零、跨单元跳变），不能当作 $u^{″}$ 使用。另一个常见疏漏是边界条件：Dirichlet 条件必须写进试探空间（$\eta \in H_0^1$），Neumann 条件则是自然边界条件（从分部积分的边界项自动出现，不需要强加）。把两类边界条件混用，或忘记在最小化时固定端点，会得到错误的解。

解决：明确所用的是弱解还是经典解，报告中说明解的正则性（在哪些点可能不光滑）。求导数时事先判断解的光滑程度：需要二阶导数时先做后处理（例如用导数恢复技术把分段常数导数光滑化），或者改用混合元方法把导数本身作为未知量求解。边界条件的处理逐条核对：Dirichlet 条件体现在试探空间与试验空间的取法上，Neumann 与 Robin 条件体现在弱形式的边界项上（当边界项不为零时把它移到右端作为载荷的一部分）。