2014數(shù)值分析MATLAB上機實驗

1、數(shù)值分析實習報告姓名：gestepoA學號：201*班級：* 班序言隨著計算機技術(shù)的迅速發(fā)展， 數(shù)值分析在工程技術(shù)領(lǐng)域中的應用越來越廣泛， 并且成為數(shù)學與計算機之間的橋梁。 要解決工程問題， 往往需要處理很多數(shù)學模 型，不僅要研究各種數(shù)學問題的數(shù)值解法， 同時也要分析所用的數(shù)值解法在理論 上的合理性，如解法所產(chǎn)生的誤差能否滿足精度要求： 解法是否穩(wěn)定、 是否收斂 及熟練的速度等。而且還能減少大量的人工計算。由于工程實際中所遇到的數(shù)學模型求解過程迭代次數(shù)很多，計算量很大，所 以需要借助如MATLAB C+, VB, JAVA的輔助軟件來解決，得到一個滿足誤差限 的解。本文所計算題目，均采用 MA

2、TLABS行編程，MATLA被稱為第四代計算機 語言，利用其豐富的函數(shù)資源，使編程人員從繁瑣的程序代碼中解放出來MATLAB最突出的特點就是簡潔，它用更直觀的、符合人們思維習慣的代碼。 它具有以下優(yōu)點：1友好的工作平臺和編程環(huán)境。MATLAB面精致，人機交互性強，操作簡單。2簡單易用的程序語言。MATLAB!個高級的矩陣/陣列語言，包含控制語 言、函數(shù)、數(shù)據(jù)結(jié)構(gòu)，具有輸入、輸出和面向?qū)ο缶幊烫攸c。用戶可以在命令窗 口中將輸入語句與執(zhí)行命令同步，也可以先編好一個較大的復雜的應用程序（ M 文件）后再一起運行。3 強大的科學計算機數(shù)據(jù)處理能力。包含大量計算算法的集合，擁有 600 多 個工程中要用

3、到的數(shù)學運算函數(shù)。4 出色的圖像處理功能， 可以方便地輸出二維圖像， 便于我們繪制函數(shù)圖像。目錄1 第一題 41.1 實驗目的 41.2 實驗原理和方法 41.3 實驗結(jié)果 51.3.1 最佳平方逼近法 51.3.2 拉格朗日插值法 71.3.3 對比 82 第二題 92.1 實驗目的. 92.2 實驗原理和方法 102.3 實驗結(jié)果 102.3.1 第一問. 102.3.2 第二問. 112.3.3 第三問 113 第三題 123.1 實驗目的. 123.2 實驗原理和方法 123.3 實驗結(jié)果 124 MATLAB?序 141第一題某過程涉及兩變量x和y,擬分別用插值多項式和多項式擬合給出

4、其對應規(guī)律的 近似多項式，已知xi與yi之間的對應數(shù)據(jù)如下：12345678910x12345678910y34.640.314.6-14.2-13.324.875.2103.597.478.2588719448721570234795743847392請用次數(shù)分別為3, 4, 5, 6的多項式擬合并給出最好近似結(jié)果f(x) 請用插值多項式給出最好近似結(jié)果。1.1實驗目的：在某個范圍內(nèi)近似代學習逼近和插值的原理和編程方法，由給出的已知點構(gòu)造多項式,替已知點所代表的函數(shù)，以便于簡化對未知函數(shù)的各種計算。1.2試驗原理和方法:實驗原理：拉格朗日插值法中先構(gòu)造插值基礎(chǔ)函數(shù):?= n?=o?初?加（?

5、= 0,12? , ?，然后構(gòu) j *k?造出拉格朗日多項式：?（?=空=。（n?=0? ?。j r ? ?最佳平方逼近中，設(shè)逼近函數(shù)？?？ ?=? + ?+ ? + ? + ?,逼近函數(shù)和真實函? 1?數(shù)之差 r = ? - ? 1 ? = ?蔦-?,即:?= ? ?根據(jù)最小二?1?乘準則令En?=o? = ?,?可以得到？?=仔的?-1 ?實驗方法：逼近法采用最佳平方逼近，依據(jù)最小二乘原則：馬?=0? = ?,?油已知條件采用離散型。插值法采用拉格朗日插值法。在逼近法中，由于是離散型的，所以法方程系數(shù)陣設(shè)計成求和。分別求出3、4、5、6次的多項式，逼近結(jié)果和真實值有一定差距， 最小二乘正是

