數(shù)值分析課程設(shè)計含代碼

1、 成 績 評 定 表學生姓名班級學號專 業(yè)信息與計算科學課程設(shè)計題目數(shù)值分析算法案例評語組長簽字：成績?nèi)掌?20 年 月 日課程設(shè)計任務(wù)書學 院理學院專 業(yè)信息與計算科學學生姓名班級學號課程設(shè)計題目數(shù)值分析算法案例實踐教學要求與任務(wù):要求：格式以學校畢業(yè)論文格式要求為準，不準粘貼圖片，尤其公式。對每個試驗，要求有：實驗基本原理，實驗?zāi)康?，實驗?nèi)容及數(shù)據(jù)來源和實驗結(jié)論。以班級為單位統(tǒng)一裝訂封皮。6月25日，十八周的周二交論文每人至少四個實驗，最少15頁任務(wù)（實驗項目）：線性方程組數(shù)值解法 參考題目：（1) 列主元Gauss消去法（2）LU分解法插值法和數(shù)據(jù)擬合 參考題目：（1）Lagrange插

2、值（2）Newton插值（3）最小二乘擬合數(shù)值積分 參考題目：（1）復(fù)化Simpon積分（2）變步長的梯形積分公式（3）龍貝格求積公式常微分方程數(shù)值解 Runge-Kutta方法數(shù)值方法實際應(yīng)用 用數(shù)值方法解決實際問題（自選）工作計劃與進度安排:線性方程組數(shù)值解法 （4學時）插值法和數(shù)據(jù)擬合 （4學時）數(shù)值積分 （4學時）常微分方程數(shù)值解 （4學時）數(shù)值方法實際應(yīng)用 （4學時）答辯 （4學時）指導(dǎo)教師： 201 年 月 日專業(yè)負責人：201 年 月 日學院教學副院長：201 年 月 日摘 要實驗方法與理論方法是推動科學技術(shù)發(fā)展的兩大基本方法，但有局限性。許多研究對象，由于空間或時間的限制，既不

3、可能用理論精確描述，也不能用實驗手段實現(xiàn)。 數(shù)值模擬或稱為科學計算突破了實驗和理論科學的局限，在科技發(fā)展中起到越來越重要的作用?？梢哉J為，科學計算已于實驗、理論一起成為科學方法上不可或缺的三個主要手段。 計算數(shù)學的研究是科學計算的主要組成部分，而數(shù)值分析則是計算數(shù)學的核心。數(shù)值計算是研究使用計算機來解決各種數(shù)學問題的近似計算方法與理論，其任務(wù)是提供在計算機上可解的、理論可靠的、計算復(fù)雜性低的各種常用算法。數(shù)值分析的主要內(nèi)容：1）、數(shù)值代數(shù)：求解線性和非線性方程組的解，分直接方法和間接方法兩大類；2）、插值、曲線擬合和數(shù)值逼近；3）、數(shù)值微分和數(shù)值積分；4）、常微分和偏微分方程數(shù)值解法。 本文主

4、要通過Matlab軟件，對數(shù)值分析中的一些問題進行求解，如列主元Gauss消去法，Lagrange插值多項式，復(fù)化Simpson公式，Runge-Kutta方法以及數(shù)值分析在實際問題中的應(yīng)用，并在求解的過程中更加熟識這門課程的主要內(nèi)容，以及加強對課程知識的掌握。在學習與設(shè)計計算方法時，從數(shù)學理論角度，學會分析方法的誤差、收斂性和穩(wěn)定性，保證計算方法的準確性；從實際應(yīng)用的角度出發(fā)，掌握計算方法的結(jié)構(gòu)與流程，能夠把計算方法轉(zhuǎn)換為可在計算機上直接處理的程序，保證算法的可用性。關(guān)鍵詞：列主元Gauss消去法；Lagrange插值；復(fù)化Simpson公式；Runge-Kutta目 錄實驗一 列主元Gau

5、ss消去法11.1 實驗?zāi)康?1.2 基本原理11.3 實驗內(nèi)容21.4 實驗結(jié)論3實驗二 拉格朗日插值多項式42.1 實驗?zāi)康?2.2 基本原理42.3 實驗內(nèi)容42.4 實驗結(jié)論9實驗三 復(fù)化Simpson求積公式103.1 實驗?zāi)康?03.2 基本原理103.3 實驗內(nèi)容103.4 實驗結(jié)論12實驗四 龍格-庫塔(Runge-Kutta)方法134.1 實驗?zāi)康?34.2 基本原理134.3 實驗內(nèi)容144.4 實驗結(jié)論15實驗五 數(shù)值方法實際應(yīng)用165.1 實驗?zāi)康?65.2 基本原理165.3 實驗內(nèi)容165.4 實驗結(jié)論22參考文獻2322實驗一 列主元Gauss消去法1.1 實驗

