數(shù)值積分及MATLAB實(shí)現(xiàn)綜述

1、數(shù)值積分及MATLAB實(shí)現(xiàn)綜述各種求積公式的MATLAB編程實(shí)現(xiàn)與應(yīng)用MATLAB是山MathWorks公式開發(fā)的一種主要用于數(shù)值計(jì)算及可視化圖 形處理的工程語言，是當(dāng)今最優(yōu)秀的科技應(yīng)用軟件之一。它將數(shù)值計(jì)算、矩陣 運(yùn)算、圖形圖像處理、信號(hào)處理和仿真等諸多強(qiáng)大的功能集成在較易使用的交 互計(jì)算機(jī)環(huán)境中，為科學(xué)研究、工程應(yīng)用提供了一種功能強(qiáng)、效率高的編程工 具。下面我們將各種求積算法通過MATLAB軟件編程實(shí)現(xiàn)，以下程序均用 MATLAB7.0 編寫，運(yùn)行壞境：1、硬件環(huán)境 CPU (intel Core i3-2310M,2.1GHz), 內(nèi)存(2GB昱聯(lián))，2、軟件環(huán)境windows7 (32

2、位)操作系統(tǒng)。以下總共編寫 了六個(gè)算法程序，部分代碼參考文獻(xiàn)10-13,為了體現(xiàn)程序的正確性，以下程 序都以為例進(jìn)行運(yùn)算。原積分的精確值為Ju A1 = 1= 0.946083070367183Jo x牛頓-科特斯求積公式的MATLAB實(shí)現(xiàn)先用M文件定義一個(gè)名為fl.m的函數(shù)：% i是要調(diào)用第兒個(gè)被積函數(shù)g(i), x是自變量function f=fl(i,x) g(l)=sqrt(x);if x=0g =1;elseg(2)=sin(x)/x;endg(3)=4/(l+xA2);f=g ；程序一：function C,g=NCotes(a,b5njn)% a, b分別為積分的上下限；% n是

3、子區(qū)間的個(gè)數(shù)；% m是調(diào)用上面第兒個(gè)被積函數(shù)；%當(dāng)n=l時(shí)計(jì)算梯形公式；當(dāng)n=2時(shí)計(jì)算辛浦生公式，以此類推; i=n;h=(b-a)/i;z=0;for j=0:ix(j+l)=a+j*h;s=l;if j=0s=s;elsefor k=l:js=s*k;endendr=l；if i-j=Or=r;elsefor k=l:(i-j)r=r*k;endif mod(i-j),2)=lq=-(i*s*r);elseq=i*s*r;endfor k=0:iif k=j y=y*(sym(t)-k)； endl=int(y,OJ)；C(j+l)=l/q；z=z+C(j+1 )*fl (m,x(j+1)

4、;endg=(b-a)*z1) 當(dāng)輸入a = O.b = , n = l9m = 2時(shí),即在MATLAB命令窗口輸入 NCotes(0,l.l,2)即可得用梯形公式的積分值和相應(yīng)科特斯系數(shù)如圖3.12) 當(dāng)輸入4 = 0上=1, n = 2,加=2時(shí)，即在MATLAB命令窗口輸入 NCotes(0,l,2,2)即可得用辛浦生公式的積分值和相應(yīng)科特斯系數(shù)如圖3.23) 當(dāng)輸入。=0上=1, n = 49m = 2B寸,即在MATLAB命令窗口輸入 NCotes(0JA2)即可得用科特斯公式的積分值和相應(yīng)科特斯系數(shù)如圖3.3圖3.1圖3.2圖3.3復(fù)化求積公式的MATLAB實(shí)現(xiàn)、復(fù)化梯形求積公式的

