椭圆积分——从不会到会 still 不会

本章又名：积分积着积着发现反函数变成椭圆函数了，很气。

前置知识

三角函数、反三角函数、一点定积分、知道 $\sqrt{1-x^2}$ 积出来是 $\arcsin$、会求导、最好见过椭圆方程。

约定俗成：下文用 $E$ 表示第二类椭圆积分、$F$ 表示第一类椭圆积分、$\Pi$ 表示第三类，虽然这仨符号长得完全不像人名但国际数学界都这么叫，忍一下。

引子：椭圆弧长

你画一个椭圆

$$\frac{x^2}{a^2}+\frac{y^2}{b^2}=1$$

想用高中那一套求个周长，然后发现弧长公式长这样：

$$s=\int_0^{2\pi}\sqrt{a^2\sin^2\theta+b^2\cos^2\theta}\,d\theta$$

化简一下，令离心率 $e=\sqrt{1-\frac{b^2}{a^2}}$，$a>b$：

$$s=4a\int_0^{\frac{\pi}{2}}\sqrt{1-e^2\sin^2\theta}\,d\theta$$

这就是第二类完全椭圆积分 $E(e)$。

所以椭圆积分最早就是算椭圆周长算出来的，名字很诚实。

第一类椭圆积分 $F(\varphi,k)$

定义：

$$F(\varphi,k)=\int_0^{\varphi}\frac{d\theta}{\sqrt{1-k^2\sin^2\theta}}$$

其中 $k\in[0,1)$ 叫模数，$\varphi$ 叫振幅。

当 $\varphi=\frac{\pi}{2}$ 时，叫完全椭圆积分，记作 $K(k)$：

$$K(k)=\int_0^{\frac{\pi}{2}}\frac{d\theta}{\sqrt{1-k^2\sin^2\theta}}$$

特殊值

• $k=0$：$K(0)=\frac{\pi}{2}$，退化回普通积分。
• $k\to 1^-$：$K(k)\to +\infty$，因为 $1-\sin^2\theta$ 在 $\theta=\frac{\pi}{2}$ 附近变成 $0$，分母爆炸。

它积出来是啥？

很不幸，第一类椭圆积分没有初等函数形式。就像 $\int e^{-x^2}dx$ 积不出来要叫 erf 一样，椭圆积分也要被赐名为 $F$、$K$。

但你可以查表，或者用 AGM 算法（后面讲）。

第二类椭圆积分 $E(\varphi,k)$

定义：

$$E(\varphi,k)=\int_0^{\varphi}\sqrt{1-k^2\sin^2\theta}\,d\theta$$

完全形式：

$$E(k)=\int_0^{\frac{\pi}{2}}\sqrt{1-k^2\sin^2\theta}\,d\theta$$

特殊值

• $k=0$：$E(0)=\frac{\pi}{2}$。
• $k=1$：$E(1)=1$（因为 $\int_0^{\frac{\pi}{2}}\cos\theta d\theta = 1$）。

椭圆周长公式

椭圆周长 $C=4aE(e)$，其中 $e=\sqrt{1-\frac{b^2}{a^2}}$。所以别再说椭圆周长没公式了，它有，只是不是初等函数。

第三类椭圆积分 $\Pi(n;\varphi,k)$

定义：

$$\Pi(n;\varphi,k)=\int_0^{\varphi}\frac{d\theta}{(1-n\sin^2\theta)\sqrt{1-k^2\sin^2\theta}}$$

多了一个参数 $n$，叫特征值。完全形式记为 $\Pi(n;k)$。

这个用得比前两类少，但物理里经常冒出来（比如电磁学、摆的大幅度周期）。

三个参数，非常不讲道理，建议先记住形式。

它们之间的关系

Legendre 关系：

