【導讀】數(shù)值積分就是解決此類問題的一種行之有效的方法。析的一個重要分支；因此，探討近似計算的數(shù)值積分方法是有著明顯的實際意義的。數(shù)值積分問題的產(chǎn)生出發(fā)，詳細介紹了一些數(shù)值積分的重要方法。分公式，即龍貝格求積公式和高斯-勒讓德求積公式。除了研究這些數(shù)值積分算法的理論外，積公式進行運算，分析比較了各種求積公式的計算誤差。

　　

【正文】 1 ( ) ( ) ( ) ( ) ( )f x dx A f x A f x A f x A f x? ? ? ? ??其代數(shù)精確度為 7 綜上，高斯 勒讓德求積公式的積分節(jié)點和積分系數(shù)表如 表 表 積分節(jié)點和積分系數(shù)表 重慶郵電大學本科 畢業(yè)設計（論文） 22 三 、 高斯 勒讓德求積公式的余項 對 n 階 Gauss Legendre? 求積公式（共有 1n? 個求積節(jié)點） 1( 2 2 ) 2111 ( ) ( )( 2 2) ! nnnR f x dxn ??? ??? ? ? () 其中， ( 1,1)??? ， 1 0 1( ) ( ) ( ) ( )nnx x x x x x x? ? ? ? ? ? 01, , , nx x x 為多項式 1()npx? 的零點 。 積分，可以計算得到： 2 3 4 ( 2 2 )32 [ ( 1 ) ! ] ()[ ( 2 2 ) ! ] ( 2 3 )n nn nRfnn ?? ????? () 四 、 一般區(qū)間上的高斯 勒讓德求積公式 前面講解的高斯 勒讓德 求積公式是計算標準區(qū)間 [1,1]? 上的定積分，對于一般的有界區(qū)間 [, ]ab 上的定積分 ，可以通過變量代換轉(zhuǎn)化為區(qū)間 [1,1]? 的定積分。 即在積分 ()ba f x dx?中，令 22b a a bxt????，當 xa? 時， 1t?? ；當 xb? 時， 1t? ；且 2badx dt?? 。這個變換為線性變換，其反變換為 2 batxb a b a????? 于是成立 11( ) ( )2 2 2bab a b a a bf x d x f t d t?? ? ????? 然后，我們可以用 [1,1]? 區(qū)間上的高斯 勒讓德求積公式 。 第 五 節(jié) 復化兩點高斯 勒讓德 求積 公式 本 節(jié) 對兩點 Gauss Legendre? 公式 11 33( ) ( ) ( )f x d x f f? ? ? ?? () 進行推廣和復化，得到了新的數(shù)值積分公式 重慶郵電大學本科 畢業(yè)設計（論文） 23 11111 0 2233( ) ( ( ) ( ) )2 6 6niiih h hf x d x f x f x?? ???? ? ? ??? () 同時分析了它的積分誤差和收斂階 定義 [8] 設 ()baI f x dx??， nI 是一種復化求積公式， 若 當 0h? 時成立 0lim , 0nhhII CCh?? ?? () 則稱求積公式 nI 是 p 階收斂的。 例如，復化的 Trapezoid 公式和復化的 Simpson 公式分別具有二階和四階收斂 性。 定理 [8] 設 (4)( ) [ , ]f x C a b? ， bah n?? ，則復化兩點高斯 勒讓德求積公式為 1110 2233( ) ( ( ) ( ) )2 6 6nba iikh h hf x d x f x f x????? ? ? ??? ) 的余項表達式為 4 ( 4 )5()[ ] ( ) , [ , ]2 135b a hR f f a b?????? () 該方法具有四階收斂性 。 考慮積分 210 xI e dx???， 準確解為 ? 用復化兩點高斯 勒讓德公式可求得 I = ，其絕對誤差為R =，與復化梯形公式，復化 Simpson公式以及原兩點 GaussLegendre 公式相比， 這里 構(gòu)造的 復化兩點高斯 勒讓德求積公式 方法在運算量和精度方面都有顯著的 改善 。 第 六 節(jié) 本章小結(jié) 本章主要介紹了精確度比較高的兩種數(shù)值求積公式，即龍貝格求積公式和高斯 型 求積公式 ， 對其進行了相應的誤差分析。 其中高斯型求積公式主要介紹了高斯 勒讓德求積公式， 并對兩點高斯 勒讓德求積公式進行了復化。 重慶郵電大學本科 畢業(yè)設計（論文） 24 第 三 章 各種求積公式的 MATLAB 編程實現(xiàn)與應用 MATLAB 是由 MathWorks 公式開發(fā)的一種主要用于數(shù)值計算及可視化圖形處理的工程語言，是當今最優(yōu)秀的科技應用軟件之一。它將數(shù)值計算、矩陣運算、圖形圖像處理、信號處理和仿真等諸多 強大 的功能集成在較易使用的交互計算機環(huán)境中， 為科學研究、工程應用提供了一種功能強、效率高的編程工具。 下面我們將各種求積算法通過 MATLAB 軟件 編程實現(xiàn)， 以下程序 均 用 編寫，運行壞境 ： 硬件環(huán)境 CPU（ intel Core i32310M,） ，內(nèi)存（ 2GB 昱聯(lián)） ,軟件環(huán)境 windows7（ 32 位） 操作系統(tǒng)。 以下總 共 編寫了六個算法程序，部分 代碼 參考文獻 [1013],為了體現(xiàn)程序的正確性， 以下程序 都以 10 sinxI dxx??為例 進行運算 。 原積分的精確值為 10s in 0 . 9 4 6 0 8 3 0 7 0 3 6 7 1 8 3xI d xx??? 第一節(jié) 幾個低次牛頓 科特斯求積公式 的 MATLAB 實現(xiàn) 先用 M 文件定義一個名為 ： % i 是要調(diào)用第幾個 被積 函數(shù) g(i)， x 是自變量 function f=f1(i,x) g(1)=sqrt(x)。 if x==0 g(2)=1。 else g(2)=sin(x)/x。 end g(3)=4/(1+x^2)。 重慶郵電大學本科 畢業(yè)設計（論文） 25 f=g(i)。 程序 一 ： function [C,g]=NCotes(a,b,n,m) % a， b 分別為積分的上下限； % n 是子區(qū)間的個數(shù)； % m是調(diào)用上面第幾個被積函數(shù)； % 當 n=1 時計算梯形公式；當 n=2 時計算辛浦生公式，以此類推； i=n。 h=(ba)/i。 z=0。 for j=0:i x(j+1)=a+j*h。 s=1。 if j==0 s=s。 else for k=1:j s=s*k。 end end r=1。 if ij==0 r=r。 else for k=1:(ij) r=r*k。 end end if mod((ij),2)==1 q=(i*s*r)。 else q=i*s*r。 end y=1。 for k=0:i if k~=j y=y*(sym(39。t39。)k)。 end end 重慶郵電大學本科 畢業(yè)設計（論文） 26 l=int(y,0,i)。 C(j+1)=l/q。 z=z+C(j+1)*f1(m,x(j+1))。 end g=(ba)*z 1） 當輸入 0, 1ab??， 1, 2nm??時，即在 MATLAB 命令窗口輸入 NCotes(0,1,1,2)即可得 用 梯形公式的積分值 和相應 科特斯系數(shù) 如圖 2） 當輸入 0, 1ab??， 2, 2nm??時，即在 MATLAB 命令窗口輸入 NCotes(0,1,2,2)即可得 用 辛浦生公式的積分值 和相應 科特斯系數(shù) 如圖 3） 當輸入 0, 1ab??， 4, 2nm??時，即在 MATLAB 命令窗口輸入 NCotes(0,1,4,2)即可得 用科特斯 公式的積分值 和相應 科特斯系數(shù) 如圖 圖 重慶郵電大學本科 畢業(yè)設計（論文） 27 圖 圖 重慶郵電大學本科 畢業(yè)設計（論文） 28 第二節(jié) 復化求積公式 的 MATLAB 實現(xiàn) 一 、 復化梯形 求積 公式的 MATLAB 實現(xiàn) 通過 ()fx的 1n? 個等步長節(jié)點逼近積分 11( ) ( ( ) ( ) ) ( )2nbka khf x d x f a f b h f x??? ? ? ?? 其中， kx a kh?? ， 0 , nx a x b??。 程序 二 ： function s=trapr1(f,a,b,n) % f 是被積函數(shù) ； % a,b 分別為積分的上下限 ； % n 是子區(qū)間的個數(shù) ； % s 是梯形總面積 ； h=(ba)/n。 s=0。 for k=1:(n1) x=a+h*k。 s=s+feval(39。f39。,x)。 end format long s=h*(feval(39。f39。,a)+feval(39。f39。,b))/2+h*s。 先用 M 文件定義一個名為 ： function y=f(x) if x==0 y=1。 else y=sin(x)/x。 end 在 MATLAB 命令窗口中輸入 trapr1(39。f39。,0,1,4) 回車得到 如圖 重慶郵電大學本科 畢業(yè)設計（論文） 29 圖 若取子區(qū)間的個數(shù) 8n? 在 MATLAB 命令窗口中輸入 trapr1(39。f39。,0,1,8) 回車得到 如圖 圖 二、 復化辛浦生 求積 公式的 MATLAB 實現(xiàn) 程序 三 ： function s=simpr1(f,a,b,n) 重慶郵電大學本科 畢業(yè)設計（論文） 30 % f 是被積函數(shù) ； % a,b 分別為積分的上下限 ； % n 是子區(qū)間的個數(shù) ； % s 是梯形總面積 ,即所求積分數(shù)值 ； h=(ba)/(2*n)。 s1=0。 s2=0。 for k=1:n x=a+h*(2*k1)。 s1=s1+feval(39。f39。,x)。 end for k=1:(n1) x=a+h*2*k。 s2=s2+feval(39。f39。,x)。 end s=h*(feval(39。f39。,a)+feval(39。f39。,b)+4*s1+2*s2)/3。 先用 M 文件定義一個名為 ： function y=f(x) if x==0 y=1。 else y=sin(x)/x。 end 在 MATLAB 命令窗口中輸入 simpr1(39。f39。,0,1,4) 回車得到 如圖 圖 重慶郵電大學本科 畢業(yè)設計（論文） 31 若取 子區(qū)間個數(shù) 8n? 時 在 MATLAB 命令窗口中輸入 simpr1(39。f39。,0,1,8) 回車得到 如圖 圖 三、 復化科特斯求積公式的 MATLAB 實現(xiàn) 程序四 ： function s=cotespr1(f,a,b,n) % f 是被積函數(shù) ； % a,b 分別為積分的上下限 ； % n 是子區(qū)間的個數(shù) ； % s 是梯形總面積 ,即所求積分數(shù)值 ； h=(ba)/n。 s1=0。 s2=0。 s3=0。 s4=0。 for k=1:n x=a+(4*k3)*h/4。 s1=s1+feval(39。f39。,x)。 end for k=1:n 重慶郵電大學本科 畢業(yè)設計（論文） 32 x=a+(4*k2)*h/4。 s2=s2+feval(39。f39。,x)。 end for k=1:n x=a+(4*k1)*h/4。 s3=s3+feval(39。f39。,x)。 end for k=1:(n1) x=a+4*k*h/4。 s4=s4+feval(39。f39。,x)。 end s=h*(7*feval(39。f39。,a)+7*feval(39。f39。,b)+32*s1+12*s2+32*s3+14*s4)/90。 同樣 先用 M 文件定義一個名為 ： function y=f(x