數(shù)值分析大作業(yè)牛頓下山法拉格朗日法切比雪夫法及Matlab程序

1、課 程 設(shè) 計(jì)課程名稱： 數(shù)值分析 設(shè)計(jì)題目： 學(xué) 號(hào)： 姓 名： 完成時(shí)間： 題目一： 解線性方程組的直接法設(shè)方程組，其中，矩陣中，由相應(yīng)的矩陣元素計(jì)算，使解向量。(1) 不變，對(duì)的元素加一個(gè)擾動(dòng)，求解方程組；(2) 不變，對(duì)的元素和分別加一個(gè)擾動(dòng)，求解方程組；(3) 對(duì)上述兩種擾動(dòng)方程組的解做誤差分析。一.數(shù)學(xué)原理： 本計(jì)算采用直接法中的列主元高斯消元法，高斯列主元消元法原理如下：1、設(shè)有n元線性方程組如下：2、第一步：如果a11!=0, 令li1= ai1/a11, I= 2,3,n用（-li1）乘第一個(gè)方程加到第i個(gè)方程上,得同解方程組：a(1)11 a(1)12 . . . a(1)

2、1n x1 b(1)1a(1)21 a(1)22 . . . a(1)2n x2 b(1)2 . . . . . . . = .a(1)n-11 a(1)n-12 . . a(1)n-1n xn-1 b(1)n-1a(1)n1 a(1)n2 . . . a(1)nn xn b(1)n簡(jiǎn)記為： A(2) x = b(2)其中 a(2)ij = a(1)ij li1 * a(1)1j , I ,j = 2,3,.,n b(2)I = b(1)I li1 * b(1)1 , I = 2,3,.,n第二步：如果a(2)22 != 0,令li2= a（2）i2/a（2）22, I= 3,n依據(jù)同樣的原理

3、，對(duì)矩陣進(jìn)行化間（省略），依次下去，直到完成！最后，得到上三角方程組： a(1)11 a(1)12 . . . a(1)1n x1 b(1)10 a(1)22 . . . a(1)2n x2 b(1)2 . . . . . . . = .0 0 . . a(n-1)n-1n xn-1 b(n-1)n-10 0 . . . a(n)nn xn b(n)n簡(jiǎn)記為： A(n) x = b(n)最后從方程組的最后一個(gè)方程進(jìn)行回代求解為： Xn = b(n) / a(n)nn Xi = ( b(k)k - a(k)kjxj ) / a(k)kk二解題過程：1.由題中所給條件可求出b。B = 6.0000

4、 7.7156 9.9299 12.7560 16.3238 20.7813（1） 不變，對(duì)的元素加一個(gè)擾動(dòng)，求解方程組。B =6.0000 7.7156 9.9299 12.7560 16.3238 20.7813+0.0001解得x =0.5997 2.6920 -1.8500 3.3917 0.0000 1.1667（2）不變，對(duì)的元素和分別加一個(gè)的擾動(dòng)，求解方程組。 A = 1.000000000000000 1.000000000000000 1.000000000000000 1.000000000000000 1.000000000000000 1.000000000000000

5、 1.000000000000000 1.500000000000000 2.250000000000000 3.375000000000000 5.062500000000000 7.593751000000000x =三、誤差分析：從上面計(jì)算結(jié)果可以看出，當(dāng)系數(shù)矩陣或右端向量發(fā)生極小的擾動(dòng)，方程組的解也會(huì)產(chǎn)生很大的誤差，產(chǎn)生的原因是范德蒙陣為變態(tài)陣。由數(shù)值計(jì)算知識(shí)可知 其中 為條件數(shù)，從上式看到，當(dāng)A的條件數(shù)很大時(shí)，解的相對(duì)誤差也很大，此時(shí)的對(duì)應(yīng)的線性方程為病態(tài)線性方程組。計(jì)算條件數(shù)時(shí)，取矩陣的無窮范數(shù)，經(jīng)計(jì)算得矩陣、受擾動(dòng)后的矩陣和的條件數(shù)為；可以看到三個(gè)矩陣的條件數(shù)非常大，即使系數(shù)矩陣

