MATLAB數(shù)值分析實(shí)驗(yàn)二(復(fù)合梯形、辛普森和龍貝格求積-以及二重積分計(jì)算等)(共5頁)

1、佛山科學(xué)技術(shù)學(xué)院實(shí) 驗(yàn) 報(bào) 告課程名稱 數(shù)值分析 實(shí)驗(yàn)項(xiàng)目 數(shù)值積分 專業(yè)班級(jí) 機(jī)械工程 姓 名 余紅杰 學(xué) 號(hào) 2111505010 指導(dǎo)教師 陳劍 成 績 日 期 月 日 一、實(shí)驗(yàn)?zāi)康?、理解如何在計(jì)算機(jī)上使用數(shù)值方法計(jì)算定積分的近似值;2、學(xué)會(huì)復(fù)合梯形、復(fù)合Simpson和龍貝格求積分公式的編程與應(yīng)用。3、探索二重積分在矩形區(qū)域的數(shù)值積分方法。二、實(shí)驗(yàn)要求（1） 按照題目要求完成實(shí)驗(yàn)內(nèi)容；（2） 寫出相應(yīng)的Matlab 程序；（3） 給出實(shí)驗(yàn)結(jié)果(可以用表格展示實(shí)驗(yàn)結(jié)果)；（4） 分析和討論實(shí)驗(yàn)結(jié)果并提出可能的優(yōu)化實(shí)驗(yàn)。（5） 寫出實(shí)驗(yàn)報(bào)告。三、實(shí)驗(yàn)步驟1、用不同數(shù)值方法計(jì)算積分（1）取

2、不同的步長，分別用復(fù)合梯形及復(fù)合辛普森求積計(jì)算積分，給出誤差中關(guān)于的函數(shù)，并與積分精確值比較兩公式的精度。（2）用龍貝格求積計(jì)算完成問題（1）。2、給出一種求矩形區(qū)域上二重積分的復(fù)化求積方法，然后計(jì)算二重積分，其中積分區(qū)域。1.%Int_t.m 復(fù)化梯形：function F = Int_t(x1,x2,n)% 復(fù)化梯形求積公式% x1,x2 為積分起點(diǎn)和中點(diǎn)%分為n個(gè)區(qū)間，沒選用步長可以防止區(qū)間數(shù)為非整數(shù)。%樣點(diǎn)矩陣及其函數(shù)值：x = linspace(x1,x2,n+1);y = f(x);m = length(x);%本題中用Matlab計(jì)算端點(diǎn)位置函數(shù)值為NaN,故化為零：y(1) =

3、 0;y(m) = 0;%算出區(qū)間長度，步長h：h = (x2 -x1)/n;a = 1 2*ones(1,m-2) 1;%計(jì)算估計(jì)的積分值：F = h/2*sum(a.*y);%f.mfunction y = f(x)y = sqrt(x).*log(x);%run11.mclc,clear;%分為10個(gè)區(qū)間，步長0.1的積分值：F = Int_t(0,1,10);F10 = F%分為100個(gè)區(qū)間F = Int_t(0,1,100);F100 = F%誤差計(jì)算W10 = abs(-4/9)-F10);W100 = abs(-4/9)-F100);W = W10 W100%復(fù)化辛普森：%Int

4、_s.mfunction F = Int_s(x1,x2,n)% 復(fù)化梯形求積公式% x1,x2 區(qū)間，分為n個(gè)區(qū)間。%樣點(diǎn)矩陣及其函數(shù)值：x = linspace(x1,x2,n+1);y = f(x);m = length(x);h = (x2 -x1)/n;y(1)=0;y(m)=0;%本題中用Matlab計(jì)算端點(diǎn)位置函數(shù)值為NaN,故化為零：F1=sum(y);xo = x + h/2;xo(m) = ;y = f(xo);F2 = sum(y);F = (h/6)*(2*F1 + 4*F2);%run112.mclc,clear;%分為10個(gè)區(qū)間，步長0.1的積分值：F = Int_