6、目的1) 理解列主元消去法的原理；2) 熟悉列主元消去法的計算步驟，能用代碼編寫；3) 解決實際問題。1.2 基本原理在順序Gauss消去法中，必須要求；否則無法進行計算。即使，但其絕對值很小，由于舍入誤差的影響，也可能會引起很大的誤差，從而使上述方法失效。為了使消元過程中減小舍入誤差和不至于中斷，可以按照不同的自然順序進行消元。在第k步消元時，增廣矩陣為 （1.1）不一定選取作為主元，而從同列中選取絕對值最大的作為主元素，即 （1.2） 若，此時矩陣不可逆，方程的解不確定，則停止計算；否則，當r>k時，則其增廣矩陣中交換第k行和第r行，即 （1.3）使成為主元。然后再按Gauss消去法

7、進行消元運算。于是就得到列主元Gauss消去法。1.3 實驗內(nèi)容1.3.1 程序來源 首先建立一個gaussMethod.m的文件，用來實現(xiàn)列主元的消去方法。文件內(nèi)容如下：function x=gaussMethod(A,b)%高斯列主元消去法，要求系數(shù)矩陣非奇異的n = size(A,1);if abs(det(A)<= 1e-8 error('系數(shù)矩陣是奇異的'); return;endfor k=1:n ak = max(abs(A(k:n,k); index = find(A(:,k)=ak); if length(index) = 0 index = find(

8、A(:,k)=-ak); end %交換列主元 temp = A(index,:); A(index,:) = A(k,:); A(k,:) = temp; temp = b(index);b(index) = b(k); b(k) = temp; %消元過程 for i=k+1:n m=A(i,k)/A(k,k); %消除列元素 A(i,k+1:n)=A(i,k+1:n)-m*A(k,k+1:n); b(i)=b(i)-m*b(k); endend %回代過程 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)&#

9、39;)/A(k,k); end; end然后調(diào)用gaussMethod函數(shù)，來實現(xiàn)列主元的高斯消去法。建立一個文件gauss，內(nèi)容如下：cleardisp('*')x=gaussMethod(input('請輸入系數(shù)矩陣：'),input('請輸入常數(shù)列：')disp('*')1.3.2 實例分析例：在Matlab上，利用列主元法求線性方程組的解: 解：運行程序，按照提示輸入方程的系數(shù)矩陣及常數(shù)列，如下所示：*請輸入系數(shù)矩陣：1 2 1 4;2 0 4 3;4 2 2 1;-3 1 3 2請輸入常數(shù)列：13;28;20;6x =

10、 3 -1 4 2* 即該方程的解為：1.4 實驗結(jié)論把向量計算得到的解帶入方程組，驗證正確性，和其他的方法比較，列主元具有一定的簡單性，比較容易實現(xiàn)。避免使用其他方法的誤差或不能進行性。而列主元也有一定的限制，要求行列式的值不為0。實驗二 拉格朗日插值多項式2.1 實驗?zāi)康?）熟悉簡單的拉格朗日插值多項式的基本概念；2）熟悉Lagrange公式及源代碼，會利用它來計算基本函數(shù)；3）能構(gòu)造出正確的插值多項式；2.2 基本原理設(shè)函數(shù)在區(qū)間a,b上有定義，且已知在點上的值若存在一個次數(shù)不超過n的多項式 （2.1）使其滿足 （2.2）則稱為的n次插值多項式，稱點為插值節(jié)點，稱條件（2.2）為插值條件

11、。包含插值節(jié)點的區(qū)間成為插值區(qū)間。通過平面上不同的兩點可以確定一條直線經(jīng)過這兩點，就是拉格朗日線性插值問題，對于不在同一直線的三點得到的插值多項式為拋物線。拉格朗日是比較基礎(chǔ)的方法，本身比較容易實現(xiàn)，容易理解。給定n+1個不同節(jié)點，構(gòu)造的n次拉格朗日插值多項式： ， （2.3）2.3 實驗內(nèi)容2.3.1 程序來源首先建立一個Lagrange.m的文件，用來實現(xiàn)Lagrange插值。文件內(nèi)容如下：%輸入：x是插值節(jié)點橫坐標向量；y是插值節(jié)點對應(yīng)縱坐標向量%輸出：C是拉格朗日插值多項式的系數(shù)矩陣；L是插值函數(shù)系數(shù)矩陣functionC,L=Lagrange(x,y)w=length(x);n=w-

12、1;L=zeros(w,w);for k=1:n+1 V=1; for j=1:n+1 if k=j V=conv(V,poly(x(j)/(x(k)-x(j); end end L(k,:)=V;endC=y*L然后調(diào)用Lagrange函數(shù)，來實現(xiàn)Lagrange插值法。建立一個文件Lg，內(nèi)容如下：cleardisp('*')x=input('請輸入已知點的橫坐標組：');y=input('請輸入已知點的縱坐標：');C,L=Lagrange(x,y);yi=polyval(C,input('請輸入需要計算得橫坐標組：')xx=

13、1.5:0.05:6.5;yy=polyval(C,xx);plot(xx,yy,x,y,'o')disp('*')2.3.2 實例分析例1 有4對數(shù)據(jù)（1.6,3.3），（2.7,4.22），（3.9,5.61），（5.6,2.94），寫出這4個數(shù)據(jù)點的Lagrange插值公式，并計算出橫坐標組xi=2.101,4.234時對應(yīng)的縱坐標值。解：4個數(shù)據(jù)點的Lagrange插值公式為： 運行程序，按照提示輸入已知點的橫坐標組、縱坐標組及需要計算得橫坐標組，如下所示：*請輸入已知點的橫坐標組：1.6,2.7,3.9,5.6請輸入已知點的縱坐標：3.3,1.22,5

14、.61,2.94C = -1.0539 11.0551 -34.4933 34.5053請輸入需要計算得橫坐標組：2.101,4.234yi = 1.0596 6.6457*即 輸出圖形：圖2.1 輸出擬合曲線例2 將區(qū)間-5，5等分5份、10份，求函數(shù)拉格朗日差值多項式，做出函數(shù)原圖像，觀察龍格現(xiàn)象。解：首先將區(qū)間五等分，取各端點坐標擬合曲線，輸入：clearx=-5:2:5;y=1./(1+x.2);C,L=Lagrange(x,y);輸出結(jié)果：C = 0.0000 0.0019 -0.0000 -0.0692 0.0000 0.5673即 然后輸出圖形，比較擬合效果。在matlab中輸入

15、：xx=-5:0.1:5;yy=polyval(C,xx);hold onplot(xx,yy,'-',x,y,'o')xp=-5:0.01:5;z=1./(1+xp.2);plot(xp,z,'r')legend('拉格朗日插值曲線','插值點','原曲線')輸出五等分插值龍格現(xiàn)象圖形：圖2.2 五等分插值龍格現(xiàn)象圖形 十等分的過程與五等分基本相似，如下輸入：clear,clfx=-5:1:5;y=1./(1+x.2);C,L=Lagrange(x,y);輸出：C = -0.0000 -0.000

16、0 0.0013 -0.0000 -0.0244 0.0000 0.1974 -0.0000 -0.6742 -0.0000 1.0000即 然后輸出圖形，比較擬合效果。在matlab中輸入：xx=-5:0.1:5;yy=polyval(C,xx);hold onplot(xx,yy,'-',x,y,'o')xp=-5:0.01:5;z=1./(1+xp.2);plot(xp,z,'r') legend('拉格朗日插值曲線','插值點','原曲線')輸出十等分插值圖形：圖2.3 十等分插值龍格現(xiàn)象圖

17、形2.4 實驗結(jié)論通過圖像可以得出：1）并不是插值節(jié)點越多，插值多項式逼近函數(shù)效果就越好；2）誤差較大的地方，是在插值區(qū)間兩端點附近出現(xiàn)；3）求插值多項式，不需要求解線性方程組，當數(shù)據(jù)比較多時，此公式才能 顯示出優(yōu)越性；4）函數(shù)值可以用符號形式表示，數(shù)據(jù)點未確定的縱坐標可以用多項式表示；實驗三 復(fù)化Simpson求積公式3.1 實驗?zāi)康?）了解單獨的Simpson公式原理及幾何意義；2）學會復(fù)化Simpson求積公式的編程與應(yīng)用；3）掌握Mat lab提供的計算積分的各種函數(shù)的使用方法；3.2 基本原理 復(fù)化Simpson公式是一種比較實用的積分方法，可以給出誤差估計。首先將區(qū)間a,b N等分

18、，子區(qū)間的長度為 （3.1）在每個子區(qū)間上采用Simpson公式，在用Simpson公式時，還要將子區(qū)間再二等分，因此有2N+1個分點。即 （3.2）經(jīng)推導(dǎo)得到， （3.3）稱為為復(fù)化Simpson值，稱（3.3）式為復(fù)化Simpson公式。3.3 實驗內(nèi)容3.3.1 程序來源 編寫復(fù)化Simpson求積函數(shù)（函數(shù)名：s_quad.m）Function I=S_quad(x,y)% 復(fù)化求積公式% x為被積函數(shù)自變量的等距節(jié)點；y為被積函數(shù)在節(jié)點處的函數(shù)值。n=length(x);m=length(y); % 積分自變量的節(jié)點數(shù)應(yīng)與它的函數(shù)值個數(shù)相同；if n=m error ('Th

19、e length of X and Y must be equal'); return;endif rem(n-1,2)=0 % 如果n-1不能被2整除，則調(diào)用復(fù)化公式 error ('節(jié)點數(shù)不滿足要求'); return;endN=(n-1)/2;h=(x(n)-x(1)/N;a=zeros (1,n);for k=1:Na(2*k-1)=a(2*k-1)+1;a(2*k)=a(2*k)+4; a(2*k+1)=a(2*k+1)+1;endI=h/6*sum(a.*y);然后調(diào)用s_quad函數(shù)，來實現(xiàn)復(fù)化Simpson公式法。建立一個文件SPS，內(nèi)容如下：clear

20、disp('*')x=input('請輸入積分上下限及點間的間隔（例如-1:0.1:1）：');y=input('請輸入被積公式：y=');I=S_quad(x,y);disp('得出積分值I=')disp(I);disp('*')3.3.2 實例分析：例1 用復(fù)化Simpson公式求積分，在積分區(qū)間中點與點之間的間隔取為 0.1。解：運行程序，按照提示輸入積分上下限、點間的間隔及被積公式，如下所示：*請輸入積分上下限及點間的間隔（例如-1:0.1:1）：-1:0.1:1請輸入被積公式：y=exp(-x.2)得出積

21、分值I= 1.4936* 真值為：1.4937例2 計算積分，將區(qū)間8等分。解：運行程序，按照提示輸入積分上下限、等分后的區(qū)間長度及被積公式，如下 所示：*請輸入積分上下限及點間的間隔（例如-1:0.1:1）：0:0.125:1請輸入被積公式：y=x./(4+x.2)得出積分值I= 0.1116* 真值為：0.1115723.4 實驗結(jié)論復(fù)化Simpson計算所得的結(jié)果誤差較小，精度較高，更適合科學計算與應(yīng)用，且公式具有收斂性，穩(wěn)定性良好。實驗四 龍格-庫塔(Runge-Kutta)方法4.1 實驗?zāi)康?）熟悉Runge-Kutta常微分方程初值問題的基本原理2）了解Runge-Kutta常微

22、分方程初值問題的計算3）能編程實現(xiàn)Runge-Kutta常微分方程初值問題4.2 基本原理 龍格-庫塔(Runge-Kutta)方法是一種在工程上應(yīng)用廣泛的高精度單步算法。由于此算法精度高，采取措施對誤差進行抑制，所以其實現(xiàn)原理也較復(fù)雜。該算法是構(gòu)建在數(shù)學支持的基礎(chǔ)之上的。對于一階精度的歐拉公式有： (4.1)當用點處的斜率近似值與右端點處的斜率的算術(shù)平均值作為平均斜率的近似值，那么就會得到二階精度的改進歐拉公式： (4.2)依次類推，如果在區(qū)間內(nèi)多預(yù)估幾個點上的斜率值，并用他們的加權(quán)平均數(shù)作為平均斜率的近似值，顯然能構(gòu)造出具有很高精度的高階計算公式。經(jīng)數(shù)學推導(dǎo)、求解，可以得四階龍格庫塔公式，

23、也就是在工程中應(yīng)用廣泛的經(jīng)典龍格庫塔算法。 四階龍格庫塔方法基本公式： (4.3)4.3 實驗內(nèi)容4.3.1 程序來源寫出四階Runge-Kutta方法源代碼的得.m文件：function R=rk4(f,a,b,ya,M)% a,b分別為左右端點，ya為初始值，M為步長幅數(shù)。h=(b-a)/M;T=zeros(1,M+1);Y=zeros(1,M+1);T=a:h:b;Y(1)=ya;for j=1:M k1=feval(f,T(j),Y(j); k2=feval(f,T(j)+h/2,Y(j)+k1*h/2); k3=feval(f,T(j)+h/2,Y(j)+k2*h/2);k4=fev

24、al(f,T(j)+h,Y(j)+h*k3); Y(j+1)=Y(j)+(k1+2*k2+2*k3+k4)/6;endR=T' Y'4.3.2 實例分析：例 用四階龍格庫塔方法求解初值問題取步長h=0.1。并將a,b,ya,M分別取定0,1,0,10。解：首先寫出執(zhí)行函數(shù)的tt.m文件，其中tt表示f，即為求導(dǎo)函數(shù)：function f=tt(t,y)f=t2+t-y;再輸入調(diào)用四階Runge-Kutta方法，執(zhí)行語句如下：>> R=rk4('tt',0,1,0,10)輸出：R = 0 0 0.10000000000000 0.05162708333

25、333 0.20000000000000 0.16847829348958 0.30000000000000 0.30751722078089 0.40000000000000 0.46666197888861 0.50000000000000 0.64581185656207 0.60000000000000 0.84496198189452 0.70000000000000 1.06411211920748 0.80000000000000 1.30326225710000 0.90000000000000 1.56241239502055 1.00000000000000 1.84156

26、253294245 4.4 實驗結(jié)論龍格-庫塔法具有精度高，收斂，穩(wěn)定（在一定條件下），計算過程中可以改變步長，不需要計算高階導(dǎo)數(shù)等優(yōu)點，但仍需計算在一些點上的值，如四階龍格-庫塔方法每計算一步需要計算四次的值，這給實際計算帶來一定的復(fù)雜性。實驗五 數(shù)值方法實際應(yīng)用5.1 實驗?zāi)康?）了解梯形公式原理及幾何意義；2）學會復(fù)化梯形求積公式的編程與應(yīng)用；3）掌握Matlab提供的計算積分的各種函數(shù)的使用方法；5.2 基本原理 在實際問題中，往往會遇到一些困難。有些函數(shù)找不到用初等函數(shù)表示的原函數(shù)，例如，對于積分 (5.1)而言，不存在用初等函數(shù)表示的原函數(shù)。而有些函數(shù)雖然能找到原函數(shù)，但計算過于復(fù)

27、雜，例如，橢圓型積分 (5.2) 而有些情況下，只能知道某些點處的函數(shù)值，并沒有函數(shù)的具體表達式。這些情況，使我們有必要研究積分的數(shù)值計算問題。下面我們就以梯形公式為例做以說明。 所謂梯形求積公式就是用梯形面積來近似曲邊梯形面積，利用梯形公式和連續(xù)增加a,b的區(qū)間數(shù)來逼近： (5.3)第j次循環(huán)在個等距節(jié)點處對采樣。5.3 實驗內(nèi)容 衛(wèi)星軌道是一個橢圓，橢圓周長計算公式是，這里a是橢圓半長軸，c是地球中心與軌道中心（橢圓中心）的距離，記h為近地點距離，H為遠地點距離，R=637km為地球半徑，則我國第一顆人造衛(wèi)星近地點距離h=439km，遠地點距離H=2384，試求衛(wèi)星軌道的周長。解: 第一步

28、：先利用Matlab畫出被積函數(shù)的圖形。 輸入程序如下：clearH=2384;h=439;R=6371;a=(2*R+H+h)/2c=(H-h)/2x=0:0.1:pi/2;y=sqrt(1-(c/a)2*(sin(x).2);plot(x,y,'-')title('梯形法則');xlabel('x');ylabel('y');輸出結(jié)果：a = 7.782500000000000e+003c = 9.725000000000000e+002輸出圖形：圖5.1 被積函數(shù)的圖形 第二步：應(yīng)用數(shù)值積分梯形公式。 首先建立一個名為tra

29、pezg.m的M文件，程序如下：function I=trapezg(f_name3,a,b,n)format long%輸出用15位數(shù)字表示n=n;h=(b-a)/n; x=a+(0:n)*h;f=feval(f_name3,x);I=h/2*(f(1)+f(n+1);if n>1 I=I+h*sum(f(2:n);endh1=(b-a)/100;xc=a+(0:100)*h1; fc=feval(f_name3,xc);plot(xc,fc,'r');hold onxlabel('x');ylabel('y');plot(x,f)title('數(shù)值積分梯形效果圖');plot(x,zeros(size(x),'.')for i=1:n;plot(x(i),x(i),0,f(i),end 然后建立一個名為f_name3.m的M文件定義函數(shù)，Matlab命令如下：function y=f_name3(x)y=sqrt(1-(