6、或右端向量發(fā)生很小的變化，也會(huì)導(dǎo)致解產(chǎn)生很大的誤差。四、收獲與體會(huì)：運(yùn)用matlab編程解決數(shù)學(xué)問題很方便，病態(tài)陣的條件數(shù)非常大，給系數(shù)矩陣或者右端向量一個(gè)很微小的擾動(dòng)，方程組的解也會(huì)產(chǎn)生很大的變動(dòng)。通過做這個(gè)題目，我對(duì)讓課上抽象病態(tài)現(xiàn)象有了直觀的認(rèn)識(shí)。題目二：多項(xiàng)式插值在區(qū)間上對(duì)龍格函數(shù)做插值，分別用等距節(jié)點(diǎn)(節(jié)點(diǎn)步長(zhǎng))和切比雪夫多項(xiàng)式的零點(diǎn)做插值節(jié)點(diǎn)，畫出原函數(shù)和兩個(gè)插值函數(shù)的圖像進(jìn)行比較，并利用對(duì)兩種插值方法做誤差分析。一.數(shù)學(xué)原理：1.拉格朗日等距節(jié)點(diǎn)插值拉格朗日插值多項(xiàng)式為 其中 余項(xiàng)為 與有關(guān)。2.拉格朗日以切比雪夫零點(diǎn)為節(jié)點(diǎn)插值切比雪夫多項(xiàng)式 在區(qū)間 上有 個(gè)零點(diǎn)，為 由上式得到

