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

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

2、下總共編寫(xiě)了六個(gè)算法程序，部分代碼參考文獻(xiàn)10-13,為了體現(xiàn)程序的正確性，以下程序都以為例進(jìn)行運(yùn)算。原積分的精確值為牛頓-科特斯求積公式的matlab實(shí)現(xiàn)先用m文件定義一個(gè)名為f1.m的函數(shù)：% i是要調(diào)用第幾個(gè)被積函數(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ū)間的個(gè)數(shù)；% m是調(diào)用上面第幾個(gè)被積函數(shù)；% 當(dāng)n=1時(shí)計(jì)算梯形公式；當(dāng)n=

3、2時(shí)計(jì)算辛浦生公式，以此類(lèi)推； 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、輸入，時(shí)，即在matlab命令窗口輸入 ncotes(0,1,1,2)即可得用梯形公式的積分值和相應(yīng)科特斯系數(shù) 如圖3.12）當(dāng)輸入，時(shí)，即在matlab命令窗口輸入 ncotes(0,1,2,2)即可得用辛浦生公式的積分值和相應(yīng)科特斯系數(shù)如圖3.23）當(dāng)輸入，時(shí)，即在matlab命令窗口輸入 ncotes(0,1,4,2)即可得用科特斯公式的積分值和相應(yīng)科特斯系數(shù)如圖3.3圖 3.1 圖 3.2 圖3.3 復(fù)化求積公式的matlab實(shí)現(xiàn)一、復(fù)化梯形求積公式的matlab實(shí)現(xiàn)通過(guò)的個(gè)等步長(zhǎng)節(jié)點(diǎn)逼近積分其中，。程序二：function s=trapr1(f,a,b,n)% f是被積函數(shù)；% a,

5、b分別為積分的上下限；% n是子區(qū)間的個(gè)數(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文件定義一個(gè)名為f.m的函數(shù)：function y=f(x)if x=0 y=1;else y=sin(x)/x;end在matlab命令窗口中輸入 trapr1(f,0,1,4) 回車(chē)得到 如圖3.4 圖 3.4若取子區(qū)間的個(gè)數(shù)在matlab命令窗口中輸入 trapr1(f,0,1,8) 回車(chē)得到 如圖3.5 圖3.5 龍貝格

6、求積公式的matlab實(shí)現(xiàn)構(gòu)造數(shù)表來(lái)逼近積分其中。表示數(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)回車(chē)得到 如圖3.6 圖3.6 高斯-勒讓德求積公式的mat

7、lab實(shí)現(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é)點(diǎn)個(gè)數(shù)；% m是調(diào)用f1.m中第幾個(gè)被積函數(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文件分別把上面兩個(gè)自定義函數(shù)定義為名為guass1.m函數(shù)和guasslegendre.m函數(shù)用m文件定義一個(gè)名為f1.m的函數(shù)function f=f1(i,x) g(1)=sqrt(x);if x=0 g(2)=1;e