MATLAB實(shí)驗(yàn)三 定積分的近似計(jì)算

1、實(shí)驗(yàn)三 定積分的近似計(jì)算一、問(wèn)題背景與實(shí)驗(yàn)?zāi)康睦门ｎD萊布尼茲公式雖然可以精確地計(jì)算定積分的值，但它僅適用于被積函數(shù)的原函數(shù)能用初等函數(shù)表達(dá)出來(lái)的情形如果這點(diǎn)辦不到或者不容易辦到，這就有必要考慮近似計(jì)算的方法在定積分的很多應(yīng)用問(wèn)題中，被積函數(shù)甚至沒(méi)有解析表達(dá)式，可能只是一條實(shí)驗(yàn)記錄曲線，或者是一組離散的采樣值，這時(shí)只能應(yīng)用近似方法去計(jì)算相應(yīng)的定積分本實(shí)驗(yàn)將主要研究定積分的三種近似計(jì)算算法：矩形法、梯形法、拋物線法對(duì)于定積分的近似數(shù)值計(jì)算，Matlab有專(zhuān)門(mén)函數(shù)可用二、相關(guān)函數(shù)（命令）及簡(jiǎn)介1 sum(a)：求數(shù)組a的和2 format long：長(zhǎng)格式，即屏幕顯示15位有效數(shù)字（注：由于本實(shí)驗(yàn)

