自适应高斯-克朗罗德积分法在振荡函数积分中的变量替换技巧
我们先来看一个具体问题：如何计算积分
\[I = \int_{0}^{1} \cos(100x) \cdot e^{-x} \, dx \]
这个积分的特点是振荡频率高（\(\omega = 100\)），在区间 \([0,1]\) 上被积函数 \(\cos(100x)e^{-x}\) 震荡剧烈，直接用等距节点求积会需要极多节点才能捕捉振荡，效率很低。今天我们就讲一种结合变量替换与自适应高斯-克朗罗德积分的方法来高效计算此类振荡积分。
问题分析
-
直接使用自适应高斯-克朗罗德积分会遇到什么困难？
- 高斯-克朗罗德（Gauss-Kronrod）积分是在给定区间上用高阶（Gauss）节点与扩展的（Kronrod）节点分别计算积分值，通过比较两者差异来估计误差，从而自适应细分区间。
- 但对于高频振荡函数，在单个子区间内函数可能振荡多次，导致即使在高阶高斯节点上采样，也可能因为相位不对齐而完全错过振荡峰值与谷值，误差估计失效，最终被迫无限细分，效率极低。
-
变量替换的思路来源
- 振荡主要来源于 \(\cos(100x)\)，如果我们做一个变量替换 \(t = \phi(x)\)，使得在新变量 \(t\) 下，被积函数的振荡频率降低或者变得更容易采样。
- 一种常见技巧是采用线性频率映射：令 \( t = \omega x = 100x\)，则积分变为
\[ I = \frac{1}{100} \int_{0}^{100} \cos(t) \cdot e^{-t/100} \, dt \]
这时振荡部分 $\cos(t)$ 的频率是1（每 $2\pi$ 周期振荡一次），但积分区间长度变为100，仍然很长，直接数值积分节点数需求仍很大。
- 更好的想法：利用振荡函数的周期性，但这里乘以了衰减因子 \(e^{-x}\)，不是纯周期函数，所以需要更巧妙的变换。
方法设计：平稳相位型变量替换
对于积分 \(\int_a^b f(x) e^{i\omega g(x)} dx\)（我们的题目是实部情形，\(g(x)=x\)），当 \(\omega\) 很大时，主要贡献来自 \(g'(x)=0\) 的点（平稳相位点）。我们这里 \(g'(x)=1\) 没有零点，但我们可以通过变量替换改变相位函数的导数，使其在区间内变化平缓，从而降低有效振荡频率。
具体步骤：
- 选择变量替换函数 \(\phi\)
我们希望新变量 \(t = \phi(x)\) 满足 \(\phi'(x) \approx 1/\omega\) 在振荡剧烈的地方，从而 \(\cos(\omega x) = \cos(\omega \phi^{-1}(t))\) 在新变量下振荡变慢。
一个实用选择：令
\[ t = \frac{\arctan(\alpha x)}{\arctan(\alpha)}, \quad \alpha > 0 \]
这个映射将 \([0,1]\) 映射到 \([0,1]\)，且导数 \(\phi'(x) = \frac{\alpha}{(1+(\alpha x)^2)\arctan(\alpha)}\)。
当 \(\alpha\) 很大时，在 \(x\) 靠近0时导数很大，在 \(x\) 靠近1时导数很小。我们可以通过调节 \(\alpha\) 使得在函数振荡剧烈区域（整个区间）导数整体降低。
- 如何选取 \(\alpha\)
我们希望 \(\phi'(x)\) 在区间平均意义上与 \(1/\omega\) 相当，从而新变量下振荡频率降为 \(O(1)\)。
近似取平均导数：
\[ \bar{\phi}' = \frac{1}{1-0} \int_0^1 \phi'(x) dx = \frac{1}{\arctan(\alpha)} \int_0^1 \frac{\alpha}{1+(\alpha x)^2} dx = \frac{\arctan(\alpha)}{\arctan(\alpha)} = 1 \]
这提示我们直接使用该映射并不能降低平均振荡频率。因此我们需要更针对性的设计。
-
针对高频振荡的专用替换
更有效的方法：令新变量 \(t\) 直接正比于相位函数 \(g(x) = x\)，但压缩振荡区间。
考虑 \(t = \frac{x}{1 + (\omega-1)x}\)，这个函数在 \(x=0\) 时导数为1，在 \(x=1\) 时导数为 \(1/\omega\)，整体导数从1递减到 \(1/\omega\)，使得在区间后半段振荡被拉伸，从而在原变量下密集的振荡在新变量下变得稀疏。
但这个替换会导致积分表达式复杂，需要求逆函数 \(x = \frac{t}{1 - (\omega-1)t}\)，并计算 Jacobi 因子 \(dx/dt\)。 -
简化实用策略：分段线性频率缩放
将区间 \([0,1]\) 分成若干段，在每一段上用线性替换使得该段内相位变化为 \(O(1)\)。例如分成 \(m\) 段，每段长度 \(\Delta x = 1/m\)，在第 \(k\) 段做变换
\[ t = \omega (x - x_{k-1}) \quad (\text{局部频率归一化}) \]
但这本质上和直接对每个小区间用高斯积分没有区别，并没有减少所需区间数。
- 采用“振荡感知”的自适应策略
关键思想：在自适应高斯-克朗罗德积分中，不仅根据误差估计细分区间，还根据区间内振荡次数判断是否需要先进行变量替换。
算法步骤：- 对于当前区间 \([a,b]\)，计算相位变化量 \(\Delta \text{phase} = \omega \cdot (b-a)\)。
- 如果 \(\Delta \text{phase} > 2\pi\)（即区间内振荡超过一个周期），则先进行区间内的变量替换：令 \(u = \omega (x - a)\)，积分变为
\[ \int_a^b f(x) \cos(\omega x) dx = \frac{1}{\omega} \int_0^{\omega(b-a)} f(a+u/\omega) \cos(u+a\omega) du \]
此时被积函数振荡频率为1，但积分上限可能很大（若 $\omega(b-a)$ 很大）。我们继续对这个长区间应用自适应高斯-克朗罗德积分，但由于振荡频率已降为1，自适应效率会提高。
- 如果 \(\Delta \text{phase} \le 2\pi\)，则直接在该区间上应用普通的高斯-克朗罗德积分。
举例计算
以 \(I = \int_0^1 \cos(100x) e^{-x} dx\) 为例，精确值可用分部积分求出：
\[I = \frac{e^{-1}(100\sin 100 + \cos 100) - 1}{10001} \approx -0.0000976238 \]
步骤1：整个区间 \([0,1]\)，\(\Delta \text{phase} = 100 > 2\pi\)，因此做变量替换 \(u = 100x\)：
\[I = \frac{1}{100} \int_0^{100} \cos(u) e^{-u/100} du \]
现在被积函数 \(g(u) = \frac{1}{100} e^{-u/100} \cos(u)\)，振荡频率为1，但积分区间长度 \(L=100\) 很大。
步骤2：对新变量 \(u\) 的区间 \([0,100]\) 应用自适应高斯-克朗罗德积分。
- 由于指数衰减因子 \(e^{-u/100}\) 在 \(u\) 大时很小，积分主要贡献来自前几个周期，因此自适应算法会在 \([0, 10\pi] \approx [0,31.4]\) 内细分较多，而对 \(u>40\) 的区间，因函数值很小而快速收敛。
步骤3：在自适应积分过程中，每个子区间若仍包含超过一个振荡周期（即长度 \(>2\pi\)），可继续类似替换，但因为频率已为1，通常不需再替换。
步骤4：最终通过自适应高斯-克朗罗德积分得到近似值，与精确值比较误差。
误差控制与优势
- 优势：变量替换将高频振荡转化为低频振荡，使得高斯-克朗罗德积分在子区间上的误差估计更可靠，避免了因振荡过快而导致的误判。
- 误差控制：自适应过程中的误差估计基于 Kronrod 扩展节点与 Gauss 节点结果的差异。替换后函数更平滑，该差异更能真实反映积分误差，从而在达到指定容差时停止细分。
- 注意事项：变量替换会引入 Jacobi 因子 \(dx/du\)，需确保该因子在区间内光滑，否则可能引入新的奇异性。在我们的例子中，\(dx/du = 1/100\) 是常数，非常理想。
总结
对于高频振荡积分，直接应用自适应高斯-克朗罗德积分往往效率低下。通过先进行振荡感知的变量替换，将高频振荡转化为低频振荡，可以显著提高自适应积分的效率与可靠性。核心是检测区间内振荡周期数，当周期数过多时，通过线性相位映射降低频率，然后在新变量下继续自适应积分。这种方法结合了解析变换（降频）与数值自适应（局部误差控制）的优点，是处理高振荡积分的有效技巧之一。