我們考慮如下積分的數值計算問題
我們假定
是有界的,因為
的快速下降性,我們知道上述積分的主要貢獻來自於
這樣的乙個鄰域(當然,
只是乙個示例。)
因此,我們可以採用復合梯形公式來數值求解積分。
取為等距節點,則復合梯形公式給出數值積分近似
下面我們以
以及 為例來計算一下,首先這個積分的精確值為
數值上,我們分別改變量值積分節點個數
來看一下誤差變化
**如下
%% 被積函式
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是積分精度 ...