2、要比較近似解法和精確求解間的誤差，需要更高的精度）3double()：若輸入的是字符則轉(zhuǎn)化為相應(yīng)的ASCII碼；若輸入的是整型數(shù)值則轉(zhuǎn)化為相應(yīng)的實(shí)型數(shù)值4 quad()：拋物線法求數(shù)值積分格式： quad(fun,a,b) ，注意此處的fun是函數(shù)，并且為數(shù)值形式的，所以使用*、/、等運(yùn)算時(shí)要在其前加上小數(shù)點(diǎn)，即 .*、./、.等例：Q = quad(1./(x.3-2*x-5),0,2);5 trapz()：梯形法求數(shù)值積分格式：trapz(x,y)其中x為帶有步長(zhǎng)的積分區(qū)間；y為數(shù)值形式的運(yùn)算（相當(dāng)于上面介紹的函數(shù)fun）例：計(jì)算x=0:pi/100:pi;y=sin(x);trapz(x

3、,y)6 dblquad()：拋物線法求二重?cái)?shù)值積分格式：dblquad(fun,xmin,xmax,ymin,ymax)，fun可以用inline定義，也可以通過(guò)某個(gè)函數(shù)文件的句柄傳遞例1：Q1 = dblquad(inline(y*sin(x), pi, 2*pi, 0, pi)順便計(jì)算下面的Q2，通過(guò)計(jì)算，比較Q1 與Q2結(jié)果（或加上手工驗(yàn)算），找出積分變量x、y的上下限的函數(shù)代入方法Q2 = dblquad(inline(y*sin(x), 0, pi, pi, 2*pi)例2：Q3 = dblquad(integrnd, pi, 2*pi, 0, pi)這時(shí)必須存在一個(gè)函數(shù)文件int

4、egrnd.m：function z = integrnd(x, y) z = y*sin(x);7fprintf（文件地址，格式，寫(xiě)入的變量）：把數(shù)據(jù)寫(xiě)入指定文件例：x = 0:.1:1;y = x; exp(x);fid = fopen(exp.txt,w); %打開(kāi)文件fprintf(fid,%6.2f %12.8fn,y); %寫(xiě)入fclose(fid) %關(guān)閉文件8 syms 變量1 變量2 ：定義變量為符號(hào)9sym(表達(dá)式)：將表達(dá)式定義為符號(hào)解釋?zhuān)篗atlab中的符號(hào)運(yùn)算事實(shí)上是借用了Maple的軟件包，所以當(dāng)在Matlab中要對(duì)符號(hào)進(jìn)行運(yùn)算時(shí)，必須先把要用到的變量定義為符號(hào)10

5、 int(f,v,a,b)：求f關(guān)于v積分，積分區(qū)間由a到b11 subs(f，x，a)：將 a 的值賦給符號(hào)表達(dá)式 f 中的 x，并計(jì)算出值若簡(jiǎn)單地使用subs(f)，則將f的所有符號(hào)變量用可能的數(shù)值代入，并計(jì)算出值三、實(shí)驗(yàn)內(nèi)容1 矩形法根據(jù)定積分的定義，每一個(gè)積分和都可以看作是定積分的一個(gè)近似值，即在幾何意義上，這是用一系列小矩形面積近似小曲邊梯形的結(jié)果，所以把這個(gè)近似計(jì)算方法稱為矩形法不過(guò)，只有當(dāng)積分區(qū)間被分割得很細(xì)時(shí)，矩形法才有一定的精確度針對(duì)不同的取法，計(jì)算結(jié)果會(huì)有不同，我們以為例（?。?，（1） 左點(diǎn)法：對(duì)等分區(qū)間，在區(qū)間上取左端點(diǎn)，即取，0.78789399673078，理論值，此

6、時(shí)計(jì)算的相對(duì)誤差（2）右點(diǎn)法：同（1）中劃分區(qū)間，在區(qū)間上取右端點(diǎn)，即取，0.78289399673078，理論值，此時(shí)計(jì)算的相對(duì)誤差（3）中點(diǎn)法：同（1）中劃分區(qū)間，在區(qū)間上取中點(diǎn)，即取，0.78540024673078，理論值，此時(shí)計(jì)算的相對(duì)誤差如果在分割的每個(gè)小區(qū)間上采用一次或二次多項(xiàng)式來(lái)近似代替被積函數(shù)，那么可以期望得到比矩形法效果好得多的近似計(jì)算公式下面介紹的梯形法和拋物線法就是這一指導(dǎo)思想的產(chǎn)物2 梯形法等分區(qū)間，相應(yīng)函數(shù)值為（）曲線上相應(yīng)的點(diǎn)為（）將曲線的每一段弧用過(guò)點(diǎn)，的弦（線性函數(shù)）來(lái)代替，這使得每個(gè)上的曲邊梯形成為真正的梯形，其面積為，于是各個(gè)小梯形面積之和就是曲邊梯形面積

7、的近似值，即，稱此式為梯形公式仍用的近似計(jì)算為例，取，0.78539399673078，理論值，此時(shí)計(jì)算的相對(duì)誤差很顯然，這個(gè)誤差要比簡(jiǎn)單的矩形左點(diǎn)法和右點(diǎn)法的計(jì)算誤差小得多3 拋物線法由梯形法求近似值，當(dāng)為凹曲線時(shí)，它就偏??；當(dāng)為凸曲線時(shí)，它就偏大若每段改用與它凸性相接近的拋物線來(lái)近似時(shí)，就可減少上述缺點(diǎn)，這就是拋物線法將積分區(qū)間作等分，分點(diǎn)依次為，對(duì)應(yīng)函數(shù)值為（），曲線上相應(yīng)點(diǎn)為（）現(xiàn)把區(qū)間上的曲線段用通過(guò)三點(diǎn)，的拋物線來(lái)近似代替，然后求函數(shù)從到的定積分：由于，代入上式整理后得同樣也有將這個(gè)積分相加即得原來(lái)所要計(jì)算的定積分的近似值：，即這就是拋物線法公式，也稱為辛卜生（Simpson）公式

8、仍用的近似計(jì)算為例，取，=0.78539816339745，理論值，此時(shí)計(jì)算的相對(duì)誤差4. 直接應(yīng)用Matlab命令計(jì)算結(jié)果（1） 數(shù)值計(jì)算方法1：int(1/(1+x2),x,0,1) (符號(hào)求積分)方法2：quad(1./(1+x.2),0,1) （拋物線法求數(shù)值積分）方法3：x=0:0.001:1; y=1./(1+x.2);trapz(x,y) （梯形法求數(shù)值積分）（2）數(shù)值計(jì)算方法1：int(int(x+y2,y,-1,1),x,0,2) （符號(hào)求積分）方法2：dblquad(inline(x+y2),0,2,-1,1) （拋物線法二重?cái)?shù)值積分）四、自己動(dòng)手1 實(shí)現(xiàn)實(shí)驗(yàn)內(nèi)容中的例子，

