椭圆积分——从不会到会 still 不会
前言(可跳过)
本文重点在于锻炼读者的"注意力"——不是让你背下更多公式,而是让你看见公式背后的结构。
我自己检验是否真正理解一个公式的方法是:你能否用自然语言精确描述它给你的感觉?
比如 ∫1−x21dx=arcsinx,你看到的是"一个公式",还是"单位圆上角度与弧长的互换"?
椭圆积分之所以难,不是因为它公式复杂,而是因为它迫使你在"角度"和"弧长"之间建立一一对应——这种对应在圆里是平凡的(弧长=半径×角度),在椭圆里却需要全新的函数来描述。
本文会给出所有关键公式的几何意义和证明思路。作者实力有限,可能无法道尽,希望知友评论区补充指正。
前置知识
- 三角函数、反三角函数
- 一点定积分
- 知道 1−x2 积出来是 arcsin
- 会求导
- 最好见过椭圆方程
其中第二条是很多人会忽略的,请读者务必最起码会使用定积分的几何意义。
约定俗成:下文用 E 表示第二类椭圆积分、F 表示第一类椭圆积分、Π 表示第三类。
引子:椭圆弧长——为什么需要新函数
你画一个椭圆
a2x2+b2y2=1
想用高中那一套求个周长,然后发现弧长公式长这样:
s=∫02πa2sin2θ+b2cos2θdθ
化简一下,令离心率 e=1−a2b2,a>b:
s=4a∫02π1−e2sin2θdθ
这就是第二类完全椭圆积分 E(e)。
(图源:AI生成,模拟手绘风格)
最重要的直觉:
- 圆的弧长公式 s=rθ 之所以简单,是因为圆的曲率处处相等——你转过的角度和走过的弧长成严格正比
- 椭圆的曲率在变化——靠近长轴端点曲率大,靠近短轴端点曲率小
- 所以"转过的角度"和"走过的弧长"不再成正比,E(e) 就是对这种非线性关系的精确编码
所以椭圆积分最早就是算椭圆周长算出来的,名字很诚实。
第一类椭圆积分 F(φ,k)
定义与几何意义
F(φ,k)=∫0φ1−k2sin2θdθ
其中 k∈[0,1) 叫模数,φ 叫振幅。
几何意义:F(φ,k) 是椭圆上的"归一化弧长参数"。
具体地,考虑椭圆 x=asinθ,y=bcosθ。从 θ=0 到 θ=φ 的弧长为
s(φ)=a∫0φ1−e2sin2θdθ=aE(φ,e)
而 F(φ,k) 与 E(φ,k) 通过反函数关系耦合:如果 u=F(φ,k),那么 φ 就是 u 的振幅,记作 φ=am(u,k)。
为什么 φ=2π 叫"完全":
因为 F(φ,k) 关于 φ 是周期性的,周期为 4K(k)(K(k) 是完全椭圆积分)。当 φ 从 0 走到 2π 时,点在椭圆上从短轴端点走到长轴端点,恰好遍历了椭圆的四分之一周长。这就是为什么 φ=2π 时叫"完全"——它对应着四分之一周期,而不是整个周期。
当 φ=2π 时,记作 K(k):
K(k)=∫02π1−k2sin2θdθ
特殊值
- k=0:K(0)=2π,退化回普通积分。
- k→1−:K(k)→+∞,因为 1−sin2θ 在 θ=2π 附近变成 0,分母爆炸。
它积出来是啥?
很不幸,第一类椭圆积分没有初等函数形式。就像 ∫e−x2dx 积不出来要叫 erf 一样,椭圆积分也要被赐名为 F、K。
但你可以查表,或者用 AGM 算法(后面讲)。
直觉:F(φ,k) 衡量的是**"角度被压缩的程度"**。k 越大,sin2θ 越"有效",分母越小,同样的 dθ 对应的 F 增量越大。
第二类椭圆积分 E(φ,k)
定义:
E(φ,k)=∫0φ1−k2sin2θdθ
完全形式:
E(k)=∫02π1−k2sin2θdθ
特殊值
- k=0:E(0)=2π。
- k=1:E(1)=1(因为 ∫02πcosθdθ=1)。
椭圆周长公式
椭圆周长 C=4aE(e),其中 e=1−a2b2。所以别再说椭圆周长没公式了,它有,只是不是初等函数。
直觉:E(φ,k) 和 F(φ,k) 是对偶的。F 的被积函数在分母,E 的在分子。k 增大时,F 发散(分母趋零),E 收敛到 1(分子趋零)。
第三类椭圆积分 Π(n;φ,k)
定义:
Π(n;φ,k)=∫0φ(1−nsin2θ)1−k2sin2θdθ
多了一个参数 n,叫特征值。完全形式记为 Π(n;k)。
这个用得比前两类少,但物理里经常冒出来(比如电磁学、摆的大幅度周期)。
三个参数,非常不讲道理,建议先记住形式。
它们之间的关系——Legendre 关系
E(k)K(k′)+E(k′)K(k)−K(k)K(k′)=2π
其中 k′=1−k2 叫互补模数。
这个式子为什么成立?
证明思路:Legendre 关系可以通过对 K(k) 和 E(k) 分别求导,利用它们的微分方程联立消元得到。
具体地,我们知道
dkdK=k(1−k2)E(k)−(1−k2)K(k),dkdE=kE(k)−K(k)
考虑函数 f(k)=E(k)K(k′)+E(k′)K(k)−K(k)K(k′)。
对 f(k) 求导,利用 k′=1−k2 和链式法则(dkdk′=−k′k),代入上述导数公式,经过一系列化简(约一页纸),可以发现 dkdf=0。
所以 f(k) 是常数。代入 k=0(此时 k′=1,K(0)=E(0)=2π,K(1)→∞ 但 E(1)=1 使乘积有限),计算得 f(0)=2π。
因此 f(k)=2π 对所有 k 成立。
直觉:Legendre 关系是椭圆积分周期性的体现。K(k) 和 E(k) 分别衡量椭圆积分的"第一类"和"第二类"行为,而 k′ 代表"正交方向"的模数。这个关系保证了无论 k 怎么变,K 和 E 的乘积组合始终保持不变——就像能量守恒一样。
怎么记住它:
- 三项轮换:EK′+E′K−KK′
- 常数项:2π(来自 k=0 时的平凡值)
- 对称性:k↔k′ 时式子不变
求导和恒等式
dkdK(k)=k(1−k2)E(k)−(1−k2)K(k)
dkdE(k)=kE(k)−K(k)
证明基本就是积分号下求导然后化简,比较机械,但草稿纸需要大一点。
直觉:K 的导数包含 E,E 的导数包含 K。它们互相定义,就像 sin 和 cos 的导数互相包含对方一样。
AGM 算法——手算 K(k)
K(k) 可以用算术-几何平均(AGM)快速算:
设 a0=1,b0=1−k2=k′,然后迭代:
an+1=2an+bn
bn+1=anbn
cn+1=2an−bn
当 n 足够大时,an≈bn≈M(a0,b0),AGM 收敛非常快,每轮有效位数翻倍。
第一类完全椭圆积分:
K(k)=2M(1,k′)π
第二类完全椭圆积分:
E(k)=2M(a0,b0)π(1−n=0∑∞2n−1cn2)
为什么 AGM 能算出椭圆积分?
证明思路:AGM 与椭圆积分的联系来自Landen 变换。Landen 变换是一种椭圆积分的变量替换,它将 F(φ,k) 转化为另一个椭圆积分,其模数 k1=1+k′1−k′ 比原来的 k 更小。反复应用 Landen 变换,模数趋近于 0,积分趋近于 2π。
AGM 正是 Landen 变换的代数形式:an+1=2an+bn 和 bn+1=anbn 分别对应 Landen 变换中的算术和几何操作。当 an 和 bn 收敛到同一个值 M 时,积分值被"压缩"成了 2Mπ。
(图源:AI生成,模拟手绘风格)
最重要的直觉:
- AGM 是算术平均和几何平均的迭代,两者收敛到同一个值
- 每轮迭代,cn(差的一半)平方级减小,所以收敛极快
- K(k) 就是 2π 除以这个收敛值——椭圆积分被"压缩"成了一个极限
代码实现:
#include<bits/stdc++.h>
using namespace std;
typedef long double ld;
const ld PI = acosl(-1.0);
// 返回 K(k), E(k)
pair<ld, ld> elliptic_K_E(ld k) {
ld a = 1.0, b = sqrtl(1 - k * k), c = k;
ld sum = 0.0;
ld pw2 = 1.0; // 2^n
for (int n = 0; n < 10; ++n) { // 10 次迭代精度足够
ld an = (a + b) / 2.0;
ld bn = sqrtl(a * b);
ld cn = (a - b) / 2.0;
sum += pw2 * cn * cn;
pw2 *= 2.0;
a = an; b = bn; c = cn;
}
ld K = PI / (2 * a);
ld E = K * (1.0 - sum / 2.0); // 注意系数写法
return {K, E};
}
int main() {
ld k = 0.5;
auto [K, E] = elliptic_K_E(k);
cout << "K(" << k << ") = " << K << "\n";
cout << "E(" << k << ") = " << E << "\n";
return 0;
}
坑点:
- k 太接近 1 时 K 会很大,double 可能不够,建议用 long double。
- 初始 b0 是 1−k2,不是 k。
- 上面的 E 公式系数容易记混,建议每次用 E(0)=2π 和 E(1)=1 验证。
椭圆函数 Jacobi sn, cn, dn
椭圆积分是 Jacobi 椭圆函数的反函数。就像 arcsin 和 sin 互为反函数一样。
定义:如果
u=F(φ,k)=∫0φ1−k2sin2θdθ
那么 Jacobi 椭圆函数:
sn(u,k)=sinφ
cn(u,k)=cosφ
dn(u,k)=1−k2sin2φ
(图源:AI生成,模拟手绘风格)
基本恒等式
sn2(u,k)+cn2(u,k)=1
k2sn2(u,k)+dn2(u,k)=1
导数
dudsn(u,k)=cn(u,k)dn(u,k)
dudcn(u,k)=−sn(u,k)dn(u,k)
duddn(u,k)=−k2sn(u,k)cn(u,k)
当 k=0 时 sn→sin,cn→cos,dn→1。
最重要的直觉:
- sn 和 cn 是 sin 和 cos 的推广——单位圆被"压扁"成了椭圆
- dn 是椭圆独有的,衡量"压扁程度"
- 导数公式里 dn 的出现,是因为椭圆上"角度"和"弧长"的关系不再像圆那样简单
例题
例 1:单摆大角度周期
单摆小角度周期 T0=2πgL。大角度精确解:
T=4gLK(sin2θ0)
其中 θ0 是最大摆角。所以第一类完全椭圆积分就是单摆的修正因子。
直觉:摆角越大,摆锤在最高点附近"停留"的时间越长(因为势能曲线变平),周期变长。K(sin2θ0) 就是对这个"停留效应"的精确量化。
例 2:椭圆周长
椭圆 a2x2+b2y2=1,a>b,离心率 e=1−a2b2。
C=4aE(e)
如果 a=5,b=3,则 e=54,C=20E(0.8)。
例 3:数值验证
用 AGM 算 K(0) 和 E(0),应该得到 2π。如果算出来不对,说明你 b0 写错了。
大坑
注意一个很神奇的问题,比如上面 AGM 代码,你可能会发现我将
ld sum = 0.0;
ld pw2 = 1.0; // 2^n
注释了个"2^n",为什么呢?
举个例子,你把它写成了 2n−1,代码变成了下面这样。
sum += pw2 * cn * cn; // pw2 = 2^n
可能乍一看没问题,毕竟系数差个 2 而已。
BUT,当你把目光放到 E(k) 的公式那部分的时候,你就发现
E(k)=2M(a0,b0)π(1−n=0∑∞2n−1cn2)
当 n=0 时,2n−1=21,所以 c02 的系数是 21,不是 1!!!!!
如果你代码里写成了 2n,那么 c02 的系数就是 1,E(k) 就会多减 21c02,直接导致结果错误!!!
正确代码里 sum += pw2 * cn * cn 之后 pw2 *= 2.0,这样第 n 次迭代时 pw2 正好是 2n,但公式里是 2n−1,所以最后要 (1.0 - sum / 2.0) 来补偿。
总结
- 椭圆积分有三类:F,E,Π,分别是 ⋯1、⋯、(1−nsin2)⋯1 的积分。
- 完全椭圆积分就是上限到 2π,对应椭圆的四分之一周期。
- 椭圆函数是椭圆积分的反函数,Jacobi sn/cn/dn 是三角函数的推广。
- AGM 能高效算 K 和 E,其原理是 Landen 变换。
- Legendre 关系是椭圆积分周期性的体现,可以通过微分方程证明。
在做微积分的时候,看到巧合式子,一定要多想。想明白为什么,会给你的数学直觉带来很大的帮助。
(图源:AI生成,模拟手绘风格)
参考与推荐阅读
- OEIS 和各种手册里的椭圆积分表
- Abramowitz & Stegun: Handbook of Mathematical Functions
- 如果你想硬刚推导:Whittaker & Watson《A Course of Modern Analysis》