數(shù)值分析試驗報告Matlab實現(xiàn)

1、學生實驗報告實驗課程名稱數(shù)值分析開課實驗室數(shù)學與統(tǒng)計學院實驗室學 院 2010 年級數(shù)學與應(yīng)用數(shù)學專業(yè)班01班學生姓名 學 號開課時 間 2012 至 2013 學年第 一 學期總成績教師簽名根據(jù)高等代數(shù)的知識，若| a|。0,上式的解存在且唯一。(1)(2)課程名稱數(shù)值分析實驗項目名 稱Gauss Mte 法實驗項目類型驗證演示綜合設(shè)計其他指導教師何光輝成 績實驗目的：(1) 高斯列主元消去法求解線性方程組的過程(2) 熟悉用迭代法求解線性方程組的過程(3) 設(shè)計出相應(yīng)的算法，編制相應(yīng)的函數(shù)子程序?qū)嶒瀮?nèi)容分別用高斯列主元消元法和直接消元法求解線性方程組:一2100- 3 1 X1 I一 10

2、1_ 3-4-1213|X251 1123-4X31-2 1149T3 一乂 一1 17 一實驗原理對于線性方程組811X1+812X2+ +anXn=b821X1+822X2+ +a2nXn=b2an1X1 +2為 * +annXn=bn常記為矩陣形式Ax = b15(1) Gauss直接消元法考慮上述線性方程組的增廣矩陣 A：b,對增廣矩陣進行行變換，將(2)式化為等價的三角形方陣，然后回代解之，這就是Gauss消元法。具體如下:a)消元令a/=aj，1，j =1,2,，n ; W)= h , i =1,2,，n對k=1到n-1，若a；：)黃0 ,進行1 一a(kk)lik /, akka

3、(k書 aik=0,(k+) a 口 ja(k) = aij -4 >akjk)=b-ik * A ,i =k +1,k +2,，ni =k +1,k +2,，ni, j =k +1,k +2,，ni =k +1,k +2,，nb)回代，若a；n #0，b；Xn =Ebnn1(i) v (i)、Xi = i (bi - ' aj Xj)aiij 占 1(2) Gauss列主元消元法設(shè)列主元消元法已完成Ax=b的第k-1 (1玄kn-1)次消元，的到方程組Ax = b A(k)x = b(k)在進行第k次消元前，先進行 2個步驟：(k)(k)(k)(k)l(k)A =0a) 在ak