5、s(0,1,10);S10 = F%分為100個(gè)區(qū)間F = Int_s(0,1,100);S100 = F%誤差計(jì)算W10 = abs(-4/9)-S10);可以明顯看出其精度高于復(fù)化梯形W100 = abs(-4/9)-S100);W = W10 W100%run113.m 擬合誤差和步長之間的三次曲線關(guān)系。clc,clear;%建立梯形誤差、辛普森誤差、步長矩陣：T=zeros(1,10);S=zeros(1,10);h=zeros(1,10);for i=1:10 F = Int_t(0,1,10*i); T(i) = -4/9 - F; F = Int_s(0,1,10*i); S(i

6、) = -4/9 - F; h(i) = 1/(10*i);end也可以明顯看出辛普森誤差曲線各項(xiàng)系數(shù)都較小的TP = polyfit(h,T,3)SP = polyfit(h,S,3)%龍貝格：%Romberg.m：function F=Romberg(x1,x2,n)%建立龍貝格推算矩陣、求最初步長：R = zeros(4);h = (x2-x1)/n;x = linspace(x1,x2,n+1);%計(jì)算矩陣第一列：復(fù)化梯形結(jié)果：for i=1:4 F = Int_t(x1,x2,n*i); R(i,1) = F;end%計(jì)算第二列：辛普森for i =1:3 R(i,2) = (4/3

7、)*R(i+1,1)-(1/3)*R(i,1);end%計(jì)算第三列：復(fù)化斯科特for i =1:2 R(i,3) = (16/15)*R(i+1,2)-(1/15)*R(i,2);endR(1,4) = (64/63)*R(2,3)-(1/63)*R(1,3);F = R(1,4);%run12.m clc,clear;F = Romberg(0,1,10);F10 = F;F = Romberg(0,1,100);F100 = F;F10,F100;-4/9-F10,-4/9-F100%右圖為初始劃分為10個(gè)區(qū)間和100個(gè)區(qū)間進(jìn)行運(yùn)算的結(jié)果，可以看出初始10次劃分的精度比辛普森和梯形結(jié)果高出

8、不少2，我選用的是類似于復(fù)化梯形的求積方法，由于是正方形的積分區(qū)域，把它劃分為n*n塊，在其中1.（x(k),y(k)）,2.(x(k+1),y(k),3.(x(k),y(k+1),4.(x(k+1),y(k+1)的小區(qū)間上ds=dxdy是已知的，而且存在一e個(gè)f(e)ds=該函數(shù)在此小區(qū)域的數(shù)值，此處就選用（f(1),f(2),f(3),f(4)）/4的值來近似，可知隨著n的增大，是會(huì)越來越逼近真實(shí)值的。程序如下：%qes2.m 被積分函數(shù)：function f=qes2(x,y)f = exp(-x*y);%Int_xy.m 積分代碼：function F = Int_xy(x1,x2,y

9、1,y2,n)%僅適合此題目的正方形區(qū)域二重積分%x1,x2,y1,y2分別為橫坐標(biāo)和縱坐標(biāo)起點(diǎn)與中點(diǎn)；%n為劃分為n*n的小正方形網(wǎng)格%建立所有網(wǎng)格交點(diǎn)處的橫坐標(biāo)和縱坐標(biāo)：x = linspace(x1,x2,n+1);y = linspace(y1,y2,n+1);%計(jì)算每個(gè)小方格的面積，dxdy：ds = (x2-x1)*(y2-y1)/(n2);%建立一個(gè)矩陣存放各個(gè)小方塊的積分值：Z = zeros(n);%初始化積分值：F = 0;%把小方塊的積分值存入矩陣Z中，第一行第一列存放最左下角方格數(shù)據(jù)，以此類推：%并依次把小方格區(qū)域的積分值累加，得到最終積分值。for i = 1:n for j = 1:n f = qes2(x(i),y(j); f1= f; f = qes2(x(i+1),y(j); f2= f; f = qes2(x(i),y(j+1); f3= f; f = qes2(x(i+1),y(j+1); f4= f; Z(i,j) =(ds * (f1+f2+f3+f4)/4; F = F + Z(i,j); endend%run2.m 運(yùn)行代碼：clc,clear;