5、MATLAB實(shí)現(xiàn)通過/(X)的H + 1個(gè)等步長節(jié)點(diǎn)逼近積分f f(x)dx (/(a) + 于)+ hf(xk)a2A=1其中，Xk=a + kh , x0 =a.xn=h o程序二function s=traprl(f,a,bji)% f是被積函數(shù)；% a,b分別為積分的上下限；% n是子區(qū)間的個(gè)數(shù)；% s是梯形總面積；h=(b-a)/n;s=0;for k=l:(n-l) x=a+h*k; s=s+fevaI(T,x)；end format long s=h*(feval(,f,a)+feval(,f,b)/2+h*s;先用M文件定義一個(gè)名為f.m的函數(shù):function y=f(x)i

6、f x=0 y=i;else y=sin(x)/x; end在MATLAB命令窗口中輸入traprl(f,0丄4)回車得到如圖3.4圖3.4若取子區(qū)間的個(gè)數(shù)77 = 8在MATLAB命令窗口中輸入 traprl(TAl,8)回車得到如圖3.5圖3.5龍貝格求積公式的MATLAB實(shí)現(xiàn)構(gòu)造T數(shù)表來逼近積分 fWdx =其中。R(JJ)表示了數(shù)表的最后一行，最后一列的值。程序五：function R.quad,errji=romber(f,a.b,n Jelta)% f是被積函數(shù)% a,b分別是積分的上下限% n+1是T數(shù)表的列數(shù)% delta是允許誤差% R是T數(shù)表% quad是所求積分值M=l;

7、h=b-a;err=lJ=0;R=zeros(4,4);R( 1J )=h*(feval(Y ,a)+feval(T ,b)/2 while (errdelta)&(Jn)l(J4)J=J+1;h 二 h/2;s=0;for p=l:Mx=a+h*(2*p-l); s=s+feval(T,x);endR(J+lJ)=R(J,l)/2+h*s;M=2*M;for K=1:JR(J+ 1,K+1)=R(J+1 ,K)+(R(J+1 ,K)-R(J,K)/(4AK-1); enderr=abs(R(JJ)-R(J+1 ,K+1);endquad 二 R(J+1 J+l)先用M文件定義一個(gè)名為f.m的

8、函數(shù)：function y=f(x) if x=0else y=sin(x)/x;end在MATLAB命令窗口中輸入 romberCf,0,h5,0.5*(10A(-8)回車得到如圖3.6Convnand Window ronber( f * , 0, 5, 0. 5*(10 (-8)err =1R =0.9207354920395000000000000000quad =CL 9460B307036718ans =0920254924030600000939703284806180. 9461458822735000094451352166539CL 9460BS93395179CL 946

9、03300406367000945690863582700 9460B3310888470 9460B3069350920 94605307039722009458502993430. 946083085384950. 946083070361380. 946083070367260. 94608307036718l3圖3. 6高斯-勒讓德求積公式的MATLAB實(shí)現(xiàn)程序六：function A,x二Guassl(N)i=N+l;f=(sym(t,)A2-l)Ai;f=diff(f,i);t=solve(f);for j=l:ifor k=l:iX(j,k)=t(k)l);endif mod(j

10、,2)=0B(j)=O;elseB(j)=2/j;endendX=inv(X);for j=l:iA(j)=0； x(j)=O； for k=l:iA(j)=A(j)+X(j,k)*B(k); x(j)=x(j)+t(j)； endx(j)=x(j)/k;endfunction g= GuassLegendre (a,b,n,m)% a,b分別是積分的上下限；% n+1為節(jié)點(diǎn)個(gè)數(shù)；% m是調(diào)用fl.m中第兒個(gè)被積函數(shù)：A,x二Guassl(n);g=o；for i=l:n+ly(i)=(b-a)/2*x(i)+(a+b)/2;f(i)=fl(m,y(i);g=g+(b-a)/2*f(i)*A(i);end用M文件分別把上面兩個(gè)自定義函數(shù)定義為名為Guassl.m函數(shù)和 GuassLegendre.m 函數(shù)用M文件定義一個(gè)名為fl.