4、k至ank這一列內(nèi)選出取大值，即禮*| = PQaX aik | ,右aik,k = 0 ,方程組無確定解，應(yīng)給出退出信息。b) 若a(kkk #0，則交換第ik行和k行，然后用Gauss消元法進行消元。四、MATLA歌件實現(xiàn)(1) 寫出Gauss消元法和列主元消元法實現(xiàn)的MA TLAB函數(shù)根據(jù)以上的算法，寫出如下程序：%Gauss元法 %function y=Gauss1(A,b)m,n=size(A);%檢查系數(shù)正確性if m=nerror( '矩陣A的行數(shù)和列數(shù)必須相同);return ;endif m=size(b)error( ' b的大小必須和A的行數(shù)或A的列數(shù)相同

5、');return ;end%再檢查方程是否存在唯一解if rank(A)=rank(A,b)error( ' A矩陣的秩和增廣矩陣的秩不相同，方程不存在唯一解');return ;end%這里采用增廣矩陣行變換的方式求解c=n+1;A(:,c)=b; %肖元過程 for k=1:n-1A(k+1:n, k:c)=A(k+1:n, k:c)-(A(k+1:n,k)/ A(k,k)*A(k, k:c); end%U代結(jié)果x=zeros(length(b),1);x(n)=A(n,c)/A(n,n);for k=n-1:-1:1x(k)=(A(k,c)-A(k,k+1:n)

6、*x(k+1:n)/A(k,k); end% 顯示計算結(jié)果%disp('x=');%disp(x);y=x;% %U主元消元法求解線，性方程組Ax=b%隅輸入矩陣系數(shù),b為方程組右端系數(shù) %方程組的解保存在x變量中function y=Gauss_line(A,b) format long ; %設(shè)置為長格式顯示，顯示15位小數(shù) m,n=size(A);% 先檢查系數(shù)正確性if m=nerror('矩陣A的行數(shù)和列數(shù)必須相同);return;end if m=size(b)error(' b的大小必須和A的行數(shù)或A的列數(shù)相同');return;end%再

7、檢查方程是否存在唯一解if rank(A)=rank(A,b)error(' A矩陣的秩和增廣矩陣的秩不相同，方程不存在唯一解');return;endc=n+1;A(:,c)=b;%(增廣)for k=1:n-1 r,m=max(abs(A(k:n,k);%i 主元m=m+k-1;%修正操作行的值if (A(m,k)=0) if (m=k)A(k m,:)=A(m k,:);瀕行endA(k+1:n, k:c)=A(k+1:n, k:c)-(A(k+1:n,k)/ A(k,k)*A(k, k:c);唯肖去end end x=zeros(length(b),1);%回代求解x(

8、n)=A(n,c)/A(n,n); for k=n-1:-1:1 x(k)=(A(k,c)-A(k,k+1:n)*x(k+1:n)/A(k,k); end y=x; format short ; %設(shè)置為默認格式顯示，顯示5位(2) 建立MATLAB界面利用MATLAB的GUI建立如下界面求解線性方程組：詳見程序。:、計算實例、數(shù)據(jù)、結(jié)果、分析下面我們對以上的結(jié)果進行測試，求解:一2100-3 | x 1一 10-3-4-1213|X251123-4X31-214149一13 一乂 一7 一G-L求君京WG求程*受輸入數(shù)據(jù)后點擊,得到如下結(jié)果：Gauss直接法求解結(jié)果Gauss列主元法求解結(jié)果

9、更改以上數(shù)據(jù)進行測試，求解如下方程組:一4得到如下結(jié)果:求解線性方程組Ax=b清按matlab中輸入矩陣方式輸入矩陣A4 33 1；3432343;1 234請按matlab中輸入矩陣方式輸入b11 -1 -1Gauss直接法求解結(jié)果-1-1.3323fl-0T62775Gf-0l7Gauss列主元法求解結(jié)果'福L術(shù)霍元-t. 3323e-016Copyright: yfg Tel、實驗中遇到的問題及解決辦法在本實驗中，遇到的問題主要有兩個：(1) 如何將上述的Gauss消元法的算法在 MATLAB中實現(xiàn)針對此問題我借鑒了網(wǎng)上以及課本上的算法的 MATLAB實

10、現(xiàn)的程序；(2) 如何將建立界面使得可以隨意輸入想要求解的相關(guān)矩陣后就可以直接求解針對此問題，我通過網(wǎng)上的一些關(guān)于MATLAB的GUI設(shè)計的相關(guān)資料，總結(jié)經(jīng)驗完成了此項任務(wù)。七、實驗結(jié)論通過以上的測試，我們發(fā)現(xiàn)以上算法和程序能夠求出線性方程組的比較精確解。八、參考文獻1 楊大地，王開榮.2006.數(shù)值分析.北京:科學出版社2 何光輝.2008.數(shù)值分析實驗.重慶大學數(shù)理學院數(shù)學實驗教學中心3 百度文庫，百度知道教師簽名年 月曰課程名稱數(shù)值分析實驗項目名 稱插值方法實驗項目類型驗證演示綜合設(shè)計其他指導教師何光輝成 績實驗目的：(1) 學會拉格朗日插值、牛頓插值等基本方法(2) 設(shè)計出相應(yīng)的算法，

11、編制相應(yīng)的函數(shù)子程序(3) 會用這些函數(shù)解決實際問題實驗內(nèi)容(1) 設(shè)計拉格朗日插值算法，編制并調(diào)試相應(yīng)的函數(shù)子程序(2) 設(shè)計牛頓插值算法，編制并調(diào)試相應(yīng)的函數(shù)子程序(3) 給定函數(shù)四個點的數(shù)據(jù)如下：X1.12.33.95.1Y3.8874.2764.6512.117試用拉格朗日插值確定函數(shù)在x=2.101, 4.234處的函數(shù)值。(4)已知1 =1,w'4 = 2,廿9 =3,用牛頓插值公式求 .5的近似值。三、實驗原理(1) 拉格朗日插值n 次拉格朗日插值多項式為：Ln(X)=y0l0(X)+y1l1(X)+y2l2(X)+- +ynln(X)n=1 時，稱為線性插值,L1(x)

12、=y0(x-x1)/(x0-x1)+ y1(x-x0)/(x1-x0)=y0+(y1-x0)(x-x0)/(x1-x0)n=2時，稱為二次插值或拋物線插值，精度相對高些1,而在其+ynln(x)L2(x)=y0(x-x1)(x-x2)/(x0-x1)/(x0-x2)+y1(x-x0)(x-x2)/(x1-x0)/(x1-x2)+y2(x-x0)(x-x1)/(x2-x0)/(x2-x1) 對節(jié)點xi(i=0,1,-,n)任一點xk(0<=k<=n)作一 n次多項式lk(xk),使它在該點上取值為余點 xi(i=0,1,-k,k+1, -,n止為 0,則插值多項式為Ln(x)=y0

13、l0(x)+y1l1(x)+y2l2(x)+-上式表明：n個點xi(i=0,1, - -k,k+1, -,n廊是lk(x)的零點。(2) 牛頓插值插商公式f(Xj)k和3*,.冬:j(Xj -X1).(Xj-Xj)(Xj -Xj).(Xj -Xk)Newton插值多項式為NJX) = f(Xi) fX,X2(X-Xi)fX,X2，X3(X-Xi)(X -X2)f Xi,X2,X3,.Xn(X X)(X X2).(X XnQ；四、MATLABC件實現(xiàn)(1) 分別寫出lagrange插值法和Newton插值法的求解函數(shù)%lagran 阿法求解函數(shù) % %x, y為初始數(shù)據(jù)，z為插值點 functi

14、on z=lagrange(x,y,a) format long ; % 顯示 15 位 n=length(x); % 取長度 %初始計算s=0;%進入公式計算for j=0:(n-1)t=1;for i=0:(n-1) if i=j t=t*(a-x(i+1)/(x(j+1)-x(i+1);endends=s+t*y(j+1);endz=s;%顯示輸出結(jié)果format short ;%NeWt姒求解函數(shù) %x, y為初始數(shù)據(jù)，z為插值點function j=Newton(x,y,z)n=max(size(x);l=1;a=y(1);B=a;s=1;%一次因子的乘積，預設(shè)為1dx=y;%li

15、商for i=1:n-1dx0=dx;for j=1:n-idx(j)=(dx0(j+1)-dx0(j)/(x(i+j)-x(j);enddf=dx(1);s=s*(z-x(i);%一次因子乘積a=a+s*df;%計算各次Newton插值的值l=l+1;B=a; %結(jié)果保存在變量 B中endj=B;(2)建立界面利用MATLAB中的GUI編程建立如下界面：請按matlab中的矩陣方式輸入Xi請按matlab中的矩陣方式輸入Yi輸入插值點x|-Lagrange插值結(jié)果:執(zhí)行插應(yīng)LNewton插值結(jié)果曲蝦邊立網(wǎng)Copyright: yfg Tel:見程序。計算實例、數(shù)據(jù)、

16、結(jié)果、分析 下面我們對以上的問題進行測試： 輸入數(shù)據(jù)：2 何光輝.2008.數(shù)值分析實驗.數(shù)理學院數(shù)學實驗教學中心3 百度文庫，百度知道教師簽名年 月曰課程名稱數(shù)值分析實驗項目名 稱數(shù)值微積分實驗項目類型驗證演示綜合設(shè)計其他指導教師何光輝成 績一、實驗目的：(1) 學會復化梯形、復化辛浦生求積公式的應(yīng)用(2) 設(shè)計出相應(yīng)的算法，編制相應(yīng)的函數(shù)子程序(3) 會用這些函數(shù)解決實際問題二、實驗內(nèi)容(1) 設(shè)計復化梯形公式求積算法，編制并調(diào)試相應(yīng)的函數(shù)子程序(2) 設(shè)計夏化辛浦生求積算法，編制并調(diào)試相應(yīng)的函數(shù)子程序(4) 分別用復化梯形公式和復化辛浦生公式計算定積分,i sin x , dx0 x三、

17、實驗原理(1)復化梯形求積公式:開始定義f (x)；給出a,b/輸入 n7, b a -hu ;xu a;T u 0nXi =a 十i h;Tu T 十 f (%);(i =1,n)Tu gf(a)+2T + f(b) 2*/輸出Tn/.，.吉束圖1復化梯形求積公式算法的流程圖Stepl給出被積函數(shù)f (x)、區(qū)問a,b端點a,b和等分數(shù)n ;b aStep2 求出 xk = kh, h=;nStep3 計算 f(a), f(b)* f(xQ ;k=0Step4 得Tn =1hf(a)+iff(xk) + f(b) 2k =0(2)復化辛普森求積公式圖2復化辛普森求積公式算法的流程圖Stepl

18、給出被積函數(shù)f(x)、區(qū)問a,b端點a,b和等分數(shù)n ；b aStep2 求出 xk = kh, h =;nStep3 計算 f(a), f(b),成 f(xQ,成 f (x)；k=sk=0k2n -1n -1Step4 得&=hf(a)+4£ f(x)+2£ f(xQ + f(b)6ek2 E四、 MATLA歌件實現(xiàn)(1) 分別寫出復化梯形和復化辛浦生求積的求解函數(shù)%«W 公式求積分值%function T=trap(f,a,b)%f為積分函數(shù) %a,b為積分區(qū)間%門是等分區(qū)間份數(shù) n=200; h=(b-a)/n;% 步長T=0; for k=1:(n

19、-1) x0=a+h*k; T=T+limit(f,x0); endT=h*(limit(f,a)+limit(f,b)/2+h*T;T=double(T);%SimpSOn求積分值 % function S=simpson(f,a,b) %f為積分函數(shù) %a,b為積分區(qū)間 %門是等分區(qū)間份數(shù) n=200;h=(b-a)/(2*n);% 步長s1=0; s2=0; for k=1:n x0=a+h*(2*k-1); s1=s1+limit(f,x0); end for k=1:(n-1) x0=a+h*2*k; s2=s2+limit(f,x0); end S=h*(limit(f,a)+li

20、mit(f,b)+4*s1+2*s2)/3; S=double(S);(2) 建立界面利用MATLAB中的GUI編程建立如下界面：在本實驗中，遇到的問題主要有兩個：(5) 如何將上述的積分算法在MATLAB中實現(xiàn)針對此問題我借鑒了網(wǎng)上以及課本上的算法的 MATLAB實現(xiàn)的程序；(6) 如何將建立界面使得可以隨意輸入想要求解的相關(guān)矩陣后就可以直接求解針對此問題，我通過網(wǎng)上的一些關(guān)于MATLAB的GUI設(shè)計的相關(guān)資料，總結(jié)經(jīng)驗完成了此項任務(wù)。七、實驗結(jié)論通過以上的測試，我們發(fā)現(xiàn)以上算法和程序能夠求出積分的比較精確解。八、參考文獻1 楊大地，王開榮.2006.數(shù)值分析.北京:科學出版社2 何光輝.2

21、008.數(shù)值分析實驗.數(shù)理學院數(shù)學實驗教學中心3 百度文庫，百度知道教師簽名年 月曰課程名稱數(shù)值分析實驗項目名 稱常微分方程的數(shù)值解法實驗項目類型驗證演示綜合設(shè)計其他指導教師何光輝成 績一、實驗目的：(1) 學會歐拉方法和四階龍格-庫塔方法的使用(2) 設(shè)計出相應(yīng)的算法，編制相應(yīng)的函數(shù)子程序(3) 會用這些函數(shù)解決實際問題二、實驗內(nèi)容用歐拉方法和四階龍格-庫塔方法求解微分方程初值1可題：y -sin(y)*x , y(0)=10,求y(1)。三、實驗原理(1)歐拉法歐拉法是解初值問題的最簡單的數(shù)值方法。從(9.2)式由于y (X0) = y°已給定，因而可以算出y'(x

22、76;) = f (xo, y°)設(shè)x = h充分小，則近似地有：y(x1)y(x0)頃xo) = f(x0,yo) h記 yi =y(xi) i =0,1，,n從而我們可以取y = y(o =hf(x°,y。)作為y (Xi)的近似值。利用yi及f (Xi, yi)又可以算出y(X2)的近似值:V2=V hf (Xi,yi)-般地，在任意點 Xn+i = (n + i)h處y(x)的近似值由下式給出yn d = yn hf ( xn，yn )這就是歐拉法的計算公式，h稱為步長。在實際計算時，可將歐拉法與梯形法則相結(jié)合，計算公式為：y：0i =yn hf(Xn,yn)k =

23、 0,i,2,.ynk：)= yn +?f(Xn,yn)+ f(Xz,yn%)】(2)四階龍格-庫塔方法四階龍格-庫塔法求解公式如下：h 一 一一 一yn# =yn + 二(ki +2k2 +2k3 +k4) 6ki = f (Xn,yn)(ih'«2 =f Xn+ 二 h, y- kiI22)(ih'k3 = f Xn +h, yn +k2 I k4 =f Xnh, ynhk3四、 MATLA歌件實現(xiàn)(i)分別寫出歐拉方法和四階龍格-庫塔方法求解微分方程的求解函數(shù)%Euiefc,初值 y(a)=c%function y=Euler(f,a,b,c)n=i000;h=

24、(b-a)/n;X=a:h:b;Y=zeros(i,n+i);Y(i)=c;for i=2:n+iX=X(i-i);y=Y(i-i);Y(i)=Y(i-i)+eval(f)*h;end%?-?a?t - ? - - - %function y=RK(f,a,b,c)i6n=1000;h=(b-a)/n;X=a:h:b;Y=zeros(1,n+1);Y(1)=c;for i=1:nx=X(i);y=Y(i);K1=h*eval(f);x=x+h/2;y=y+K1/2;K2=h*eval(f);x=x;y=Y(i)+K2/2;K3=h*eval(f);x=X(i)+h;y=Y(i)+K3;K4=h

25、*eval(f);Y(i+1)=Y(i)+(K1+2*K2+2*K3+K4)/6;end(2)建立界面利用MATLAB中的GUI編程建立如下界面:微分方程數(shù)值解(V=f(x,v)請輸入函數(shù)f 3v)請輸入初值點X。和所求點X1請輸入初值點X。的函數(shù)值y(X0)EulerS求解結(jié)果四階R-K公式求解結(jié)果詳見程序。97卸3Eu®法求解結(jié)果Ejler求程四階R-K公式求解結(jié)果六、實驗中遇到的問題及解決辦法在本實驗中，遇到的問題主要有兩個：(7) 如何將上述的求解微分方程的算法在MATLAB中實現(xiàn)針對此問題我借鑒了網(wǎng)上以及課本上的算法的 MATLAB實現(xiàn)的程序；(8) 如何將建立界面使得可以

26、隨意輸入想要求解的相關(guān)矩陣后就可以直接求解針對此問題，我通過網(wǎng)上的一些關(guān)于MATLAB的GUI設(shè)計的相關(guān)資料，總結(jié)經(jīng)驗完成了此項任務(wù)。七、實驗結(jié)論通過以上的測試，我們發(fā)現(xiàn)以上算法和程序能夠求出微分方程組的比較精確解。八、參考文獻1 楊大地，王開榮.2006.數(shù)值分析.北京:科學出版社2 何光輝.2008.數(shù)值分析實驗.數(shù)理學院數(shù)學實驗教學中心3 百度文庫，百度知道教師簽名年 月曰課程名稱數(shù)值分析實驗項目名 稱估計水塔的水流量實驗項目類型驗證演示綜合設(shè)計其他指導教師何光輝成 績一、實驗目的：(1) 學會對實際問題的分析方法(2) 學會利用所學的知識解決實際問題(3) 設(shè)計出相應(yīng)的算法，編制相應(yīng)的

27、應(yīng)用程序二、實驗內(nèi)容某居民區(qū)，其自來水是有一個圓柱形水塔提供，水塔高12.2m，塔的直徑為17.4m,水塔是由水泵根據(jù)水塔中的水位自動加水，一般水泵每天工作兩次。按照設(shè)計，當水塔中的水位降低至最低水 位，約8.2m時，水泵自動啟動加水。當水位升至最高水位，約 10.8m時，水泵停止工作。表略。三、實驗原理計算中將流量定義為單位時間流出的水的高度乘以水塔橫截面積。把時間分成5段：第1未供水段、水泵開啟第1段、第2未供水段、水泵開啟第2段、第3未 供水段。先直接對第1、2、3未供水段進行5次曲線擬合。再對得到的曲線分別求 導，取得流速(即單位時間內(nèi)流出的水的高度)。水泵開啟第1、2段，分別在兩端

28、各取兩個點，用時刻流速進行擬合得到這兩段的流速。流速乘以水塔橫截面積就得 到任何時刻的水流量。對其進行分段積分，求和得到一天的總水流量。四、MATLA歌件實現(xiàn)(1)程序function pushbutton1_Callback(hObject, eventdata, handles)% hObject handle to pushbutton1 (see GCBO)% eventdata reserved - to be defined in a future version of MATLAB% handles structure with handles and user data (se

29、e GUIDATA)figure(1);x=0,3316,6635,10619,13937,17921,21240,25223,28543,32284,39435,43318,46636,4995;y=31.75,31.10,30.54,29.94,29.55,28.92,28.50,27.87,27.52,26.97,35.50,34.45,33.50,32.67,31.56,30.81,30.12,29.27,28.42,27.67,26.97,34.75,33.89,33.40;t=x/3600;%M間單位為小時h=y/3.281;%水位高度為米x1=t(1:10);y1=h(1:10)

30、;f1=polyfit(x1,y1,5);t1=0:0.01:t(10);h1=polyval(f1,t1);plot(x1,y1, 'o' ,t1,h1,'k');xlabel( '時間(h)');ylabel( '水位(m)');title('第一階段供水的時間水位圖 ')function pushbutton2_Callback(hObject, eventdata, handles)% hObject handle to pushbutton2 (see GCBO)% eventdata reserved

31、- to be defined in a future version of MATLAB% handles structure with handles and user data (see GUIDATA)figure(2);x=0,3316,6635,10619,13937,17921,21240,25223,28543,32284,39435,43318,46636,4995;y=31.75,31.10,30.54,29.94,29.55,28.92,28.50,27.87,27.52,26.97,35.50,34.45,33.50,32.67,31.56,30.81,30.12,29

32、.27,28.42,27.67,26.97,34.75,33.89,33.40;t=x/3600;%M間單位為小時h=y/3.281;%水位高度為米x2=t(11:21);y2=h(11:21);f2=polyfit(x2,y2,5);t2=t(11):0.01:t(21);h2=polyval(f2,t2);plot(x2,y2,'o' ,t2,h2,'r');xlabel( '時間(h)');ylabel( '水位(m)');title('第二階段供水時段水位圖')% - Executes on button

33、press in pushbutton4.function pushbutton4_Callback(hObject, eventdata, handles)% hObject handle to pushbutton4 (see GCBO)% eventdata reserved - to be defined in a future version of MATLAB% handles structure with handles and user data (see GUIDATA)figure(3);x=0,3316,6635,10619,13937,17921,21240,25223

34、,28543,32284,39435,43318,46636,4995;y=31.75,31.10,30.54,29.94,29.55,28.92,28.50,27.87,27.52,26.97,35.50,34.45,33.50,32.67,31.56,30.81,30.12,29.27,28.42,27.67,26.97,34.75,33.89,33.40;t=x/3600;h=y/3.281;x3=t(22:24);y3=h(22:24);f3=polyfit(x3,y3,5);t3=t(22):0.01:t(24);h3=polyval(f3,t3);plot(x3,y3, 'o' ,t3,h3,'r');xlabel( '時間(h)');ylabel( '水位(m)');title( '第三階段的時間水位圖')% - Executes on button press in pushbutton7.functionpus