6、讓這些差距達到最小，理論上多項式次數(shù)越高結(jié)果和真實值差距越小。拉格朗日插值法中“l(fā)a=la*(p-x(j)/(x(k)-x(j)”語句實現(xiàn)的是我們通常書寫的連乘形式拉格朗日插值多項式，但是表示不方便，而如果用“ s=collect(s)”函數(shù)將其展開成降幕排列多項式以后，由于余項問題結(jié)果會和原本的多項式有偏差，這種偏差隨著x的增大而增大。求出多項式后和題目中給出的參考點進行比較。最后，選擇六次最佳平方逼近多項式和拉格朗日插值多項式(九次)進行比較，選取xi=a+ih=1+0.2*i(i=0,1, ? ,45)，分別繪制兩者的圖像進行比較。1.3試驗結(jié)果1.3.1最佳平方逼近法三次多項式：-1.

7、033*xA3 + 19.33*xA2 - 94.48*x + 131.8擬合結(jié)果:12345678910x12345678910y55.61711.896-5.561-2.95213.52537.67263.29184.18494.15387.0000000000000四次多項式：-0.3818*xA4 + 7.368*xA3 - 42.14*乂人2 + 73.53*x+ 0.745 擬合結(jié)果：12345678910x12345678910y39.12132.08010.085-5.563-2.73021.56061.117100.588115.45772.045222802222017 I

8、10D83B-l102?五次多項式：0.09807*L5 - 3.079P4 + 34.5*tA3 -163.5*2 + 304.7*t -139.5 擬合結(jié)果：12345678910x12345678910y33.21945.7749.032-16.500-8.90626.90870.983100.07399.41674.5001203335840六次多項式：0.01936*tA6 - 0.5408*tA5 + 5.114十4 - 16.93 - 0.8672 + 66.38*t - 18.7 擬合結(jié)果：12345678910x12345678910y34.50541.14913.270-1

9、3.94-12.2223.71173.950105.1696.34578.4006408650409460對比可知，六次多項式擬合結(jié)果最好。1.3.2拉格朗日插值法插值多項式 5.353*10A( -5)*xA9 - 0.003088*xA8 + 0.07229*xA7 - 0.8792*xA6 + 5.932*xA5-22.41*xA4 + 50.1計3 - 86.47*乂人2 + 113.5*x - 25.2注：此多項式為拉格朗日多項式的近似式，當x=10的時候偏差可以達到 23以上。對比數(shù)據(jù)：123456789x1.50001.90002.30002.70003.10003.50003.

10、90004.30004.7000y42.149841.46235.11824.38511.273-1.7813-12.300-18.1566-17.9060222691011121314151617x5.10005.50005.90006.30006.70007.10007.50007.9000y-11.0222.028419.85440.36261.08479.56893.7700102.367696087插值結(jié)果:123456789x1.50001.90002.30002.70003.10003.50003.90004.30004.7000y42.384041.49435.07424.36

11、011.279-1.7683-12.297-18.1626-17.9117212781011121314151617x5.10005.50005.90006.30006.70007.10007.50007.9000y-11.0212.033319.85640.35861.07979.57093.7788102.371054493其中紅點表示參考點。1.3.3比較選取xi=a+ih=1+0.2*i(i=0,1, ? ,45)，分別繪制六次多項式擬合和拉格朗日插值結(jié)果圖:120其中綠線表示拉格朗日插值多項式圖像，藍線表示六次多項式擬合圖像。兩者效果近似但后者比前者低三次。2第二題用雅格比法與高斯-

12、賽德爾迭代法解下列方程組 Ax=b1或Ax=b2，研究其收斂性。 上機驗證理論分析是否正確，比較它們的收斂速度，觀察右端項對迭代收斂有無 影響。(1) A行分別為 A仁6,2, - 1, A2=1,4, - 2, A3=- 3,1,4 ;b1=- 3,2,4 T;b2=100,-200,345 T。(2) A行分別為 A1=1,0,8,0.8,A2=0.8,1,0.8,A3=0.8,0.8,1; b1=3,2,1T;b2=5,0, - 10T。 A行分別為 A仁1,3, A2=-7,1 ; b1=4,6 T。2.1試驗目的學習jacobi迭代法和GuassSeidel迭代法的原理和編程方法，研

13、究方程組系數(shù)陣和右邊 項對方程的解及其收斂性的影響，判斷迭代法的收斂條件。2.2 實驗原理和方法實驗原理：將方程組系數(shù)陣 A分解為A= D+ L+ U,其中D為對角陣，L為減去D的下三角陣，U 為減去 D 的上三角陣。Jacobi 迭代法中構(gòu)造如下迭代公式：?( ?+1) = -? -1 (?+ ?) ?( ?) + ?-1 ? 而 Gauss-Seidel 迭代法的迭代公式為：?(?+1) = - (?+ ?)-1 ?(?) + (?+ ?)-1 ?/ f初始值直接選取為0。在判斷其收斂性時，分別求解其迭代矩陣的譜半徑p?, P? =/?為迭代矩陣的特征值。實驗方法：分別編寫 jacobi

14、迭代及其收斂判別函數(shù)和 Seidel 迭代及其收斂判別函數(shù)。 如果在初試迭 代步數(shù)之內(nèi)還未收斂就進行收斂判別，收斂判別的依據(jù)是迭代矩陣的譜半徑是否小于1。比較同一方程組的 jacobi 迭代法和 Seidel 迭代法的結(jié)果是否相同， 在達到精度要求后比 較兩種方法的迭代次數(shù)，比較哪一個的效率更高。比較方程組系數(shù)陣和等號右邊的變化會對方程的解和收斂速度造成什么影響。 如果迭代不收斂，那么考慮為什么不收斂，如果把方程組系數(shù)陣進行強對角占優(yōu)處理， 是否會收斂。2.3 實驗結(jié)果規(guī)定誤差界： 1e-42.3.1 第一問62-1-3 ?= 14-2 ,?= 2 .-3144-0.7273由 jacobi

15、迭代法求得?= 0.8081，設(shè)定迭代 20 次，實際迭代 16 次，精度為 9.4022e-005。0.2525-0.7272由 seidel 迭代法求得 ?= 0.8081 ，設(shè)定迭代 20 次，實際迭代 10 次，精度為：0.25259.0769e-00562-1100 ?= 14-2 ，?= -200 -31434536.3636由 jacobi 迭代法求得 ?= -2.0707 ，設(shè)定迭代 40 次，實際迭代 23 次，精度為114.04046.9948e-005。36.3637由 Seidel 迭代法求得 ?= -2.0707 ，設(shè)定迭代 20 次，實際迭代 15 次，精度為114

16、.04048.6384e-005通過對比可知： 1、 Seidel 迭代的收斂速度明顯高于 jacobi 迭代。 2、 b 矩陣對收斂速 度和誤差精度有影響， b 中元素較大時會放慢收斂速度并加大誤差。2.3.2 第二問1 ?= 0.80.80.8 0.831 0.8，?= 20.811由 jacobi 迭代法求解， 100 次迭代尚且不能達到精度。此時調(diào)用 jacobi 迭代法的收斂判別函數(shù)，求得特征值為：1 = -1.6 ,2、3 = 0.8 , P (?= 1.6 1，迭代不收斂。5.7691由于 Seidel 迭代法求解 ?= 0.7693 ，迭代次數(shù) 31，精度為 8.7826e-0

17、05-4.230710.80.85 ?= 0.810.8，?= 0 0.80.81-10Jacobi 迭代不收斂。32.6922Seldel 迭代法求得 ?= 7.6922 ，迭代次數(shù) 38，精度為 8.4552e-005。-42.3076比較得知 A矩陣元素如果相差很小，迭代次數(shù)會大幅增加，綜合比較可知b矩陣元素如果相差很大會增加迭代次數(shù)。2.3.3 第三問試驗結(jié)果和討論?= 1-731，?= 4?= 6此時 p （第=4.5826 1 ,p(?3= 21 1 , jacobi迭代法和Seidel迭代法都不收斂。如果交換 A 中行的順序， 得到 -71；，用jacobi迭代計算，迭代8次，解

18、得x = I；：64 。用Seidel迭代法計算，只需迭代5次，得到x = -0.6364 ，精度為2.6326e-005。此時p (?= 1.54550.2182 p(?) = 0.0476 ，從此可以看出收斂速度的快慢。3第三題1給定函數(shù)f(x)r，5x5,及節(jié)點Xi 5 i, i 0,1,L 10,求其三次樣1 x2條插值多項式(可取I型或II型邊界條件)，并畫圖及與f(x)的圖形進行比較分析。注：涉及到線性方程組求解問題需采用適當?shù)那蠼馑惴ā?.1實驗目的學習三彎矩法的原理和編程方法，對比原函數(shù)和三次樣條插值的結(jié)果。3.2實驗原理和方法實驗原理:記 h?= ? ?_1， ?尸?+1CC

19、6?+1? ?_1，人丄十 /+?= 1 - ? ?=廠（一 -?），給出插值?+? ?+1?+1re? (?) = ?) ?(? = ?(? = ?。組成遞推式：?_1 + 2?+ ?+1= ?由于系數(shù)陣按行對角占優(yōu)，方程組存在唯一確定解，可以使用高斯列主元消去法來解 方程。最后將各個參數(shù)帶入樣條函數(shù)(?_ ?3(? ?_1)3?_1 2?= ?_1 二 + ? + (?菊1_廠？?)6? ? 6? ? 6即可求得樣條函數(shù)。兩端處的二階導數(shù)?(?)=? ? 2-7+ （?毎-?）?6?- ?_1?實驗方法由于在本題中x(i+1)_x(i)=1,所以h(i)=1。在編程中直接將 h設(shè)置成常數(shù)，

20、簡化了運算。首先求解、g，然后列出方陣求解 M(i)，在求解方程組的過程中采用列主元素高斯消去法，分為消元和回代兩個過程，編寫這兩個函數(shù)，解出除了兩端的M，而兩端點的 M值等于兩端點的函數(shù)二階導數(shù)值。編寫函數(shù)求出樣條函數(shù)的系數(shù)，然后求出方程，對于三彎矩法三次樣條函數(shù)，如果有n個點，則有n_1個樣條函數(shù)，除了兩端需要求解n_2個M值，即解n_2階方程組。在表達樣條函數(shù)的時候采用if語句，對不同的區(qū)間進行劃分， 然后細分_5,5這個區(qū)間，間隔0.1將其分為100份，這樣可以體現(xiàn)出連續(xù)性，此時繪圖對比三次樣條函數(shù)和原函數(shù)。本題中原函數(shù)的二階導不是計算機解出的沒有編寫相關(guān)程序。本題中求解的樣條函數(shù)，M

21、ATLAB系統(tǒng)自動的將公因式提取，并且合并同類項，所以表達出的函數(shù)并不整齊規(guī)律，為了更好地體現(xiàn)三次樣條函數(shù)的結(jié)構(gòu)和性質(zhì)，我專門手寫了規(guī)整的樣條函數(shù)。3.3實驗結(jié)果三彎矩法求解代入 x=-5-4-3-2-10123451 y = f(?=21 + ?0.0588? ( - 5 ) = ? (5)0.00842y=0.03850.05880.10000.20000.50001.00000.50000.20000.10000.0385S=M=0.01410.05990.09930.7431-1.87150.74310.09930.05990.01410.00140.00240.03710.05650

22、.00240.01000.05650.09000.01000.01650.09000.18350.01650.12380.18350.37620.1238-0.31190.37621.3119-0.31190.12381.31190.37620.12380.01650.37620.18350.01650.01000.18350.09000.01000.00240.09000.05650.00240.00140.05650.0371求得樣條函數(shù)系數(shù)陣為：樣條函數(shù)：(97*X)/5000 - (7*(X + 4)A3)/5000 + (3*(X + 5)A3)/1250 + 1341/10000(

23、67*X)/2000 - (3*(X + 3)A3)/1250 + (X + 4)人3/100 + 381/2000 (187*X)/2000 - (X + 2)A3/100 + (33*(X + 3)人3)/2000 + 741/2000 (1927*X)/10000 - (33*(X + 1)人3)/2000 + (619*(X + 2)人3)/5000 + 5689/10000 (9357*X)/10000 - (3119*(X + 1)人3)/10000 - (619*乂人3)/5000 + 13119/10000 (3199*(X - 1)A3)/10000 - (9357*X)/1

24、0000 + (619*乂人3)/5000 + 13119/10000 (33*(X - 1)A3)/2000 - (1927*X)/10000 - (619*(X - 2)人3)/5000 + 5689/10000 (X - 2)A3/100 - (187*X)/2000 - (33*(X - 3)人3)/2000 + 741/2000 (3*(X - 3)A3)/1250 - (67*X)/2000 - (X - 4)人3/100 + 381/2000 (7*(X - 4)A3)/5000 - (97*X)/5000 - (3*(X - 5)人3)/1250 + 1341/10000 寫出

25、標準形式的樣條函數(shù)：-4-x)+0.0565*(x+5)-3-x)+0.0900*(x+4)-2-x)+0.1835*(x+3)-1 -x)+0.3762*(x+2)x-5,-4x-4,-3x-3,-2x-2,-1 x-1,0 x0,1S1(x)= 0.0014*( -4-x)A3+0.0024*(x+5)A3+0.0371*(S2(x)= 0.0024*( -3-x)A3+0.0100*(x+4)A3+0.0565*(S3(x)= 0.0100*( -2-x)A3+0.0165*(x+3)A3+0.0900*(S4(x)= 0.0165*( -1 -x)A3+0.1238*(x+2)A3+0

26、.1835*(S5(x)= 0.1238*( -x)A3-0.3119*(x+1)A3+0.3762*( -x)+1.3119(x+1)S6(x)= -0.3119*(1 -x)A3+0.1238*xA3+1.3119*(1 -x)+0.3762*xS7(x)= 0.1238*(2-x)A3+0.0165*(x-1)人3+0.3762*(2-x)+0.1835*(x-1)x1,2S8(x)= 0.0165*(3-x)A3+0.0100*(x-2)人3+0.1835*(3-x)+0.0900*(x-2)x2,3S9(x)= 0.0100*(4-x)A3+0.0024*(x-3)人3+0.0900

27、*(4-x)+0.0565*(x-3)x3,4S10(x)= 0.0024*(5 -x)A3+0.0014*(x -4)人3+0.0565*(5 -x)+0.0371*(x -4)x4,5輸入：a=-5:0.1:5;for i=1:le ngth(a) b(i)=f(a(i);end求得原函數(shù)和樣條函數(shù)的對比圖:4 MATLAB 程序第一題：離散型最佳平方逼近函數(shù)fun cti onf=squar_approx_ls(x, y,n)syms t ;N=le ngth(x);P=zeros( n+1);Q=zeros( n+1,1);a=0;for i=0: nfor j=0:nfor k=1:Na=a+x(kF(i+j);endP(i+1,j+1)=a;a=0;endendb=0;for i=0: nfor k=1:Nb=b+y(k)*x(k)AiendQ(i+1,1)=b;b=0; end s=inv(P)*Q; f=0; for i=1:n+1f=f+s(i)*tq-1);simplify(f); end f=collect(f); f=vpa(f,4);拉格朗日插值函數(shù) function s=Lagrange(x,y,x0) syms p;n=length(x);s=0;for (k=1: