MATLAB實驗三 定積分的近似計算

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

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

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

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

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

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

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

8、仍用的近似計算為例，取，=0.78539816339745，理論值，此時計算的相對誤差4. 直接應(yīng)用Matlab命令計算結(jié)果（1） 數(shù)值計算方法1：int(1/(1+x2),x,0,1) (符號求積分)方法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ù)值計算方法1：int(int(x+y2,y,-1,1),x,0,2) （符號求積分）方法2：dblquad(inline(x+y2),0,2,-1,1) （拋物線法二重數(shù)值積分）四、自己動手1 實現(xiàn)實驗內(nèi)容中的例子，

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

10、序設(shè)計方法，嘗試用函數(shù) sum 改寫附錄1和附錄3的程序，避免for 循環(huán)五、附錄附錄1：矩形法（左點法、右點法、中點法）（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; %左點xi=a+i*(b-a)/n; %右點fxj=subs(fx,x,xj); %左點值fxi=subs(fx,x,xi); %右點值fxij=subs(fx,x,(xi+xj)/2); %中點值 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; %所有左點的數(shù)組xi=a+i*(b-a)/n; %所有

14、右點的數(shù)組fxj=subs(fx,x,xj); %所有左點值fxi=subs(fx,x,xi); %所有右點值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; %左點xi=a+i*(b-a)/n; %右點xk=(xi+xj)/2; %中點fxj=subs(fx,x,xj); fxi=su