$$E(k)K(k')+E(k')K(k)-K(k)K(k')=\frac{\pi}{2}$$

其中 $k'=\sqrt{1-k^2}$ 叫互补模数。

这个式子看起来对称，实际上可以拿来当恒等式用，也可以拿来检查 AGM 算得对不对。

求导和恒等式

$$\frac{dK(k)}{dk}=\frac{E(k)-(1-k^2)K(k)}{k(1-k^2)}$$

$$\frac{dE(k)}{dk}=\frac{E(k)-K(k)}{k}$$

证明基本就是积分号下求导然后化简，比较机械，但草稿纸需要大一点。

AGM 算法——手算 $K(k)$

$K(k)$ 可以用算术-几何平均（AGM）快速算：

设 $a_0=1$，$b_0=\sqrt{1-k^2}=k'$，然后迭代：

$$a_{n+1}=\frac{a_n+b_n}{2}$$

$$b_{n+1}=\sqrt{a_nb_n}$$

$$c_{n+1}=\frac{a_n-b_n}{2}$$

当 $n$ 足够大时，$a_n\approx b_n\approx M(a_0,b_0)$，AGM 收敛非常快，每轮有效位数翻倍。

第一类完全椭圆积分：

$$K(k)=\frac{\pi}{2M(1,k')}$$

第二类完全椭圆积分：

$$E(k)=\frac{\pi}{2M(a_0,b_0)}\left(1-\sum_{n=0}^{\infty}2^{n-1}c_n^2\right)$$

每次看到这个求和都感叹，到底是谁先想到的……

代码实现：

#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;
}

坑点：

1. $k$ 太接近 $1$ 时 $K$ 会很大，double 可能不够，建议用 long double。
2. 初始 $b_0$ 是 $\sqrt{1-k^2}$，不是 $k$。
3. 上面的 $E$ 公式系数容易记混，建议每次用 $E(0)=\frac{\pi}{2}$ 和 $E(1)=1$ 验证。

椭圆函数 Jacobi sn, cn, dn

椭圆积分是 Jacobi 椭圆函数的反函数。就像 $\arcsin$ 和 $\sin$ 互为反函数一样。

定义：如果

$$u=F(\varphi,k)=\int_0^{\varphi}\frac{d\theta}{\sqrt{1-k^2\sin^2\theta}}$$

那么 Jacobi 椭圆函数：

$$\operatorname{sn}(u,k)=\sin\varphi$$ $$\operatorname{cn}(u,k)=\cos\varphi$$ $$\operatorname{dn}(u,k)=\sqrt{1-k^2\sin^2\varphi}$$

基本恒等式

$$\operatorname{sn}^2(u,k)+\operatorname{cn}^2(u,k)=1$$

$$k^2\operatorname{sn}^2(u,k)+\operatorname{dn}^2(u,k)=1$$

导数

$$\frac{d}{du}\operatorname{sn}(u,k)=\operatorname{cn}(u,k)\operatorname{dn}(u,k)$$

$$\frac{d}{du}\operatorname{cn}(u,k)=-\operatorname{sn}(u,k)\operatorname{dn}(u,k)$$

$$\frac{d}{du}\operatorname{dn}(u,k)=-k^2\operatorname{sn}(u,k)\operatorname{cn}(u,k)$$

是不是很像三角函数？当 $k=0$ 时 $\operatorname{sn}\to\sin$，$\operatorname{cn}\to\cos$，$\operatorname{dn}\to 1$。

例题

例 1：单摆大角度周期

单摆小角度周期 $T_0=2\pi\sqrt{\frac{L}{g}}$。大角度精确解：

$$T=4\sqrt{\frac{L}{g}}K\left(\sin\frac{\theta_0}{2}\right)$$

其中 $\theta_0$ 是最大摆角。所以第一类完全椭圆积分就是单摆的修正因子。

例 2：椭圆周长

椭圆 $\frac{x^2}{a^2}+\frac{y^2}{b^2}=1$，$a>b$，离心率 $e=\sqrt{1-\frac{b^2}{a^2}}$。

$$C=4aE(e)$$

如果 $a=5,b=3$，则 $e=\frac{4}{5}$，$C=20E(0.8)$。

例 3：数值验证

用 AGM 算 $K(0)$ 和 $E(0)$，应该得到 $\frac{\pi}{2}$。如果算出来不对，说明你 $b_0$ 写错了。

总结

1. 椭圆积分有三类：$F,E,\Pi$，分别是 $\frac{1}{\sqrt{\cdots}}$、$\sqrt{\cdots}$、$\frac{1}{(1-n\sin^2)\sqrt{\cdots}}$ 的积分。
2. 完全椭圆积分就是上限到 $\frac{\pi}{2}$。
3. 椭圆函数是椭圆积分的反函数，Jacobi sn/cn/dn 是三角函数的推广。
4. AGM 能高效算 $K$ 和 $E$。
5. 椭圆积分无处不在：椭圆周长、单摆、物理、密码学（某些椭圆曲线相关的东西）。

参考与推荐阅读

• OEIS 和各种手册里的椭圆积分表
• Abramowitz & Stegun: Handbook of Mathematical Functions
• 如果你想硬刚推导：Whittaker & Watson《A Course of Modern Analysis》

最后，如果你问椭圆积分有什么用，我只能说：它至少能算出椭圆到底有多长，而圆可以直接用 $2\pi r$，所以椭圆比圆更高级，没毛病。