郭德泉數值分析

1、青島農業(yè)大學 本科生課程論文題 目： 數值分析課程設計 姓 名： 郭德泉 學 院： 理學與信息科學學院 專 業(yè)： 信息與計算科學 班 級： 12級01班 學 號： 20125990 指導教師： 常桂娟 完成時間： 2014年11月30日 二一四年十二月一日 課 程 論 文 任 務 書學生姓名 郭德泉 指導教師 常桂娟 論文題目 數值分析課程設計 論文內容（需明確列出研究的問題）： 運用數值分析的知識 來解決數學研究問題中的計算問題。（1）列主元Gauss消去法解方程組。（2）雅克比迭代法和高斯迭代法求線性方程組的解。（3）求解拉格朗日插值多項式和牛頓插值多項式。 資料、數據、技術水平等方面的要

2、求：論文要符合一般學術論文的寫作規(guī)范，具備學術性、科學性和一定的創(chuàng)造性。文字要流暢、語言要準確、論點要清楚、論據要準確、論證要完整、嚴密，有獨立的觀點和見解。內容要理論聯系實際，計算數據要求準確，涉及到他人的觀點、統(tǒng)計數據或計算公式等要標明出處，結論要寫的概括簡短。參考文獻的書寫按論文中引用的先后順序連續(xù)編碼。（根據情況修改） 發(fā)出任務書日期2014.11.20 完成論文（設計）日期 2014.11.30 學科組或教研室意見（簽字） 院、系（系）主任意見（簽字） 注：此表裝訂在課程論文之前。 4 目錄【摘要】·······&

3、#183;·················································&

4、#183;·······1【關鍵詞】·········································&

5、#183;······················1 【前言】··························&

6、#183;······································2一、列主元Gauss消去法·········

7、;··········································3（一）理論知識：······&#

8、183;················································3（二）題目：&

9、#183;·················································&

10、#183;·······4（三）程序清單：········································

11、3;·············4 （四）程序運行操作過程與輸出結果：·································

12、83;···5 （五）對計算過程與結果的分析：·········································7（六）在課程設計中的心得體會：

13、·········································7二、雅克比迭代法和高斯-賽德爾迭代法······

14、;·································8（一）理論知識：···············&#

15、183;·······································8（二）題目：·········&

16、#183;················································8（三）程序清

17、單：·················································

18、3;···10（四）程序運行操作過程與輸出結果：····································10（五）對計算過程與結果的分析：····

19、3;···································13 （六）在課程設計中的心得體會············

20、;·····························13三、求拉格朗日插值多項式和牛頓插值多項式·················&

21、#183;···············14（一）基礎知識：································

22、83;·····················14（二）題目：···························&

23、#183;······························15（三）程序清單：·················

24、83;···································16（四）程序運行操作過程與輸出結果：···········&#

25、183;························18 （五）對計算過程與結果的分析：······················&#

26、183;·················18 （六）在課程設計中的心得體會：·····························&#

27、183;··········19四、參考文獻······································

28、·····················19五、課程設計自我評價···························

29、·························20 數值分析課程設計 信息與計算科學專業(yè) 郭德泉 指導教師 常桂娟【摘要】 數值分析是研究各種數學問題求解的數值計算方法，是數學中計算數學分支的重要內容。近幾十年來，隨著計算機的飛速發(fā)展，數值計算方法的學習與研究越來越離不開計算機。實際計算中遇到的數值問題只有與計算機相結合，算法與程序密切聯系，形

30、成切實可靠的數值軟件才能為社會創(chuàng)造更大的社會財富。本文主要描述運用數值分析的知識來解決數學研究問題中的計算問題，包括運用列主元高斯消去法求方程組的的解；運用雅可比迭代法和高斯-塞德爾迭代法求線性方程組的解；求四次拉格朗日插值多項式及牛頓插值多項式?！娟P鍵詞】 gauss消去法 雅克比迭代法 高斯-賽德爾迭代法 lagrange插值 Newton插值 前言 本文主要介紹要MATLAB軟件完成對各種數學問題的求解，數值分析的問題求解過程會有一定的復雜性，通過MATLAB軟件計算，提高了運算速度，避開了復雜的運算過程，對我們的學習有很大的幫助。數值分析是一門基礎課,它象通常的數學課程一樣有自身嚴密的

