fft 相位谱_数值积分——使用FFT来降低计算量

问题引入：

我们考虑如下积分的数值计算问题

我们假定

是有界的，因为

的快速下降性，我们知道上述积分的主要贡献来自于

这样的一个邻域（当然，

只是一个示例。）

因此，我们可以采用复合梯形公式来数值求解积分。

取

为等距节点，则复合梯形公式给出数值积分近似

下面我们以

以及

为例来计算一下，首先这个积分的精确值为

数值上，我们分别改变数值积分节点个数

来看一下误差变化

代码如下

%% 被积函数
func = @(x) (x+1).^2.*exp(-x.^2);

m = 10;
err = zeros(1,m);
num_set = 4+3*(1:m);
for ii = 1:m
    N = num_set(ii);
    x = linspace(-10,10,N+1);
    % 步长
    h = mean(diff(x));
    Num_int = (func(x(1))/2+sum(func(x(2:end-1)))+func(x(end))/2)*h;
    err(ii) = abs(Num_int - 3/2* sqrt(pi));
end

figure; 
semilogy(num_set,err,'-o','LineWidth',2,'MarkerSize',10)

可以看到，当积分节点个数

增长的时候，误差在极速下降——这是因为对于光滑的周期函数，复合梯形公式收敛是谱精度的。

问题：如果我们不仅仅是对一个

进行计算，而是对一系列都要计算。

那么每次计算都需要对 f 进行重新的采样，如果函数 f 的计算很复杂的话，这会带来很大的计算量。

我们希望关于 f 的计算，只用做一次。这就可以利用 Fourier transform 重写积分形式。

进化版本

利用Fourier变换以及Paserval等式，我们有

其中

首先我们看一下

的傅里叶变换

因此

重复一下“优点”

利用Fourier变换得到的积分的形式相比于之前版本的好处：

在我们要对很多不同的

进行积分求解的时候，之前版本需要对每个

都要计算许多

的值

但是新的形式绕过了这一点。

这一点尤其在

的计算量很大的时候显得格外重要。

实验

和之前的版本相比，高斯函数的波包变宽了，变成了

我们仍旧近似使用

作为支撑集。

只不过这里不能再直接使用上面

的例子了，它不是

可积的，不好求傅里叶变换。

不过我们可以在最开始处理的时候，写成

给

配上一个指数衰减的因子进行调制。不过这还是会对我们后续的论述带来麻烦（Fourier变换不是那么简单）。

因此我们选取

则

上述结果对

是个复数也成立，使用留数定理证明。

新问题

我们遇到一个新问题，在对下面的积分数值求解的过程中

我们会像之前使用复合梯形公式那样，采用等距节点

那么问题来了，我们如何计算

通过 Matlab利用fft计算Fourier变换可以知道这可以通过 fft 以及一个相位因子实现。

实验

因为

所以我们在计算它的Fourier变换时可以将区间截断成

足以。

这样，我们通过Matlab利用fft计算Fourier变换中的方法计算

在

的值。

我们选取

，并且取

作为

在数值积分时的节点

我们将比较一下数值积分和精确值

在不同

时候的误差。

当然，本身

在增大的时候，精确值就变得更小了。

不过这里的重点在于

• 即使
  较小的时候，误差也有 1e-6 的量级
• 该计算方法不过分依赖于 f 的取值运算

值得注意的一点是，在利用 fft 计算出Fourier变换在一些点上的值之后，使用样条插值的结果要远好于线性插值的结果。

代码如下：

g = @(x) exp(-x.^2);
c = -20; d = 20;
N = 80;
x = linspace(c,d,N+1);
x = x(1:end-1);

xx = 2*pi/(d-c)*([0:N/2 -N/2+1:-1]);
xxshift = [xx(N/2+2:end) xx(1:N/2+1)];

%% 计算g的傅里叶变换在对应点的值的近似
gxx = exp(-1i*c*xx)*(d-c)/N.*fft(g(x));
gxxshift = [gxx(N/2+2:end) gxx(1:N/2+1)];

% 积分节点和插值近似
xx_inter = linspace(-10,10,256);
h = mean(diff(xx_inter));
gxx_inter = interp1(xxshift,gxxshift,xx_inter,'spline','extrap');

alpha_set = 0:20;
err = zeros(size(alpha_set));
for ii = 1:length(alpha_set)
    alpha = alpha_set(ii);
    phase = exp(1i*alpha*xx_inter);
    integrand = phase.*gxx_inter.*exp(-xx_inter.^2/4);
    exact_val = exp(-alpha^2/2)*sqrt(pi/2);
    % 梯形公式
    Num_val = ((integrand(1)+integrand(end))/2+sum(integrand(2:end-1)))*h/(2*sqrt(pi));
    err(ii) = abs(Num_val-exact_val);
end

figure; 
semilogy(alpha_set,err,'-o','LineWidth',2,'MarkerSize',8);