问题引入:
我们考虑如下积分的数值计算问题
我们假定
这样的一个邻域(当然,
因此,我们可以采用复合梯形公式来数值求解积分。
![36f4c9d3652bdd2591ec3694f6a3db1a.png](https://i-blog.csdnimg.cn/blog_migrate/65ef05d19157f28f126ad83d794f0c94.png)
取
为等距节点,则复合梯形公式给出数值积分近似
下面我们以
数值上,我们分别改变数值积分节点个数
![3339f85d313ad1825b66e78894ae0f6d.png](https://i-blog.csdnimg.cn/blog_migrate/c30949237b03cd0ab5d740ab3e0bc0e5.jpeg)
代码如下
%% 被积函数
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变换得到的积分的形式相比于之前版本的好处:
在我们要对很多不同的
但是新的形式绕过了这一点。
这一点尤其在
实验
和之前的版本相比,高斯函数的波包变宽了,变成了
我们仍旧近似使用
只不过这里不能再直接使用上面
不过我们可以在最开始处理的时候,写成
给
因此我们选取
则
上述结果对
新问题
我们遇到一个新问题,在对下面的积分数值求解的过程中
我们会像之前使用复合梯形公式那样,采用等距节点
那么问题来了,我们如何计算
通过 Matlab利用fft计算Fourier变换可以知道这可以通过 fft
以及一个相位因子实现。
实验
因为
所以我们在计算它的Fourier变换时可以将区间截断成
这样,我们通过Matlab利用fft计算Fourier变换中的方法计算
的值。
我们选取
在数值积分时的节点
我们将比较一下数值积分和精确值
在不同
![fb2f68c2c321e9b311ead737afa20ac1.png](https://i-blog.csdnimg.cn/blog_migrate/e4952b53241da13d142c6127faa7c079.jpeg)
当然,本身
不过这里的重点在于
- 即使
较小的时候,误差也有 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);