數值分區間 數值積分 使用FFT來降低計算量

2021-10-16 04:22:04 字數 3228 閱讀 6578

我們考慮如下積分的數值計算問題

我們假定

是有界的,因為

的快速下降性,我們知道上述積分的主要貢獻來自於

這樣的乙個鄰域(當然,

只是乙個示例。)

因此,我們可以採用復合梯形公式來數值求解積分。

取為等距節點,則復合梯形公式給出數值積分近似

下面我們以

以及 為例來計算一下,首先這個積分的精確值為

數值上,我們分別改變量值積分節點個數

來看一下誤差變化

**如下

%% 被積函式

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));

endfigure;

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變換中的方法計算

在的值。我們選取

,並且取

作為在數值積分時的節點

我們將比較一下數值積分和精確值

在不同

時候的誤差。

當然,本身

在增大的時候,精確值就變得更小了。

不過這裡的重點在於

值得注意的一點是,在利用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);

endfigure;

semilogy(alpha_set,err,'-o','linewidth',2,'markersize',8);

python做數值積分 積分區間為小數 數值積分

一.組合梯形公式 數值積分 問題描述 組合梯形公式求函式f x 2 sin 2 sqrt x 的積分近似值。輸入形式 在螢幕上依次輸入積分上限 下限和等距子區間個數。輸出形式 輸出使用組合梯形公式求得的積分近似值。樣例1輸入 1 6 10 樣例1輸出 8.19385457 樣例1說明 輸入 積分上限...

數值分區間 置信區間與參考值範圍

誤差永遠存在,而且不可避免。即使實驗條件再精確也無法完全避免隨機干擾的影響,所以做科學實驗往往要測量多次,用取平均值之類的統計手段去得出結果。多次測量,是乙個排除偶然因素的好辦法。如國足輸掉比賽之後經常抱怨偶然因素,有時候是因為裁判不公,有時候是因為主力不在,有時候是因為不適應客場氣候,關鍵是如果你...

matlab進行數值積分的主要函式使用方法

matlab進行數值積分的主要函式 1.trapz 梯形法求解積分 x 0 pi 10 pi y sin x trapz x,y 2.quad 基於變步長simpso法求積分 q quad fun,a,b,tol 其中fun是被積函式檔名或函式控制代碼,a,b是積分下限和積分上限,tol是積分精度 ...