7、切比雪夫多項(xiàng)式的個(gè)零點(diǎn)，并作為插值節(jié)點(diǎn)做拉格朗日插值，得到多項(xiàng)式。一.計(jì)算過程：原函數(shù)為y=1./(1+x.2)，通過matlab繪圖得到原函數(shù)圖像為：1. 利用拉格朗日差分(等距節(jié)點(diǎn))：x =(-5 -4 -3 -2 -1 0 1 2 3 4 5)y=1./(1+x.2)= (0.0385 0.0588 0.1000 0.2000 0.5000 1.0000 0.5000 0.2000 0.1000 0.0588 0.0385)拉格朗日插值多項(xiàng)式為：=(t*(t/8 + 5/8)*(t - 1)*(t - 2)*(t + 2)*(t - 3)*(t + 3)*(t - 4)*(t + 4)*

8、(t - 5)/4320 - (t*(t/15 + 1/3)*(t - 1)*(t + 1)*(t - 2)*(t - 3)*(t + 3)*(t - 4)*(t + 4)*(t - 5)/10080 - (t/5 + 1)*(t - 1)*(t + 1)*(t - 2)*(t + 2)*(t - 3)*(t + 3)*(t - 4)*(t + 4)*(t - 5)/2880 + (t*(t/20 + 1/4)*(t - 1)*(t + 1)*(t - 2)*(t + 2)*(t - 3)*(t - 4)*(t + 4)*(t - 5)/40320 + (t*(t/12 + 5/12)*(t

9、+ 1)*(t - 2)*(t + 2)*(t - 3)*(t + 3)*(t - 4)*(t + 4)*(t - 5)/2880 - (t*(t/17 + 5/17)*(t - 1)*(t + 1)*(t - 2)*(t + 2)*(t - 3)*(t + 3)*(t - 4)*(t - 5)/362880 + (t*(t/26 + 2/13)*(t - 1)*(t + 1)*(t - 2)*(t + 2)*(t - 3)*(t + 3)*(t - 4)*(t - 5)/3628800 - (t*(t/35 + 1/7)*(t - 1)*(t + 1)*(t + 2)*(t - 3)*(t

10、+ 3)*(t - 4)*(t + 4)*(t - 5)/4320 + (t*(t/80 + 1/16)*(t - 1)*(t + 1)*(t - 2)*(t + 2)*(t + 3)*(t - 4)*(t + 4)*(t - 5)/10080 - (t*(t/153 + 5/153)*(t - 1)*(t + 1)*(t - 2)*(t + 2)*(t - 3)*(t + 3)*(t + 4)*(t - 5)/40320 + (t*(t/260 + 1/52)*(t - 1)*(t + 1)*(t - 2)*(t + 2)*(t - 3)*(t + 3)*(t - 4)*(t + 4)/36

11、2880注：Matlab中用t表示變量x通過matlab繪圖的插值后多項(xiàng)式曲線為：2切比雪夫零點(diǎn)插值：在-1,1上有11個(gè)不同零點(diǎn)： ，=（0.9898 0.9096 0.7557 0.5406 0.2817 0.0000 -0.2817 -0.5406 -0.7557 -0.9096 -0.9898）將-1,1區(qū)間到-5,5區(qū)間進(jìn)行轉(zhuǎn)化得：x=5*=（4.9491 4.5482 3.7787 2.7032 1.4087 0.0000 -1.4087 -2.7032 -3.7787 -4.5482 -4.9491）y=1./(1+x.2)=（ 0.0392 0.0461 0.0654 0.12

12、04 0.3351 1.0000 0.3351 0.1204 0.0654 0.0461 0.0392）由于插值后的多項(xiàng)式過于復(fù)雜，這里沒有給出具體形式，直接給出函數(shù)曲線為：原函數(shù)、拉格朗日插值、切比雪夫零點(diǎn)插值后的曲線對(duì)比如下：3.誤差分析：從上圖可以看出，在區(qū)間兩端高次等間距的朗格朗日插值與原函數(shù)之間的偏差很大，切比雪夫零點(diǎn)插值多項(xiàng)式在整個(gè)區(qū)間上與原函數(shù)之間的偏差都很小，所以直觀上可以看出切比雪夫插值多項(xiàng)式更趨近原函數(shù)。下面利用無窮范數(shù)對(duì)兩種插值方法做誤差分析，計(jì)算結(jié)果如下：對(duì)比上面計(jì)算結(jié)果易知用切比雪夫多項(xiàng)式的零點(diǎn)做插值節(jié)點(diǎn)得到的多項(xiàng)式于原函數(shù)的擬合度更高。4.收獲與體會(huì)：利用matla

13、b的強(qiáng)大數(shù)學(xué)計(jì)算功能，加深了對(duì)拉格朗日插值方法的認(rèn)識(shí)。深刻體會(huì)到對(duì)于一些函數(shù)，等間距節(jié)點(diǎn)的高次插值多項(xiàng)式的偏差會(huì)很大，通過計(jì)算可以看到，利用相應(yīng)的切比雪夫多項(xiàng)式零點(diǎn)做拉格朗日插值會(huì)消弱這一“龍格”現(xiàn)象，而且效果非常明顯。題目三：非線性方程求根利用改進(jìn)Newton法求解方程，其中迭代條件分別為：(1) (2) (3) 對(duì)不同條件下獲得的結(jié)果進(jìn)行分析比較。一.數(shù)學(xué)原理：牛頓法解非線性方程的近似解，已知的近似解，通過下式得到更精確地零點(diǎn)近似根。 該題的迭代方式采用的牛頓下山法，是基于牛頓法的改進(jìn)，即在下山法保證函數(shù)值穩(wěn)定下降的前提下，用牛頓法加快收斂速度，得出下面的迭代計(jì)算公式 其中稱為下山因子。由

14、數(shù)值計(jì)算知識(shí)可知，雖然(線性方程的重根數(shù))大于1但上式對(duì)也是二階收斂的，所以計(jì)算時(shí)的初始值取，迭代終止條件為。二.計(jì)算過程：（1），；程序運(yùn)行結(jié)果如下：x=0.50175381y=0.50175381 k=10001（k為迭代次數(shù)）由于設(shè)定了迭代次數(shù)最大為10001，因此迭代10001次時(shí)還沒有使，。（2）程序運(yùn)行結(jié)果如下：x=0.50000000y=0.50000000 k=25可見迭代至25次即找到了最優(yōu)解。(3) x=0.99999999y=0.99999999 k=12可見迭代至12次即找到了最優(yōu)解。三結(jié)果分析：由方程式可知其根為，。(1)、(2)中的初始值靠近，迭代時(shí)向收斂，但（2）

15、中r=1.5可知函數(shù)在(2)初始條件下得到的解更精確； (2)迭代25次滿足精度要求，可知函數(shù)在(2)初始條件下收斂的更快。(3)的初始值為0.85靠近，迭代時(shí)向收斂，迭代12次滿足精度要求，對(duì)比可知在(3)初始條件下，函數(shù)的收斂速度最快，得到的解最精確。四收獲與體會(huì)：牛頓下山法，不僅收斂速度快而且精度高。在所給三個(gè)不同初始條件下，迭代的次數(shù)和近似解的精確度均不一樣，初始重根數(shù)越接近1，收斂速度更快。若非線性方程有不同根，給的初始零點(diǎn)值不同時(shí)，函數(shù)迭代過程會(huì)趨向靠近的真實(shí)根。數(shù)值計(jì)算是與計(jì)算機(jī)密切相關(guān)的一門課程，通過書本與matlab相結(jié)合很好地理解了其中的原理與方法，為以后學(xué)習(xí)提供了強(qiáng)有力地

16、工具和方法。附錄 程序設(shè)計(jì)題目一： 解線性方程組的直接法 clear A=1 1 1 1 1 1;1 1.1 1.12 1.13 1.14 1.15;1 1.2 1.22 1.23 1.24 1.25;1 1.3 1.32 1.33 1.34 1.35;1 1.4 1.42 1.43 1.44 1.45;1 1.5 1.52 1.53 1.54 1.55A = 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.1000 1.2100 1.3310 1.4641 1.6105 1.0000 1.2000 1.4400 1.7280 2.0736

17、 2.4883 1.0000 1.3000 1.6900 2.1970 2.8561 3.7129 1.0000 1.4000 1.9600 2.7440 3.8416 5.3782 1.0000 1.5000 2.2500 3.3750 5.0625 7.5938 X=1;1;1;1;1;1X = 1 1 1 1 1 1 B=A*XB = 6.0000 7.7156 9.9299 12.7560 16.3238 20.7813 Aug=A B n,m=size(Aug)n = 6m = 7 for k = 1:n-1 piv,r = max(abs(Aug(k:n,k); %找列主元所在子矩

18、陣的行r r = r + k - 1; % 列主元所在大矩陣的行 if rk temp=Aug(k,:); Aug(k,:)=Aug(r,:); Aug(r,:)=temp; end if Aug(k,k)=0, error(對(duì)角元出現(xiàn)0), end % 把增廣矩陣消元成為上三角 for p = k+1:n Aug(p,:)=Aug(p,:)-Aug(k,:)*Aug(p,k)/Aug(k,k); end end AugAug = 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 6.0000 0 0.5000 1.2500 2.3750 4.0625 6.5

19、938 14.7813 0 0 -0.0600 -0.2220 -0.5514 -1.1492 -1.9826 0 0 0 -0.0080 -0.0408 -0.1306 -0.1794 0 0 0 0 -0.0012 -0.0074 -0.0086 0 0 0 0 0 0.0001 0.0001 % 解上三角方程組 A = Aug(:,1:n); b = Aug(:,n+1); x(n) = b(n)/A(n,n); for k = n-1:-1:1 x(k)=b(k); for p=n:-1:k+1 x(k) = x(k)-A(k,p)*x(p); end x(k)=x(k)/A(k,k)

20、; end xx = 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000%* %（1)A不變，對(duì)B的元素b6加一個(gè)擾動(dòng)10e-4，求解方程組； A=1 1 1 1 1 1;1 1.1 1.12 1.13 1.14 1.15;1 1.2 1.22 1.23 1.24 1.25;1 1.3 1.32 1.33 1.34 1.35;1 1.4 1.42 1.43 1.44 1.45;1 1.5 1.52 1.53 1.54 1.55A = 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.1000 1.2100 1.

21、3310 1.4641 1.6105 1.0000 1.2000 1.4400 1.7280 2.0736 2.4883 1.0000 1.3000 1.6900 2.1970 2.8561 3.7129 1.0000 1.4000 1.9600 2.7440 3.8416 5.3782 1.0000 1.5000 2.2500 3.3750 5.0625 7.5938 B =6.0000 7.7156 9.9299 12.7560 16.3238 20.7813+0.0001B = 6.0000 7.7156 9.9299 12.7560 16.3238 20.7814 Aug=A B n,

22、m=size(Aug)n = 6m = 7 for k = 1:n-1 piv,r = max(abs(Aug(k:n,k); %找列主元所在子矩陣的行r r = r + k - 1; % 列主元所在大矩陣的行 if rk temp=Aug(k,:); Aug(k,:)=Aug(r,:); Aug(r,:)=temp; end if Aug(k,k)=0, error(對(duì)角元出現(xiàn)0), end % 把增廣矩陣消元成為上三角 for p = k+1:n Aug(p,:)=Aug(p,:)-Aug(k,:)*Aug(p,k)/Aug(k,k); end end AugAug = 1.0000 1.

23、0000 1.0000 1.0000 1.0000 1.0000 6.0000 0 0.5000 1.2500 2.3750 4.0625 6.5938 14.7814 0 0 -0.0600 -0.2220 -0.5514 -1.1492 -1.9827 0 0 0 -0.0080 -0.0408 -0.1306 -0.1795 0 0 0 0 -0.0012 -0.0074 -0.0087 0 0 0 0 0 0.0001 0.0001 % 解上三角方程組 A = Aug(:,1:n); b = Aug(:,n+1); x(n) = b(n)/A(n,n); for k = n-1:-1:

24、1 x(k)=b(k); for p=n:-1:k+1 x(k) = x(k)-A(k,p)*x(p); end x(k)=x(k)/A(k,k); end xx = 0.5997 2.6920 -1.8500 3.3917 0.0000 1.1667%* %(2)B不變，對(duì)A的元素a22和a66分別加一個(gè)擾動(dòng)10e-6，求解方程組； clear format long A=1 1 1 1 1 1;1 1.1+0.000001 1.12 1.13 1.14 1.15;1 1.2 1.22 1.23 1.24 1.25;1 1.3 1.32 1.33 1.34 1.35;1 1.4 1.42 1

25、.43 1.44 1.45;1 1.5 1.52 1.53 1.54 1.55+0.000001 B =6.000000 7.715600 9.929900 12.756000 16.323800 20.781300B = 6.000000000000000 9.929900000000000 12.756000000000000 16.323799999999999 Aug=A B n,m=size(Aug); for k = 1:n-1 piv,r = max(abs(Aug(k:n,k); %找列主元所在子矩陣的行r r = r + k - 1; % 列主元所在大矩陣的行 if rk t

26、emp=Aug(k,:); Aug(k,:)=Aug(r,:); Aug(r,:)=temp; end if Aug(k,k)=0, error(對(duì)角元出現(xiàn)0), end % 把增廣矩陣消元成為上三角 for p = k+1:n Aug(p,:)=Aug(p,:)-Aug(k,:)*Aug(p,k)/Aug(k,k); end % 解上三角方程組 A = Aug(:,1:n); b = Aug(:,n+1); x(n) = b(n)/A(n,n); for k = n-1:-1:1 x(k)=b(k); for p=n:-1:k+1 x(k) = x(k)-A(k,p)*x(p); end x

27、(k)=x(k)/A(k,k); end xx =題目二：多項(xiàng)式插值在區(qū)間上對(duì)龍格函數(shù)做插值，分別用等距節(jié)點(diǎn)(節(jié)點(diǎn)步長(zhǎng))和切比雪夫多項(xiàng)式的零點(diǎn)做插值節(jié)點(diǎn)，畫出原函數(shù)和兩個(gè)插值函數(shù)的圖像進(jìn)行比較，并利用對(duì)兩種插值方法做誤差分析。原函數(shù)： x=-5:0.1:5; y=1./(1+x.2); plot(x,y)2. 利用拉格朗日差分(等距節(jié)點(diǎn))： x=-5:1:5x = -5 -4 -3 -2 -1 0 1 2 3 4 5 y=1./(1+x.2)y =0.0385 0.0588 0.1000 0.2000 0.5000 1.0000 0.5000 0.2000 0.1000 0.0588 0.03

28、85 syms t l;if(length(x) = length(y) n = length(x);else disp(x和y的維數(shù)不相等！); return; %檢錯(cuò)end p=sym(0);for (i=1:n) l=sym(y(i); for(k=1:i-1) l=l*(t-x(k)/(x(i)-x(k); end; for(k=i+1:n) l=l*(t-x(k)/(x(i)-x(k); end; p=p+l;end pp =(t*(t/8 + 5/8)*(t - 1)*(t - 2)*(t + 2)*(t - 3)*(t + 3)*(t - 4)*(t + 4)*(t - 5)/4

29、320 - (t*(t/15 + 1/3)*(t - 1)*(t + 1)*(t - 2)*(t - 3)*(t + 3)*(t - 4)*(t + 4)*(t - 5)/10080 - (t/5 + 1)*(t - 1)*(t + 1)*(t - 2)*(t + 2)*(t - 3)*(t + 3)*(t - 4)*(t + 4)*(t - 5)/2880 + (t*(t/20 + 1/4)*(t - 1)*(t + 1)*(t - 2)*(t + 2)*(t - 3)*(t - 4)*(t + 4)*(t - 5)/40320 + (t*(t/12 + 5/12)*(t + 1)*(t -

30、 2)*(t + 2)*(t - 3)*(t + 3)*(t - 4)*(t + 4)*(t - 5)/2880 - (t*(t/17 + 5/17)*(t - 1)*(t + 1)*(t - 2)*(t + 2)*(t - 3)*(t + 3)*(t - 4)*(t - 5)/362880 + (t*(t/26 + 2/13)*(t - 1)*(t + 1)*(t - 2)*(t + 2)*(t - 3)*(t + 3)*(t - 4)*(t - 5)/3628800 - (t*(t/35 + 1/7)*(t - 1)*(t + 1)*(t + 2)*(t - 3)*(t + 3)*(t -

31、 4)*(t + 4)*(t - 5)/4320 + (t*(t/80 + 1/16)*(t - 1)*(t + 1)*(t - 2)*(t + 2)*(t + 3)*(t - 4)*(t + 4)*(t - 5)/10080 - (t*(t/153 + 5/153)*(t - 1)*(t + 1)*(t - 2)*(t + 2)*(t - 3)*(t + 3)*(t + 4)*(t - 5)/40320 + (t*(t/260 + 1/52)*(t - 1)*(t + 1)*(t - 2)*(t + 2)*(t - 3)*(t + 3)*(t - 4)*(t + 4)/362880 x0=-5:0.1:5; f = subs (p,t,x0); %計(jì)算插值點(diǎn)的函數(shù)值 plot(x0,f)2切比雪夫零點(diǎn)插值：在-1,1上有11個(gè)不同零點(diǎn)： ， k=1:11k = 1 2 3 4 5 6 7 8 9 10 11 s=cos(2.*k-1).*pi./22)s =0.9898 0.9096 0