9、即分別采用矩形法、梯形法、拋物線法計(jì)算，取，并比較三種方法的精確程度2 分別用梯形法與拋物線法，計(jì)算，取并嘗試直接使用函數(shù)trapz()、quad()進(jìn)行計(jì)算求解，比較結(jié)果的差異3 試計(jì)算定積分（注意：可以運(yùn)用trapz()、quad()或附錄程序求解嗎？為什么？）4 將的近似計(jì)算結(jié)果與Matlab中各命令的計(jì)算結(jié)果相比較，試猜測(cè)Matlab中的數(shù)值積分命令最可能采用了哪一種近似計(jì)算方法？并找出其他例子支持你的觀點(diǎn)5 通過(guò)整個(gè)實(shí)驗(yàn)內(nèi)容及練習(xí)，你能否作出一些理論上的小結(jié)，即針對(duì)什么類(lèi)型的函數(shù)（具有某種單調(diào)特性或凹凸特性），用某種近似計(jì)算方法所得結(jié)果更接近于實(shí)際值？6 學(xué)習(xí)fulu2sum.m的程

10、序設(shè)計(jì)方法，嘗試用函數(shù) sum 改寫(xiě)附錄1和附錄3的程序，避免for 循環(huán)五、附錄附錄1：矩形法（左點(diǎn)法、右點(diǎn)法、中點(diǎn)法）（fulu1.m）format longn=100;a=0;b=1;inum1=0;inum2=0;inum3=0;syms x fxfx=1/(1+x2);for i=1:nxj=a+(i-1)*(b-a)/n; %左點(diǎn)xi=a+i*(b-a)/n; %右點(diǎn)fxj=subs(fx,x,xj); %左點(diǎn)值fxi=subs(fx,x,xi); %右點(diǎn)值fxij=subs(fx,x,(xi+xj)/2); %中點(diǎn)值 inum1=inum1+fxj*(b-a)/n; inum2=

11、inum2+fxi*(b-a)/n; inum3=inum3+fxij*(b-a)/n;endinum1inum2inum3integrate=int(fx,0,1)integrate=double(integrate)fprintf(The relative error between inum1 and real-value is about: %dnn,. abs(inum1-integrate)/integrate)fprintf(The relative error between inum2 and real-value is about: %dnn,.abs(inum2-inte

12、grate)/integrate) fprintf(The relative error between inum3 and real-value is about: %dnn,. abs(inum3-integrate)/integrate)附錄2：梯形法（fulu2.m）format longn=100;a=0;b=1;inum=0;syms x fxfx=1/(1+x2);for i=1:n xj=a+(i-1)*(b-a)/n; xi=a+i*(b-a)/n; fxj=subs(fx,x,xj); fxi=subs(fx,x,xi); inum=inum+(fxj+fxi)*(b-a)

13、/(2*n);endinumintegrate=int(fx,0,1)integrate=double(integrate)fprintf(The relative error between inum and real-value is about: %dnn,.abs(inum-integrate)/integrate)附錄2sum：梯形法（fulu2sum.m），利用求和函數(shù)，避免for 循環(huán)format longn=100;a=0;b=1;syms x fxfx=1/(1+x2);i=1:n;xj=a+(i-1)*(b-a)/n; %所有左點(diǎn)的數(shù)組xi=a+i*(b-a)/n; %所有

14、右點(diǎn)的數(shù)組fxj=subs(fx,x,xj); %所有左點(diǎn)值fxi=subs(fx,x,xi); %所有右點(diǎn)值f=(fxi+fxj)/2*(b-a)/n; %梯形面積inum=sum(f) %加和梯形面積求解integrate=int(fx,0,1)integrate=double(integrate)fprintf(The relative error between inum and real-value is about: %dnn,. abs(inum-integrate)/integrate)附錄3：拋物線法（fulu3.m）format longn=100;a=0;b=1;inum=0;syms x fxfx=1/(1+x2);for i=1:nxj=a+(i-1)*(b-a)/n; %左點(diǎn)xi=a+i*(b-a)/n; %右點(diǎn)xk=(xi+xj)/2; %中點(diǎn)fxj=subs(fx,x,xj); fxi=su