椭圆积分——从不会到会 still 不会 index.md

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

前言(可跳过)

本文重点在于锻炼读者的"注意力"——不是让你背下更多公式,而是让你看见公式背后的结构

我自己检验是否真正理解一个公式的方法是:你能否用自然语言精确描述它给你的感觉?

比如 11x2dx=arcsinx\int \frac{1}{\sqrt{1-x^2}}\mathrm{d}x=\arcsin x,你看到的是"一个公式",还是"单位圆上角度与弧长的互换"?

椭圆积分之所以难,不是因为它公式复杂,而是因为它迫使你在"角度"和"弧长"之间建立一一对应——这种对应在圆里是平凡的(弧长=半径×角度),在椭圆里却需要全新的函数来描述。

本文会给出所有关键公式的几何意义证明思路。作者实力有限,可能无法道尽,希望知友评论区补充指正。

前置知识

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

其中第二条是很多人会忽略的,请读者务必最起码会使用定积分的几何意义。

约定俗成:下文用 EE 表示第二类椭圆积分、FF 表示第一类椭圆积分、Π\Pi 表示第三类。

引子:椭圆弧长——为什么需要新函数

你画一个椭圆

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

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

s=02πa2sin2θ+b2cos2θdθs=\int_0^{2\pi}\sqrt{a^2\sin^2\theta+b^2\cos^2\theta}\,\mathrm{d}\theta

化简一下,令离心率 e=1b2a2e=\sqrt{1-\frac{b^2}{a^2}}a>ba>b

s=4a0π21e2sin2θdθs=4a\int_0^{\frac{\pi}{2}}\sqrt{1-e^2\sin^2\theta}\,\mathrm{d}\theta

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

椭圆弧长示意图 (图源:AI生成,模拟手绘风格)

最重要的直觉

  1. 圆的弧长公式 s=rθs=r\theta 之所以简单,是因为圆的曲率处处相等——你转过的角度和走过的弧长成严格正比
  2. 椭圆的曲率在变化——靠近长轴端点曲率大,靠近短轴端点曲率小
  3. 所以"转过的角度"和"走过的弧长"不再成正比,E(e)E(e) 就是对这种非线性关系的精确编码

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

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

定义与几何意义

F(φ,k)=0φdθ1k2sin2θF(\varphi,k)=\int_0^{\varphi}\frac{\mathrm{d}\theta}{\sqrt{1-k^2\sin^2\theta}}

其中 k[0,1)k\in[0,1)模数φ\varphi振幅

几何意义F(φ,k)F(\varphi,k)椭圆上的"归一化弧长参数"

具体地,考虑椭圆 x=asinθx=a\sin\thetay=bcosθy=b\cos\theta。从 θ=0\theta=0θ=φ\theta=\varphi 的弧长为

s(φ)=a0φ1e2sin2θdθ=aE(φ,e)s(\varphi)=a\int_0^{\varphi}\sqrt{1-e^2\sin^2\theta}\,\mathrm{d}\theta=aE(\varphi,e)

F(φ,k)F(\varphi,k)E(φ,k)E(\varphi,k) 通过反函数关系耦合:如果 u=F(φ,k)u=F(\varphi,k),那么 φ\varphi 就是 uu振幅,记作 φ=am(u,k)\varphi=\operatorname{am}(u,k)

为什么 φ=π2\varphi=\frac{\pi}{2} 叫"完全"

因为 F(φ,k)F(\varphi,k) 关于 φ\varphi周期性的,周期为 4K(k)4K(k)K(k)K(k) 是完全椭圆积分)。当 φ\varphi00 走到 π2\frac{\pi}{2} 时,点在椭圆上从短轴端点走到长轴端点,恰好遍历了椭圆的四分之一周长。这就是为什么 φ=π2\varphi=\frac{\pi}{2} 时叫"完全"——它对应着四分之一周期,而不是整个周期。

φ=π2\varphi=\frac{\pi}{2} 时,记作 K(k)K(k)

K(k)=0π2dθ1k2sin2θK(k)=\int_0^{\frac{\pi}{2}}\frac{\mathrm{d}\theta}{\sqrt{1-k^2\sin^2\theta}}

特殊值

  1. k=0k=0K(0)=π2K(0)=\frac{\pi}{2},退化回普通积分。
  2. k1k\to 1^-K(k)+K(k)\to +\infty,因为 1sin2θ1-\sin^2\thetaθ=π2\theta=\frac{\pi}{2} 附近变成 00,分母爆炸。

它积出来是啥?

很不幸,第一类椭圆积分没有初等函数形式。就像 ex2dx\int e^{-x^2}\,\mathrm{d}x 积不出来要叫 erf\operatorname{erf} 一样,椭圆积分也要被赐名为 FFKK

但你可以查表,或者用 AGM 算法(后面讲)。

直觉F(φ,k)F(\varphi,k) 衡量的是**"角度被压缩的程度"**。kk 越大,sin2θ\sin^2\theta 越"有效",分母越小,同样的 dθ\mathrm{d}\theta 对应的 FF 增量越大。

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

