1、數值線性代數實驗題目數值線性代數專業：信息與計算科學班級：班姓名:山東科技大學2013年 1月 16日實驗報告說明學院：信息學院 專業：信息 班級10-2姓名:一、主要參考資料：（1）Matlab數值計算-案例分析北京航空出版（2）Matlab 數值分析機械工業出版二、課程設計應解決的主要問題：（1） （2）QR方法（3）最小二乘法三、應用軟件：(1)Matlab7.0（2）數學公式編輯器四、發出日期： 課程設計完成日期： 指導教師簽字： 系主任簽字： 指導教師對課程設計的評語指導教師簽字：年月日、問題描述先用你所熟悉的計算機語言將平方根和改進的平方根法編成寫通用的子程序，然后用你編寫的程序求

2、解對稱正定方程組.-.X = b，其中(1)b隨機的選取，系數矩陣位100階矩陣10 11101 10110 110 11 10(2)系數矩陣為40階Hilbert矩陣,引-十，向量b的第i個分量為即系數矩陣A的第i行第j列元素為n 1b =送。 j -1、分析與程序1. 平方根法函數程序如下：fun cti on x,b=p ingfan gge nfa(A,b) n=size(A);n=n (1);x=AA-1*b; disp(Matlab 自帶解即為 x);for k=1: nA(k,k)=sqrt(A (k,k);A(k+1: n,k)=A(k+1: n,k)/A (k, k); fo

3、r j=k+1: n;A(j: n,j)=A(j: n,j)-A(j: n,k)*A(j,k); endendfor j=1: n-1b(j)=b(j)/A(j,j);b(j+1:n)=b(j+1:n)-b(j)*A(j+1:n,j);endb(n)=b(n)/A(n,n);A=A;for j=n:-1:2b(j)=b(j)/A(j,j);b(1:j-1)=b(1:j-1)-b(j)*A(1:j-1,j); end b(1)=b(1)/A(1,1);disp(平方根法的解即為b);endfunction x=ave(A,b,n)求解 Ax=bL=zeros(n,n);D=diag(n,0);%

4、L 的主對角元素均S=L*D;for i=1:n為1L(i,i)=1;end for i=1:nfor j=1:nif (eig(A) In pingfanggenfa at 4In qiujie at 10 Matlab 自帶解即為 x 平方根法的解即為 bx =1.60358.96850.85621.01950.9375-50.2500-3.0000-16.000024.0000-49.5000-30.000039.000022.0000-64.0000-12.00002.000010.2500 -10.5000-1.0000-10.875083.000046.0000-98.000012

5、.0000-69.000068.000021.000017.0000-50.7188-8.7500-8.0000112.00006.0000 -68.750022.000044.0000-28.00008.0000-44.000012.00001.0e+007 *0.0000-0.00000.0001-0.0004-0.00140.0424-0.29801.1419-2.73354.2539-4.30182.7733-1.19890.5406-0.36880.3285-0.44380.4621-0.25130.05650.0000-0.00510.0071-0.0027-0.00310.003

6、6-0.00190.00090.0002-0.0002-0.00060.00040.0001-0.00020.00010.0000-0.00000.0000-0.0000-0.0000改進平方根法解得的解即為 b1.0e+024 *0.0000-0.00000.0001-0.00120.0139-0.09540.4208-1.21012.0624-1.0394-3.33436.2567-0.2463-7.45942.80303.69900.7277 -1.7484 -0.4854 -3.60100.25325.1862 -2.12991.44100.8738 -4.56541.04224.09

7、20-2.7764-2.2148-0.89530.36654.89671.04160.1281-4.3387-1.1902-2.83348.4610-3.6008、問題描述先用你所熟悉的計算機語言將算法 2.5.1編成寫通用的子程序，然后用 你編寫的程序完成下面兩個計算任務：(1) 估計5到20階Hilbert矩陣的范數條件數；(2) 設1001+aa-1 :An =-01cR-1 -1 1 1-1 -1 -1 1一先隨機的選取x Rn，并計算出b =AnX ;然后再用列主元Gauss消去法求 解該方程組，假定計算解為?。試對n從5到30估計計算解X的精度，并 且與真實相對誤差做比較。、分析與

8、程序Hilbert 矩陣:Hilbert 矩陣的分量滿足 H(i,j)=1/(i+j-1) 比如3階Hilbert矩陣是1/1 1/2 1/31/2 1/3 1/41/3 1/4 1/5程序1.用MathCAD計算方法如下:H(n) 一(n)1i j -1=normi(H(n) normiH(n)其中H (n)表示n階Hilbert矩陣，(n)即是n階Hilbert矩陣的x范數條件數2.比較估計精度與真實相對誤差用MathCAD計算方法如下：F(n)=for i 三 1 ,nA. .1i JA.11i nfor i -= 2 , ,nfor j 三 1 ,j -1A.、1 i jX(n)=fo

9、r i 三 1 ,nL(n)二Xj -1.1 9 9 9 8xfor i 三 1 ,nLifor i 三 2 ,nb(n) =F(n) X(n)U(n) =L(n)“F(n)y( n) =L( n)b (n)x(n) =U(n) 1 y(n)r(n) =b(n) F(n) x(n)_ T -1 normiF (i) maxf (i) normi(F(i)maXb (i)T(n)二for i 三 5 ,nmax(X(i) _x(i) ti 1 ”maXX(i)其中F (n)表示n階矩陣An , X (n)是任選的x Rn , b(n)即是b = Anx ,L(n)和U(n)分別是對An進行LU分

10、解得到的下三角矩陣和上三角矩陣, x(n)是用列主元Gauss消去法求得的解，P(n)是n階矩陣A計算解的精度,T(n)是n階矩陣An計算解的真實相對誤差。3估計范數條件數用MathCAD處理計算的結果如下:當矩陣階數m ：=5 .20時，,m)即是m階Hilbert矩陣的范數條件數。以矩陣階數m為橫坐標，條件范數(m)為縱坐標，其圖像如下：1 1018.(m)5 10175101520m條件數在一定程度上刻畫了擾動對方程組解的影響程度。通常，若線性方程組的系數矩陣A的條件數KA)很大，則A是病態的；反之，A是良第15頁態的。可以看出，Hilbert矩陣的條件數很大4.比較估計精度與真實相對誤

11、差 用MathCAD處理計算的結果如下:Pi表示i階矩陣A計算解的精度，在圖中以紅色的點表示，其中i =5,30 ； tj表示j階矩陣Aj計算解的真實相對誤差，在圖中以藍色的點表示，j =5 ,30。以i，j為橫坐標，Pi， tj為縱坐標，在同一圖像上表示如下:Pi3 10 72 10一710一70102030這一方法給出了計算解相對誤差的相當好的估計，可以看出，真實相 對誤差不大于計算解的精度。5.結論條件數在一定程度上刻畫了擾動對方程組解的影響程度。通常，若線 性方程組的系數矩陣A的條件數KA）很大，則A是病態的；反之，A是良 態的?？梢钥闯?，Hilbert矩陣的條件數很大，故 Hilbe

12、rt矩陣是十分病態 的。一、問題描述用你所熟悉的計算機語言編制利用 QR分解求解線性方程組和線性最 小二乘問題的通用子程序，并用你編寫的程序完成下面兩個計算任務：（1）求解第一章上機練習題中的三個線性方程組，并將所得的計算結果與 前面的計算結果相比較說明各方法的優劣；2（2）求一個二次多項式y =at bt c使在殘向量的2范數最小意義下擬合F面的數據ti-1-0.75-0.500.250.50.75yi1.000.81250.751.001.31251.752.3125（3 ）在房產估價的線性模型y =x0 awa2x2 亠亠anXn中，q,a2,a3,an分別表示稅、浴室數目占地面積、居住

13、面積、車庫數目、房屋數目、居室數目房齡、建筑類型、戶型及壁爐數目y代表房屋價格。先根據如下的28組數據求出模型中的參數的最小二乘結果。y25.929.527.925.929.929.930.928.984.982.935.931.531.030.930.028.936.941.940.543.937.537.944.537.938.936.945.841.0分析及程序利用QR分解求解程序如下:1、求a的QR分解；2、計算Cl= Qfb.3、求解上三角方程Rx =Cl得x;調試過程及實驗結果： t=-1 -0.75 -0.5 00.25 0.5 0.75; y= 1.00 0.8125 0.75

14、 1.00 1.3125 1.75 2.3125; plot(t,y,r*); lege nd(實驗數據(ti,yi); xlabel(t), ylabel(y); title(二次多項式擬合的數據點(ti,yi)的散點圖);編寫下列MATLAB程序計算f(x)在區，力)處的函數值，即輸入程序 syms a b c t=-1 -0.75 -0.5 0 0.25 0.5 0.75; fi=a.*t.A2+ b.*t+c%運行后屏幕顯示關于a，b，c的線性方程組fi =a-b+c,9/16*a-3/4*b+c,1/4*a-1/2*b+c,c,1/16*a+1/4*b+c,1/4*a+1/2*b+c

15、,9/ 16*a+3/4*b+c編寫構造殘向量2范數的MATLAB程序 y= 1.00 0.8125 0.75 1.00 1.3125 1.75 2.3125; y= 1.00 0.8125 0.75 1.00 1.3125 1.75 2.3125; fy=fi-y; fy2=fy.A2; J=sum(fy.A2);運行后屏幕顯示誤差平方和如下J=(a-b+c-1)A2+(9/16*a-3/4*b+c-13/16)A2+(1/4*a-1/2*b+c-3/4)A2+(c-1)A2+(1/16*a+1/4*b+c-21/16)A2+(1/4*a+1/2*b+c-7/4)A2+(9/16*a+3/4

16、*b+c-37/16)A2為求a,b,c使j達到最小，只需利用極值的必要條件-0-:b=0j-：c，得到關于a，b，c的線性方程組，這可以由下面的MATLAB程序完成，即輸入程序 Ja 仁diff(J,a); Ja2=diff(J,b);Ja3=diff(J,c); Jal 仁simple(Jal), Ja2仁simple(Ja2), Ja31=simple(Ja3)運行后屏幕顯示J分別對a，b，c的偏導數如下Ja11 =451/128*a-63/32*b+43/8*c-887/128Ja21 =-63/32*a+43/8*b-3/2*c-61/32Ja31 =43/8*a-3/2*b+14*

17、c-143/8解線性方程組Ja11 ，Ja21 =，Ja31 =，輸入下列程序 A=451/128, -63/32, -3/2 ;-63/32,43/8,-3/2;43/8,-3/2,14; B=887/128,61/32,143/8; C=B/A, f=poly2sym(C)運行后屏幕顯示擬合函數f及其系數C如下C =0.30810.85871.4018f =924/2999*xT+10301/11996*x+4204/2999故所求的擬合曲線為2f(x)=0.3081x0.8581x1.4018源程序： t=-1 -0.75 -0.5 0 0.25 0.5 0.75; y= 1.00 0.

18、8125 0.75 1.00 1.3125 1.75 2.3125; plot(t,y,r*); lege nd(實驗數據(ti,yi); xlabel(t), ylabel(y); title(二次多項式擬合的數據點(ti,yi)的散點圖); syms a b c t=-1 -0.75 -0.5 0 0.25 0.5 0.75; fi=a.*t.A2+ b.*t+cfi = a-b+c, 9/16*a-3/4*b+c, 1/4*a-1/2*b+c, c, 1/16*a+1/4*b+c, 1/4*a+1/2*b+c, 9/16*a+3/4*b+c y= 1.00 0.8125 0.75 1.00 1.3125 1.75 2.3125; y= 1.00 0.8125 0.75 1.00 1.3125 1.75 2.3125; fy=fi-y; fy2=fy.A2; J=sum(fy42)J =(a-b+c-1)A2+(9/16*a-3/4*b+c-13/16)A2+(1/4*a-1/2*b+c-3/4)A2+(c-1)A2+( 1/16*a+1/4*b+c-21/16)A2+(1/4*a+1/2*b+c-7/4)A2+(9/16*a+3/4*b+c-37/16)A2 Ja1=diff(J,a); Ja2=diff(J,b); Ja3=diff(J,c); Ja11=simple