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

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

2、下總共編寫了六個算法程序，部分代碼參考文獻(xiàn)10-13,為了體現(xiàn)程序的正確性，以下程序都以為例進(jìn)行運算。原積分的精確值為牛頓-科特斯求積公式的matlab實現(xiàn)先用m文件定義一個名為f1.m的函數(shù)：% i是要調(diào)用第幾個被積函數(shù)g(i)，x是自變量function f=f1(i,x) g(1)=sqrt(x);if x=0 g(2)=1;elseg(2)=sin(x)/x;endg(3)=4/(1+x2);f=g(i);程序一：function c,g=ncotes(a,b,n,m)% a，b分別為積分的上下限；% n是子區(qū)間的個數(shù)；% m是調(diào)用上面第幾個被積函數(shù)；% 當(dāng)n=1時計算梯形公式；當(dāng)n=

3、2時計算辛浦生公式，以此類推； i=n; h=(b-a)/i; z=0;for j=0:i x(j+1)=a+j*h; s=1; if j=0 s=s; elsefor k=1:j s=s*k;endendr=1;if i-j=0 r=r;elsefor k=1:(i-j) r=r*k;endendif mod(i-j),2)=1 q=-(i*s*r);else q=i*s*r;endy=1;for k=0:i if k=j y=y*(sym(t)-k); endendl=int(y,0,i);c(j+1)=l/q; z=z+c(j+1)*f1(m,x(j+1);endg=(b-a)*z1）當(dāng)

4、輸入，時，即在matlab命令窗口輸入 ncotes(0,1,1,2)即可得用梯形公式的積分值和相應(yīng)科特斯系數(shù) 如圖3.12）當(dāng)輸入，時，即在matlab命令窗口輸入 ncotes(0,1,2,2)即可得用辛浦生公式的積分值和相應(yīng)科特斯系數(shù)如圖3.23）當(dāng)輸入，時，即在matlab命令窗口輸入 ncotes(0,1,4,2)即可得用科特斯公式的積分值和相應(yīng)科特斯系數(shù)如圖3.3圖 3.1 圖 3.2 圖3.3 復(fù)化求積公式的matlab實現(xiàn)一、復(fù)化梯形求積公式的matlab實現(xiàn)通過的個等步長節(jié)點逼近積分其中，。程序二：function s=trapr1(f,a,b,n)% f是被積函數(shù)；% a,

5、b分別為積分的上下限；% n是子區(qū)間的個數(shù)；% s是梯形總面積；h=(b-a)/n;s=0;for k=1:(n-1) x=a+h*k; s=s+feval(f,x);endformat long s=h*(feval(f,a)+feval(f,b)/2+h*s;先用m文件定義一個名為f.m的函數(shù)：function y=f(x)if x=0 y=1;else y=sin(x)/x;end在matlab命令窗口中輸入 trapr1(f,0,1,4) 回車得到 如圖3.4 圖 3.4若取子區(qū)間的個數(shù)在matlab命令窗口中輸入 trapr1(f,0,1,8) 回車得到 如圖3.5 圖3.5 龍貝格

6、求積公式的matlab實現(xiàn)構(gòu)造數(shù)表來逼近積分其中。表示數(shù)表的最后一行，最后一列的值。程序五：function r,quad,err,h=romber(f,a,b,n,delta)% f是被積函數(shù)% a,b分別是積分的上下限% n+1是t數(shù)表的列數(shù)% delta是允許誤差% r是t數(shù)表% quad是所求積分值m=1;h=b-a;err=1j=0;r=zeros(4,4);r(1,1)=h*(feval(f,a)+feval(f,b)/2while (errdelta)&(jn)|(j romber(f,0,1,5,0.5*(10(-8)回車得到 如圖3.6 圖3.6 高斯-勒讓德求積公式的mat

7、lab實現(xiàn)程序六：function a,x=guass1(n) i=n+1; f=(sym(t)2-1)i; f=diff(f,i); t=solve(f); for j=1:i for k=1:i x(j,k)=t(k)(j-1); end if mod(j,2)=0 b(j)=0; else b(j)=2/j; end end x=inv(x); for j=1:i a(j)=0; x(j)=0; for k=1:i a(j)=a(j)+x(j,k)*b(k); x(j)=x(j)+t(j); end x(j)=x(j)/k; endfunction g= guasslegendre (a,b,n,m)% a,b分別是積分的上下限；% n+1為節(jié)點個數(shù)；% m是調(diào)用f1.m中第幾個被積函數(shù)； a,x=guass1(n);g=0;for i=1:n+1 y(i)=(b-a)/2*x(i)+(a+b)/2; f(i)=f1(m,y(i); g=g+(b-a)/2*f(i)*a(i);end用m文件分別把上面兩個自定義函數(shù)定義為名為guass1.m函數(shù)和guasslegendre.m函數(shù)用m文件定義一個名為f1.m的函數(shù)function f=f1(i,x) g(1)=sqrt(x);if x=0 g(2)=1;e