31、科學體系,但它又是一門應用性很強的課程,目的是使學生能夠用本課程的理論在計算機上實現有關的科學與工程計算,計算能力的培養(yǎng)對工科各專業(yè)的學生都是十分重要的,數值分析課程的課程設計環(huán)節(jié)正是為了適應這種需要而設置的,使得學生在學習各種算法時,能在課程設計環(huán)節(jié)中深刻體會算法的內涵、物理背景和實際意義,同時也能提高學生學習該課程的興趣, 而更重要的是能提高學生綜合運用知識的能力和培養(yǎng)學生利用計算機解決實際問題的能力。 現在，讓我們一起學習用MATLAB軟件解決一下 數學問題吧。 一、列主元Gauss消去法（一）理論知識： 1.列主元Gauss消去法就是通過矩陣的行變換達到消元的目的，從而將方程組的系數矩

32、陣由對稱矩陣變?yōu)槿蔷仃?，最后獲得方程組的解。定理：當一個方陣A的1-（n-1）階順序主子式都不為0時，A存在唯一的doolittle(杜利特爾)分解。對于滿足條件的這個方陣，采用直接三角分解法可以得到這個唯一的L、U矩陣。但事實上，matlab的lu函數并不是采用直接三角分解法，因為直接三角分解法存在著引入很大舍入誤差的可能（直接分解法過程中，當對角元為零或很小時）。matlab的lu函數采用的是列主元三角分解法，又稱為PLU分解，實現LU=PA。該分解法要求的前提條件是A為非奇異矩陣 2.使用matlab自帶命令lu函數. 將系數矩陣A轉變成等價兩個矩陣L和U的乘積 ，其中L和U分別是下三

33、角和上三角矩陣。當A的所有順序主子式都不為0時，矩陣A可以分解為A=LU，其中L是單位下三角矩陣，U是上三角矩陣。L，U=lu(X)：產生一個上三角矩陣U和一個“心理上下三角矩陣”L（即由下三角矩陣和置換矩陣構成），使得X=L*U。L,U,P=lu(X);產生一個上三角矩陣U和一個下三角矩陣L，以及一個置換矩陣P，使得P*X=L*U 。其中X必須是方陣（2） 題目： 1. 用列主元Gauss消去法解方程組（3） 程序清單： 1.列主元Gauss消去法MATLAB程序如下：%magauss2.mfunction x=magauss2(A,b,flag)%用途：列主元Gauss消去法解線性方程組A

34、x=b%格式：x=magauss(A,b,flag),A為系數矩陣，b為右端項，若flag=0，則不顯示中間過程，否則顯示中間過程，默認為0，x為解向量if nargin<3,flag=0;endn=length(b);for k=1:(n-1) % 選主元 ap,p=max(abs(A(k:n,k); p=p+k-1; if p>k t=A(k,:);A(k,:)=A(p,:);A(p,:)=t; t=b(k);b(k)=b(p);b(p)=t; end %消元 m=A(k+1:n,k)/A(k,k); A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-m*A(k,k+

35、1:n); b(k+1:n)=b(k+1:n)-m*b(k);A(k+1:n,k)=zeros(n-k,1); if flag=0, Ab=A,b,endend%回代x=zeros(n,1);x(n)=b(n)/A(n,n);for k=n-1:-1:1 x(k)=(b(k)-A(k,k+1:n)*x(k+1:n)/A(k,k);end 2.輸入>> A=2 2 3;4 7 7;-2 4 5; %用數組形式將A的值輸入>> b=3 1 -7' %用數組形式將b的值輸入>> L,U,P=lu(A) %產生一個上三角矩陣U和一個下三角矩陣L（4） 程序運

36、行操作過程及輸出結果 （1）在MATLAB命令窗口執(zhí)行>> A=2 2 3;4 7 7;-2 4 5;>> b=3 1 -7'>> x=magauss2(A,b); x'得到計算結果：ans =2.0000 -2.0000 1.0000 （2）>> A=2 2 3;4 7 7;-2 4 5; %用數組形式將A的值輸入>> b=3 1 -7' %用數組形式將b的值輸入>> L,U,P=lu(A) %產生一個上三角矩陣U和一個下三角矩陣LL = %產生一個下三角矩陣L 1.0000 0 0 -0.500

37、0 1.0000 0 0.5000 -0.2000 1.0000U = %產生一個上三角矩陣U 4.0000 7.0000 7.0000 0 7.5000 8.5000 0 0 1.2000P = 0 1 0 0 0 1 1 0 0 >> y=L(P*b)y = 1.0000 -6.5000 1.2000>> x=Uyx = 2.0000 -2.0000 1.0000（五）對計算過程及結果的分析：通過檢驗可以得出結果是準確的，正確反應了實驗程序編寫的正確性。并且Gauss消去法的消元過程不至于中斷，避免了普通高斯消元法中遇到某個主元為零或者當主元絕對值很小時，計算將會停

38、止或求出的結果將與其實際結果相差很遠的情況，并且舍入誤差較小。另外利用matlab自帶命令lu函數，來幫助解決此問題會更加快捷。（6） 在課程設計中的心得體會 矩陣的轉置符號要用英文輸入法下的單引號代替，否則將會出現錯誤，運行不出結果。 二、雅可比迭代法和高斯-塞德爾迭代法（1） 理論知識： 1、.diag(A)是提取出矩陣A的主對角線元素,得到的是一維向量,diag(diag(A)是一個對角矩陣.如A=1,2,3;4,5,6;7,8,9,則diag(A)=1 5 9; 而diag(diag(A)=1,0,0;0,5,0;0,0,9 2. tril(A)產生A矩陣的下三角矩陣，其余元素補0 3

39、. triu(A)產生A矩陣的上三角矩陣，其余元素補0 4. inv(A)求A矩陣的逆矩陣 5. norm(X,p) 功能:norm函數可計算幾種不同類型的向量范數,根據p的不同可得到不同的范數 （二）題目： 用雅可比迭代法和高斯-塞德爾迭代法求線性方程組的解，初始迭代值，求近似解使其滿足（3） 程序清單： 1、Jacobi程序清單i=1;a=0.9889 0.0005 0.0002 0.3372;-0.0046 0.9946 0.0077 0.5381;-0.0002 0.0092 0.9941 0.1936;-0.8002 0.5671 0.1849 0.3891;d=diag(diag(

40、a);l=d-tril(a);u=d-triu(a);d0=inv(d);b=1;0;1;0;x0=zeros(4,1);B=d0*(l+u);f=d0*b;x=B*x0+f;while norm(x-x0,inf)>=1e-6 %10的負6次方 x0=x; x=B*x0+f; i=i+1;endxi 2、Gauss-Seidel程序清單i=1;a=0.9889 0.0005 0.0002 0.3372;-0.0046 0.9946 0.0077 0.5381;-0.0002 0.0092 0.9941 0.1936;-0.8002 0.5671 0.1849 0.3891;d=diag

41、(diag(a);l=d-tril(a);u=d-triu(a);b=1;0;1;0;x0=zeros(4,1);B=inv(d-l)*u;f=inv(d-l)*b;x=B*x0+f; while norm(x-x0,inf)>=1e-6 x0=x; x=B*x0+f; i=i+1;endxI （四）程序運行操作過程與輸出結果 1、在MATLAB中輸入>> i=1;a=0.9889 0.0005 0.0002 0.3372;-0.0046 0.9946 0.0077 0.5381;-0.0002 0.0092 0.9941 0.1936;-0.8002 0.5671 0.18

42、49 0.3891;d=diag(diag(a);l=d-tril(a);u=d-triu(a);d0=inv(d);b=1;0;1;0;x0=zeros(4,1);B=d0*(l+u);f=d0*b;x=B*x0+f;while norm(x-x0,inf)>=1e-6 %10的負6次方 x0=x; x=B*x0+f; i=i+1;endxi點擊Enter結果如下x = 0.3454 -1.0604 0.6353 1.9539i =192、 在MATLAB中輸入 >> i=1;a=0.9889 0.0005 0.0002 0.3372;-0.0046 0.9946 0.00

43、77 0.5381;-0.0002 0.0092 0.9941 0.1936;-0.8002 0.5671 0.1849 0.3891;d=diag(diag(a);l=d-tril(a);u=d-triu(a);b=1;0;1;0;x0=zeros(4,1);B=inv(d-l)*u;f=inv(d-l)*b;x=B*x0+f; while norm(x-x0,inf)>=1e-6 x0=x; x=B*x0+f; i=i+1;endxi按Enter鍵得出結果x = 0.3454 -1.0604 0.6353 1.9539i =10 （五）對計算過程及結果的分析： 該問題的數據都在小數點

44、后面五位了，如果用紙筆計算的話運算起來將會相當麻煩，這次就更顯示出了，運用MATLAB軟件計算的有效性，而且雖然是大數據運算起來準確性還是可以保證的。比較兩種方法，雅可比迭代法在解決此問題上會較為精確些。 （六）在課程設計中的心得體會 問題中的數據比較多，且比較長要細心輸入，否則容易出錯 。第二，第一次看到，兩種方法的運行結果不一致時，感覺有些詫異，但一分析這數據的精確度問題，和方法問題，明白這誤差是可以理解的。 三、求拉格朗日插值多項式及牛頓插值多項式（1） 理論知識： 定義1： 設為區(qū)間上函數，為上個互不相同的點，為給定的某一函數類，若上有函數，滿足. (1)則稱為關于節(jié)點在上的插值函數，

45、稱點為插值節(jié)點；稱為插值型值點，簡稱型值點或插值點;稱為被插函數.定義2： 已知函數在區(qū)間上的個點的值，即已知，尋求一個解析形式的函數，使之滿足 (2)則稱為插值結點,為被插值函數，為插值函數，稱條件為插值條件，若為次數不超過的多項式，即， . (3)其中為實數，則稱為插值多項式.定理1:在個相異結點滿足插值條件而次數不高于的多項式是唯一的.拉格朗日插值多項式定義：對某個多項式函數，已知有給定的k + 1個取值點：其中xj對應著自變量的位置，而yj對應著函數在這個位置的取值。假設任意兩個不同的xj都互不相同，那么應用拉格朗日插值公式所得到的拉格朗日插值多項式為：其中每個為拉

46、格朗日基本多項式（或稱插值基函數），其表達式為：拉格朗日基本多項式的特點是在xj 上取值為1，在其它的點 上取值0。求拉格朗日插值多項式的算法：輸入：數據對應的坐標點構成的向量x，y和要求的點x0；輸出：拉格朗日插值多項式已經特值對應的y值Step1: 對于i=1:n 對于j=1：n 如果j不等于i，執(zhí)行step2； 執(zhí)行step3Step2: l=l*(t-x(j)/(x(i)-x(j);/這就是計算出的拉格朗日基本多項式step3：f=f+l;/計算插值多項式（2） 題目：根據下表所列的數據點求出四次拉格朗日插值多項式及牛頓插值多項式，并計算當x=2.0時的值。0.20.40.60.81.

47、00.980.920.810.640.38（3） 程序清單：1、拉格朗日插值function f=Lagrange(x,y,x0)% 求已知數據點的拉格朗日插值多項式% 已知數據點的x坐標向量:x% 已知數據點的y坐標向量:y%插值的x坐標:x0%求得的拉格朗日插值多項式在x0處的插值:fsyms t; %定義符號變量x=0.2 0.4 0.6 0.8 1.0;y=0.98 0.92 0.81 0.64 0.38;if(length(x)=length(y) n=length(x);else disp('x和y的維數不相等!'); return;end %檢驗錯誤f=0.0;for(i=1:n) l=y(i); for(j=1:i-1) l=l*(t-x(j)/(x(i)-x(j); end; for(j=i+1:n) l=l*(t-x(j)/(x(i)-x(j);%計算拉格朗日基函數 end; f=f+l; %計算拉格朗日插值函數 simplify(f); %化簡 if(i=n) if(nargin=3) % nargin是用來判斷輸入變量個數的函數 f=subs(f,'t',x0); %計算插值點的函數值,把xo賦給t帶入f中計算,結果賦給f else f=collect(f); %將