定义:

E(φ,k)=0φ1k2sin2θdθE(\varphi,k)=\int_0^{\varphi}\sqrt{1-k^2\sin^2\theta}\,\mathrm{d}\theta

完全形式:

E(k)=0π21k2sin2θdθE(k)=\int_0^{\frac{\pi}{2}}\sqrt{1-k^2\sin^2\theta}\,\mathrm{d}\theta

特殊值

  1. k=0k=0E(0)=π2E(0)=\frac{\pi}{2}
  2. k=1k=1E(1)=1E(1)=1(因为 0π2cosθdθ=1\int_0^{\frac{\pi}{2}}\cos\theta\,\mathrm{d}\theta = 1)。

椭圆周长公式

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

直觉E(φ,k)E(\varphi,k)F(φ,k)F(\varphi,k)对偶的。FF 的被积函数在分母,EE 的在分子。kk 增大时,FF 发散(分母趋零),EE 收敛到 11(分子趋零)。

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

定义:

Π(n;φ,k)=0φdθ(1nsin2θ)1k2sin2θ\Pi(n;\varphi,k)=\int_0^{\varphi}\frac{\mathrm{d}\theta}{(1-n\sin^2\theta)\sqrt{1-k^2\sin^2\theta}}

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

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

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

它们之间的关系——Legendre 关系

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

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

这个式子为什么成立?

证明思路:Legendre 关系可以通过对 K(k)K(k)E(k)E(k) 分别求导,利用它们的微分方程联立消元得到。

具体地,我们知道 dKdk=E(k)(1k2)K(k)k(1k2),dEdk=E(k)K(k)k\frac{\mathrm{d}K}{\mathrm{d}k}=\frac{E(k)-(1-k^2)K(k)}{k(1-k^2)},\qquad\frac{\mathrm{d}E}{\mathrm{d}k}=\frac{E(k)-K(k)}{k}

