龍貝格積分算法實驗

1、實驗題目2 Romberg積分法摘要 考慮積分欲求其近似值，可以采用如下公式： （復(fù)化）梯形公式 （復(fù)化）辛卜生公式 （復(fù)化）柯特斯公式 這里，梯形公式顯得算法簡單，具有如下遞推關(guān)系因此，很容易實現(xiàn)從低階的計算結(jié)果推算出高階的近似值，而只需要花費較少的附加函數(shù)計算。但是，由于梯形公式收斂階較低，收斂速度緩慢。所以，如何提高收斂速度，自然是人們極為關(guān)心的課題。為此，記為將區(qū)間進行等份的復(fù)化梯形積分結(jié)果，為將區(qū)間進行等份的復(fù)化辛卜生積分結(jié)果，為將區(qū)間進行等份的復(fù)化柯特斯積分結(jié)果。根據(jù)李查遜（Richardson）外推加速方法，可得到 可以證明，如果充分光滑，則有 （固定） 這是一個收斂速度更快的一

2、個數(shù)值求積公式，我們稱為龍貝格積分法。 該方法的計算可按下表進行 很明顯，龍貝格計算過程在計算機上實現(xiàn)時，只需開辟一個一維數(shù)組，即每次計算的結(jié)果，可存放在位置上，其最終結(jié)果是存放在位置上。前言利用龍貝格積分法計算積分程序設(shè)計流程： 1準(zhǔn)備初值，計算且（為等份次數(shù)） 2按梯形公式的遞推關(guān)系，計算 3按龍貝格公式計算加速值 4精度控制。對給定的精度，若則終止計算，并取作為所求結(jié)果；否則，重復(fù)24步，直到滿足精度為止。問題1（1） 程序運行如下：I = Romberginterg(inline(x.2.*exp(x),0,1,25,1e-6)I = 0.7183（2） 程序運行如下：I = Romb

3、erginterg(inline(exp(x).*sin(x),1,3,25,1e-6)I = 10.9502（3） 程序運行如下：I = Romberginterg(inline(4./(1+x.2),0,1,25,1e-6)I = 3.1416（4） 程序運行如下：I = Romberginterg(inline(1./(1+x),0,1,25,1e-6)I = 0.6931問題2（1） 程序運行如下：I = Romberginterg(inline(cos(x)./sqrt(1-x),0,1,25,1e-6)I = NaN因為函數(shù)在x=0處出現(xiàn)了0/0的情況，極限為1，所以Matlab的

4、結(jié)果顯示為NaN非數(shù)，解決方法是把下限0改為一個小數(shù)eps。I = Romberginterg(inline(sin(x)./x),eps,1,25,1e-6)I = 0.9461（2） 程序運行如下：I = Romberginterg(inline(cos(x)./sqrt(1-x),0,1,25,1e-6)I = NaN與（1）的原因相同，函數(shù)在x=1處出現(xiàn)了0/0的情況，結(jié)果為無窮，此時應(yīng)選擇變換，將積分變?yōu)?，再進行變換，將積分變?yōu)椋儞Q后輸入命令：I = RombergInterg(inline(2*sin(x).*cos(sin(x).2),0,pi/2,25,1e-6)I = 1.

5、4996（3） 程序運行如下：I = Romberginterg(inline(cos(x)./sqrt(x),0,1,25,1e-6)I = NaN函數(shù)在x=0處出現(xiàn)了1/0的情況，結(jié)果為無窮。先將積分變?yōu)?，再做變換，I = 2*Romberginterg(inline(cos(x.2),0,1,25,1e-6)I = 1.8090（4） 程序運行如下：I = Romberginterg(inline(x.*sin(x)./(1-x.2),-1,1,25,1e-6)I = NaN函數(shù)在x=1,-1處出現(xiàn)了sin1/0的情況，結(jié)果為無窮。被積函數(shù)為偶函數(shù)，做變換，積分變?yōu)镮 = 2*Rombe

6、rginterg(inline(sin(x).*sin(sin(x),0,pi/2,25,1e-6)I = 1.3825或者使用Gauss-Chebyshev求積公式：I = GaussInterg(inline(x.*sin(x),Chebyshev,-1,1,1e-6)I = 1.3825由于Gauss-Chebyshev求積公式求的是在的積分因此此處的為。問題3（1） 程序運行如下：I = Romberginterg(inline(exp(-x.2),-10,10,25,1e-6)I = 1.7725由于積分區(qū)間無限，而函數(shù)在-10,10外很小，可以用函數(shù)在-10,10上的積分近似這個無

7、窮積分?；蛘呤褂肎auss-Hermite求積公式：I = GaussInterg(inline(x+1-x),Hermite,0,0,1e-6)I = 1.7725由于Gauss-Hermite求積公式求的是在的積分因此此處的為1。（2） 程序運行如下：先做變換，則原積分變?yōu)镮 = 2*Romberginterg(inline(1./(1+x.2),0,1,25,1e-6)I = 1.5708（3） 程序運行如下：使用Gauss-Hermite求積公式：I = GaussInterg(inline(cos(x).3),Hermite,0,0,1e-6)I = 1.0820由于Gauss-He

8、rmite求積公式求的是在的積分因此此處的為。（4） 程序運行如下：使用Gauss-Laguerre求積公式：I = GaussInterg(inline(sin(x).2),Laguerre,0,0,1e-6)I = 0.4000由于Gauss-Laguerre求積公式求的是在的積分因此此處的為。使用的函數(shù)function I = RombergInterg(fun, a, b, npanel, tol, flag)% RombergInterg 用Romberg方法求積分% Synopsis: I = RombergInteg(fun,a,b,n,tol)% Input: fun = (s

9、tring) 被積函數(shù)的函數(shù)名% a, b = 積分下限和積分上限% npanel = (optional) 將積分區(qū)間平分的段數(shù)，默認(rèn)為25% tol = (optional) 計算誤差上限，默認(rèn)為5e-9% flag = (optional) 是否顯示Romberg-T數(shù)表，默認(rèn)為0為不顯示% Output: I = 通過Romberg方法求積分的近似值 if nargin 4 npanel = 25; end if nargin 5 tol = 5e-9; end if nargin = tol T(m,1) = TrapezoidInteg(fun, a, b, 2m*npanel);

10、 %計算第m行T0(h/2m) T(m,m) = 0; %將矩陣T變成m*m for n = 2:m T(m,n) = ( 4(n-1)*T(m,n-1) - T(m-1,n-1) / (4(n-1) - 1); %Tm(h/2k)與Tm-1(h/2(k+1)和Tm-1(h/2k)的遞推關(guān)系 end err = abs( T(m,m) - T(m-1,m-1) ); %計算誤差值 m = m + 1; %計算下一行 end I = T(m-1,m-1); if flag = 0 disp(T); endfunction I = TrapezoidInteg(fun, a, b, npanel)

11、% TrapezoidInteg 用復(fù)化梯形公式求積分% Synopsis: I = TrapezoidInteg(fun,a,b,n)% Input: fun = (string) 被積函數(shù)的函數(shù)名% a, b = 積分下限和積分上限% npanel = (optional) 將積分區(qū)間平分的段數(shù)，默認(rèn)為25% Output: I = 通過復(fù)化梯形公式求積分的近似值 if nargin 4 npanel = 25; end nnode = npanel + 1; %節(jié)點數(shù) = 段數(shù) + 1 h = (b-a)/(nnode-1); %步長 x = a:h:b; %將積分區(qū)間分段 f = fe

12、val(fun,x); %求節(jié)點處被積函數(shù)的值 I = h * ( 0.5*f(1) + sum(f(2:nnode-1) + 0.5*f(nnode) );function I = GaussInterg(fun, type, a, b, tol)% GaussInterg 用Gauss型求積公式求積分,具體形式由使用者選取% Synopsis: I = GaussInterg(fun, type, a, b)% I = GaussInterg(fun, type, a, b, tol)% Input: fun = (string) 被積函數(shù)的函數(shù)名% type = (string) 具體G

13、auss求積公式的形式% a,b = 積分上下限，Laguerre只計算0到inf，Hermite只計算-inf到inf，所以a,b對這兩種形式無效% tol = (optional) 誤差容忍限度，默認(rèn)為5e-5% Output: I = 通過Gauss型求積公式求積分的近似值 if nargin = tol switch type case Legendre %計算無奇點被積函數(shù)在-1到1的積分 I = GaussLegendreInterg(fun, a, b, n); case Chebyshev %計算被積函數(shù)形如 f(x)/sqrt(1-x2)在-1到1的積分 I = GaussC

14、hebyshevInterg(fun, a, b, n); case Laguerre %計算被積函數(shù)形如 exp(-x)*f(x)的在0到inf的積分 I = GaussLaguerreInterg(fun, n); case Hermite %計算被積函數(shù)形如 exp(-x2)*f(x)的在-inf到inf的積分 I = GaussHermiteInterg(fun, n); otherwise error(No such type!); end err = abs(I - IOld); %計算誤差 IOld = I; %把IOld賦值為I進行下次迭代 n = n+1; %多項式節(jié)點遞增

15、endfunction I = GaussLegendreInterg(fun, a, b, n)% GaussLegendreInterg 用Gauss-Legendre型求積公式計算無奇點被積函數(shù)在-1到1的積分% Synopsis: I = GaussLegendreInterg(fun, a, b)% Input: fun = (string) 被積函數(shù)的函數(shù)名% a,b = 積分下限和積分上限% n = (optional) 高斯節(jié)點的個數(shù)，默認(rèn)為7% Output: I = 通過Gauss型求積公式求積分的近似值 if nargin 4 n = 7; end if round(n)

16、 = n | n 1 error(Gauss節(jié)點個數(shù)必須為正整數(shù)); end p = LegendreIter(n); %構(gòu)造n階勒讓德多項式 p = p/polyval(p,1); %使Pn(1) = 1; for i = 1:n+1 %構(gòu)造p的導(dǎo)數(shù)p1 p1(i) = (n-i+1)*p(i); end p1(n+1) = ; r = roots(p); %計算勒讓德多項式的根 L = polyval(p1,r); %計算與勒讓德多項式的根匹配的高斯積分系數(shù) A = 2./(1-r.2)./L.2; xi = ( (b-a)*r + (a+b) ) / 2; %找出勒讓德多項式的根在a,b

17、上對應(yīng)的數(shù)值xi yi = feval(fun, xi); %計算f(xi) I = (b-a)/2 * (A*yi); %I = sum(Ai*yi),前面的系數(shù)是積分區(qū)間改變是增加的常數(shù) function I = GaussChebyshevInterg(fun, a, b, n)% GaussChebyshevInterg 用Gauss-Chebyshev型求積公式求積分計算被積函數(shù)形如 f(x)/sqrt(1-x2)在-1到1的積分% Synopsis: I = GaussChebyshevInterg(fun, a, b)% Input: fun = (string) 被積函數(shù)的函數(shù)

18、名% a,b = 積分下限和積分上限% n = (optional) 高斯節(jié)點的個數(shù)，默認(rèn)為7% Output: I = 通過Gauss型求積公式求積分的近似值 if nargin 4 n = 7; end if round(n) = n | n 1 error(Gauss節(jié)點個數(shù)必須為正整數(shù)); end p = ChebyshevIter(n); %構(gòu)造n階切比雪夫多項式 A = pi/n * ones(n,1); %計算與切比雪夫多項式的根匹配的高斯積分系數(shù) r = roots(p); %計算切比雪夫多項式的根 xi = ( (b-a)*r + (a+b) ) / 2; %找出切比雪夫多項

19、式的根在a,b上對應(yīng)的數(shù)值xi yi = feval(fun, xi); %計算f(xi) I = (b-a)/2 * (A*yi); %I = sum(Ai*yi),前面的系數(shù)是積分區(qū)間改變是增加的常數(shù) function I = GaussLaguerreInterg(fun, n)% GaussLaguerreInterg 用Gauss-Laguerre型求積公式求積分計算被積函數(shù)形如 exp(-x)*f(x)的在0到inf的積分% Synopsis: I = GaussLaguerreInterg(fun)% I = GaussLaguerreInterg(fun, n)% Input:

20、 fun = (string) 被積函數(shù)的函數(shù)名% n = (optional) 高斯節(jié)點的個數(shù)，默認(rèn)為7% Output: I = 通過Gauss型求積公式求積分的近似值 if nargin 2 n = 7; end if round(n) = n | n 1 error(Gauss節(jié)點個數(shù)必須為正整數(shù)); end p = LaguerreIter(n); %構(gòu)造n階拉蓋爾多項式 p1 = LaguerreIter(n+1); %構(gòu)造n+1階拉蓋爾多項式 xi = roots(p); %計算拉蓋爾多項式的根 L = polyval(p1,xi); %計算與拉蓋爾多項式的根匹配的高斯積分系數(shù)

21、A = xi./(n+1)2*L.2); yi = feval(fun, xi); %計算f(xi) I = A*yi; %I = sum(Ai*yi),前面的系數(shù)是積分區(qū)間改變是增加的常數(shù) function I = GaussHermiteInterg(fun, n)% GaussHermiteInterg 用Gauss-Hermite型求積公式求積分計算被積函數(shù)形如 exp(-x2)*f(x)的在-inf到inf的積分% Synopsis: I = GaussHermiteInterg(fun)% I = GaussHermiteInterg(fun, n)% Input: fun = (

22、string) 被積函數(shù)的函數(shù)名% n = (optional) 高斯節(jié)點的個數(shù)，默認(rèn)為7% Output: I = 通過Gauss型求積公式求積分的近似值 if nargin 2 n = 7; end if round(n) = n | n 1 error(Gauss節(jié)點個數(shù)必須為正整數(shù)); end p = HermiteIter(n); %構(gòu)造n階埃爾米特多項式 p1 = HermiteIter(n-1); %構(gòu)造n-1階埃爾米特多項式 xi = roots(p); %計算埃爾米特多項式的根 L = polyval(p1(1:n),xi); %計算與埃爾米特多項式的根匹配的高斯積分系數(shù) A

23、 = 2(n-1)*sqrt(pi)*factorial(n) ./ (n2 * L.2); yi = feval(fun, xi); %計算f(xi) I = A*yi; %I = sum(Ai*yi),前面的系數(shù)是積分區(qū)間改變是增加的常數(shù) function p = LegendreIter(n)% LegendreIter 用遞推的方法計算n次勒讓德多項式的系數(shù)向量 Pn+2(x) = (2*i+3)/(i+2) * x*Pn+1(x) - (i+1)/(i+2) * Pn(x)% Synopsis: p = LegendreIter(n)% Input: n = 勒讓德多項式的次數(shù)% O

24、utput: p = n次勒讓德多項式的系數(shù)向量 if round(n) = n | n 0 error(n必須是一個非負(fù)整數(shù)); end if n = 0 %P0(x) = 1 p = 1; return; elseif n = 1 %P1(x) = x p = 1 0; return; end pBk = 1; %初始化三項遞推公式后項為P0 pMid = 1 0; %初始化三項遞推公式中項為P1 for i = 0:n-2 pMidCal = zeros(1,i+3); %構(gòu)造用于計算的x*Pn+1 pMidCal(1:i+2) = pMid; pBkCal = zeros(1,i+3)

25、; %構(gòu)造用于計算的Pn pBkCal(3:i+3) = pBk; pFwd = (2*i+3)/(i+2) * pMidCal - (i+1)/(i+2) * pBkCal; %勒讓德多項式三項遞推公式Pn+2(x) = (2*i+3)/(i+2) * x*Pn+1(x) - (i+1)/(i+2) * Pn(x) pBk = pMid; %把中項變?yōu)楹箜椷M行下次迭代 pMid = pFwd; %把前項變?yōu)橹许椷M行下次迭代 end p = pFwd;function p = ChebyshevIter(n)% ChebyshevIter 用遞推的方法計算n次勒讓德-切比雪夫多項式的系數(shù)向量

26、Tn+2(x) = 2*x*Tn+1(x) - Tn(x)% Synopsis: p = ChebyshevIter(n)% Input: n = 勒讓德-切比雪夫多項式的次數(shù)% Output: p = n次勒讓德-切比雪夫多項式的系數(shù)向量 if round(n) = n | n 0 error(n必須是一個非負(fù)整數(shù)); end if n = 0 %T0(x) = 1 p = 1; return; elseif n = 1 %T1(x) = x p = 1 0; return; end pBk = 1; %初始化三項遞推公式后項為T0 pMid = 1 0; %初始化三項遞推公式中項為T1 f

27、or i = 0:n-2 pMidCal = zeros(1,i+3); %構(gòu)造用于計算的x*Tn+1 pMidCal(1:i+2) = pMid; pBkCal = zeros(1,i+3); %構(gòu)造用于計算的Pn pBkCal(3:i+3) = pBk; pFwd = 2*pMidCal - pBkCal; %勒讓德-切比雪夫多項式三項遞推公式Tn+2(x) = 2*x*Tn+1(x) - Tn(x) pBk = pMid; %把中項變?yōu)楹箜椷M行下次迭代 pMid = pFwd; %把前項變?yōu)橹许椷M行下次迭代 end p = pFwd;function p = LaguerreIter(n

28、)% LaguerreIter 用遞推的方法計算n次拉蓋爾多項式的系數(shù)向量 Ln+2(x) = (2*n+3-x)*Ln+1(x) - (n+1)*Ln(x)% Synopsis: p = LaguerreIter(n)% Input: n = 拉蓋爾多項式的次數(shù)% Output: p = n次拉蓋爾多項式的系數(shù)向量 if round(n) = n | n 0 error(n必須是一個非負(fù)整數(shù)); end if n = 0 %L0(x) = 1 p = 1; return; elseif n = 1 %L1(x) = -x+1 p = -1 1; return; end pBk = 1; %初

29、始化三項遞推公式后項為L0 pMid = -1 1; %初始化三項遞推公式中項為L1 for i = 0:n-2 pMidCal1 = zeros(1,i+3); %構(gòu)造用于計算的x*Ln+1(x) pMidCal1(1:i+2) = pMid; pMidCal2 = zeros(1,i+3); %構(gòu)造用于計算的Ln+1(x) pMidCal2(2:i+3) = pMid; pBkCal = zeros(1,i+3); %構(gòu)造用于計算的Ln(x) pBkCal(3:i+3) = pBk; pFwd =( (2*i+3)*pMidCal2 - pMidCal1 - (i+1)*pBkCal )/ (i+2); %拉蓋爾多項式三項遞推公式Ln+2(x) = (2*n+3-x)*Ln+1(x) - (n+1)2*Ln(x) pBk = pMid; %把中項變?yōu)楹箜椷M行下次迭代 pMid = pFwd; %把前項變?yōu)橹许椷M行下次迭代 end p = pFwd;function p = HermiteIter(n)