考虑函数 f(k)=E(k)K(k)+E(k)K(k)K(k)K(k)f(k)=E(k)K(k')+E(k')K(k)-K(k)K(k')

f(k)f(k) 求导,利用 k=1k2k'=\sqrt{1-k^2} 和链式法则(dkdk=kk\frac{\mathrm{d}k'}{\mathrm{d}k}=-\frac{k}{k'}),代入上述导数公式,经过一系列化简(约一页纸),可以发现 dfdk=0\frac{\mathrm{d}f}{\mathrm{d}k}=0

所以 f(k)f(k) 是常数。代入 k=0k=0(此时 k=1k'=1K(0)=E(0)=π2K(0)=E(0)=\frac{\pi}{2}K(1)K(1)\to\inftyE(1)=1E(1)=1 使乘积有限),计算得 f(0)=π2f(0)=\frac{\pi}{2}

因此 f(k)=π2f(k)=\frac{\pi}{2} 对所有 kk 成立。

直觉:Legendre 关系是椭圆积分周期性的体现K(k)K(k)E(k)E(k) 分别衡量椭圆积分的"第一类"和"第二类"行为,而 kk' 代表"正交方向"的模数。这个关系保证了无论 kk 怎么变,KKEE 的乘积组合始终保持不变——就像能量守恒一样。

怎么记住它

  1. 三项轮换:EK+EKKKE K' + E' K - K K'
  2. 常数项:π2\frac{\pi}{2}(来自 k=0k=0 时的平凡值)
  3. 对称性:kkk\leftrightarrow k' 时式子不变

求导和恒等式

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

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

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

直觉KK 的导数包含 EEEE 的导数包含 KK。它们互相定义,就像 sin\sincos\cos 的导数互相包含对方一样。

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

K(k)K(k) 可以用算术-几何平均(AGM)快速算:

a0=1a_0=1b0=1k2=kb_0=\sqrt{1-k^2}=k',然后迭代:

an+1=an+bn2a_{n+1}=\frac{a_n+b_n}{2}

bn+1=anbnb_{n+1}=\sqrt{a_nb_n}

cn+1=anbn2c_{n+1}=\frac{a_n-b_n}{2}

nn 足够大时,anbnM(a0,b0)a_n\approx b_n\approx M(a_0,b_0),AGM 收敛非常快,每轮有效位数翻倍。

第一类完全椭圆积分:

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

第二类完全椭圆积分:

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

为什么 AGM 能算出椭圆积分?

证明思路:AGM 与椭圆积分的联系来自Landen 变换。Landen 变换是一种椭圆积分的变量替换,它将 F(φ,k)F(\varphi,k) 转化为另一个椭圆积分,其模数 k1=1k1+kk_1=\frac{1-k'}{1+k'} 比原来的 kk 更小。反复应用 Landen 变换,模数趋近于 00,积分趋近于 π2\frac{\pi}{2}

AGM 正是 Landen 变换的代数形式:an+1=an+bn2a_{n+1}=\frac{a_n+b_n}{2}bn+1=anbnb_{n+1}=\sqrt{a_nb_n} 分别对应 Landen 变换中的算术和几何操作。当 ana_nbnb_n 收敛到同一个值 MM 时,积分值被"压缩"成了 π2M\frac{\pi}{2M}

AGM 迭代与单摆 (图源:AI生成,模拟手绘风格)

最重要的直觉

  1. AGM 是算术平均和几何平均的迭代,两者收敛到同一个值
  2. 每轮迭代,cnc_n(差的一半)平方级减小,所以收敛极快
  3. K(k)K(k) 就是 π2\frac{\pi}{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;
}

坑点:

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

椭圆函数 Jacobi sn, cn, dn

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

定义:如果

u=F(φ,k)=0φdθ1k2sin2θu=F(\varphi,k)=\int_0^{\varphi}\frac{\mathrm{d}\theta}{\sqrt{1-k^2\sin^2\theta}}

那么 Jacobi 椭圆函数:

sn(u,k)=sinφ\operatorname{sn}(u,k)=\sin\varphi

cn(u,k)=cosφ\operatorname{cn}(u,k)=\cos\varphi

dn(u,k)=1k2sin2φ\operatorname{dn}(u,k)=\sqrt{1-k^2\sin^2\varphi}

椭圆函数与三角函数的关系 (图源:AI生成,模拟手绘风格)

基本恒等式

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

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

导数

ddusn(u,k)=cn(u,k)dn(u,k)\frac{\mathrm{d}}{\mathrm{d}u}\operatorname{sn}(u,k)=\operatorname{cn}(u,k)\operatorname{dn}(u,k)

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

ddudn(u,k)=k2sn(u,k)cn(u,k)\frac{\mathrm{d}}{\mathrm{d}u}\operatorname{dn}(u,k)=-k^2\operatorname{sn}(u,k)\operatorname{cn}(u,k)

k=0k=0snsin\operatorname{sn}\to\sincncos\operatorname{cn}\to\cosdn1\operatorname{dn}\to 1

最重要的直觉

  1. sn\operatorname{sn}cn\operatorname{cn}sin\sincos\cos推广——单位圆被"压扁"成了椭圆
  2. dn\operatorname{dn} 是椭圆独有的,衡量"压扁程度"
  3. 导数公式里 dn\operatorname{dn} 的出现,是因为椭圆上"角度"和"弧长"的关系不再像圆那样简单

例题

例 1:单摆大角度周期

单摆小角度周期 T0=2πLgT_0=2\pi\sqrt{\frac{L}{g}}。大角度精确解:

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

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

直觉:摆角越大,摆锤在最高点附近"停留"的时间越长(因为势能曲线变平),周期变长。K(sinθ02)K(\sin\frac{\theta_0}{2}) 就是对这个"停留效应"的精确量化。

例 2:椭圆周长

椭圆 x2a2+y2b2=1\frac{x^2}{a^2}+\frac{y^2}{b^2}=1a>ba>b,离心率 e=1b2a2e=\sqrt{1-\frac{b^2}{a^2}}

C=4aE(e)C=4aE(e)

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

例 3:数值验证

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

大坑

注意一个很神奇的问题,比如上面 AGM 代码,你可能会发现我将

ld sum = 0.0;
ld pw2 = 1.0; // 2^n

注释了个"2^n",为什么呢?

举个例子,你把它写成了 2n12^{n-1},代码变成了下面这样。

sum += pw2 * cn * cn; // pw2 = 2^n

可能乍一看没问题,毕竟系数差个 2 而已。

BUT,当你把目光放到 E(k)E(k) 的公式那部分的时候,你就发现

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

n=0n=0 时,2n1=122^{n-1}=\frac{1}{2},所以 c02c_0^2 的系数是 12\frac{1}{2},不是 11!!!!!

如果你代码里写成了 2n2^n,那么 c02c_0^2 的系数就是 11E(k)E(k) 就会多减 12c02\frac{1}{2}c_0^2,直接导致结果错误!!!

正确代码sum += pw2 * cn * cn 之后 pw2 *= 2.0,这样第 nn 次迭代时 pw2 正好是 2n2^n,但公式里是 2n12^{n-1},所以最后要 (1.0 - sum / 2.0) 来补偿。

总结

  1. 椭圆积分有三类:F,E,ΠF,E,\Pi,分别是 1\frac{1}{\sqrt{\cdots}}\sqrt{\cdots}1(1nsin2)\frac{1}{(1-n\sin^2)\sqrt{\cdots}} 的积分。
  2. 完全椭圆积分就是上限到 π2\frac{\pi}{2},对应椭圆的四分之一周期。
  3. 椭圆函数是椭圆积分的反函数,Jacobi sn/cn/dn 是三角函数的推广。
  4. AGM 能高效算 KKEE,其原理是 Landen 变换。
  5. Legendre 关系是椭圆积分周期性的体现,可以通过微分方程证明。

在做微积分的时候,看到巧合式子,一定要多想。想明白为什么,会给你的数学直觉带来很大的帮助。

注意力与本质 (图源:AI生成,模拟手绘风格)

参考与推荐阅读