2603606882數(shù)值線性代數(shù)第二版徐樹方高立張平文上機習題第三章實

1、- 1 -第三章上機習題用你所熟悉的的計算機語言編制利用qr分解求解線性方程組和線性最小二乘問題的通用子程序，并用你編制的子程序完成下面的計算任務(wù)：（1） 求解第一章上機習題中的三個線性方程組，并將所得的計算結(jié)果與前面的結(jié)果相比較，說明各方法的優(yōu)劣；（2） 求一個二次多項式，使得在殘向量的2范數(shù)下最小的意義下擬合表3.2中的數(shù)據(jù)；表 3.2ti-1-0.75 -0.5 0 0.250.50.75yi10.81250.7511.31251.752.3125（3） 在房產(chǎn)估價的線性模型中，分別表示稅、浴室數(shù)目、占地面積、車庫數(shù)目、房屋數(shù)目、居室數(shù)目、房齡、建筑類型、戶型及壁爐數(shù)目，代表房屋價格?，F(xiàn)

2、根據(jù)表3.3和表3.4給出的28組數(shù)據(jù)，求出模型中參數(shù)的最小二乘結(jié)果。（表3.3和表3.4見課本p99-100）解 分析：（1） 計算一個householder變換h：由于，則計算一個householder變換h等價于計算相應的。其中。在實際計算中，為避免出現(xiàn)兩個相近的數(shù)出現(xiàn)的情形，當時，令；為便于儲存，將規(guī)格化為，相應的，變?yōu)闉榉乐挂绯霈F(xiàn)象，用代替（2） qr分解：利用householder變換逐步將轉(zhuǎn)化為上三角矩陣，則有- 8 -，其中，。在實際計算中，從，若，依次計算對應的即對應的，將儲存到，儲存到，迭代結(jié)束后再次計算，有，（時）（3） 求解線性方程組或最小二乘問題的步驟為i 計算的qr

3、分解；ii 計算，其中iii 利用回代法求解上三角方程組（4）對第一章第一個線性方程組，由于r的結(jié)果最后一行為零，故使用前代法時不計最后一行，而用運行結(jié)果計算。運算matlab程序為 1 計算householder變換 v,belta=house(x)function v,belta=house(x) n=length(x); x=x/norm(x,inf); sigma=x(2:n)*x(2:n); v=zeros(n,1); v(2:n,1)=x(2:n); if sigma=0 belta=0; else alpha=sqrt(x(1)2+sigma); if x(1)=0 v(1)=x

4、(1)-alpha; else v(1)=-sigma/(x(1)+alpha); end belta=2*v(1)2/(sigma+v(1)2); v=v/v(1,1); endend2 計算的qr分解 q,r=qrfenjie(a)function q,r=qrfenjie(a)m,n=size(a);q=eye(m);for j=1:n if jm v,belta=house(a(j:m,j); h=eye(m-j+1)-belta*v*v; a(j:m,j:n)=h*a(j:m,j:n); d(j)=belta; a(j+1:m,j)=v(2:m-j+1); endendr=triu(

5、a(1:n,:);for j=1:n if jm h=eye(m); temp=1;a(j+1:m,j); h(j:m,j:m)=h(j:m,j:m)-d(j)*temp*temp; q=q*h; endendend3 解下三角形方程組的前代法 x=qiandaifa(l,b)function x=qiandaifa(l,b)n=length(b);for j=1:n-1 b(j)=b(j)/l(j,j); b(j+1:n)=b(j+1:n)-b(j)*l(j+1:n,j);endb(n)=b(n)/l(n,n);x=b; end4 求解第一章上機習題中的三個線性方程組 ex3_1clear;

6、clc;%第一題a=6*eye(84)+diag(8*ones(1,83),-1)+diag(ones(1,83),1);b=7;15*ones(82,1);14;n=length(a);%qr分解q,r=qrfenjie(a);c=q*b;x1=huidaifa(r(1:n-1,1:n-1),c(1:n-1);x1(n)=c(n)-r(n,1:n-1)*x1;%不選主元gauss消去法l,u=gaussla(a);x1_1=gauss(a,b,l,u);%列主元gauss消去法l,u,p=gausscol(a);x1_2=gauss(a,b,l,u,p);%解的比較figure(1);sub

7、plot(1,3,1);plot(1:n,x1);title(qr分解);subplot(1,3,2);plot(1:84,x1_1);title(gauss);subplot(1,3,3);plot(1:84,x1_2);title(pgauss);%第二題第一問a=10*eye(100)+diag(ones(1,99),-1)+diag(ones(1,99),1);b=round(100*rand(100,1);n=length(a);%qr分解tic;q,r=qrfenjie(a);c=q*b;x2=huidaifa(r,c);toc;%不選主元gauss消去法tic;l,u=gauss

8、la(a);x2_1=gauss(a,b,l,u);toc;%列主元gauss消去法tic;l,u,p=gausscol(a);x2_2=gauss(a,b,l,u,p);toc;%平方根法tic;l=cholesky(a);x2_3=gauss(a,b,l,l);toc;%改進的平方根法tic;l,d=ldlt(a);x2_4=gauss(a,b,l,d*l);toc;%解的比較figure(2);subplot(1,5,1);plot(1:n,x2);title(qr分解);subplot(1,5,2);plot(1:n,x2_1);title(gauss);subplot(1,5,3);

9、plot(1:n,x2_2);title(pgauss);subplot(1,5,4);plot(1:n,x2_3);title(平方根法);subplot(1,5,5);plot(1:n,x2_4);title(改進的平方根法);%第二題第二問a=hilb(40);b=sum(a);b=b;n=length(a);q,r=qrfenjie(a);c=q*b;x3=huidaifa(r,c);%不選主元gauss消去法l,u=gaussla(a);x3_1=gauss(a,b,l,u);%列主元gauss消去法l,u,p=gausscol(a);x3_2=gauss(a,b,l,u,p);%平

10、方根法l=cholesky(a);x3_3=gauss(a,b,l,l);%改進的平方根法l,d=ldlt(a);x3_4=gauss(a,b,l,d*l);%解的比較figure(3);subplot(1,5,1);plot(1:n,x3);title(qr分解);subplot(1,5,2);plot(1:n,x3_1);title(gauss);subplot(1,5,3);plot(1:n,x3_2);title(pgauss);subplot(1,5,4);plot(1:n,x3_3);title(平方根法);subplot(1,5,5);plot(1:n,x3_4);title(改

11、進的平方根法);5 求解二次多項式 ex3_2clear;clc;t=-1 -0.75 -0.5 0 0.25 0.5 0.75;y=1 0.8125 0.75 1 1.3125 1.75 2.3125;a=ones(7,3);a(:,1)=t.2;a(:,2)=t;q,r=qrfenjie(a);q1=q(:,1:3);c=q1*y;x=huidaifa(r,c)6 求解房產(chǎn)估價的線性模型 ex3_3clear;clc;a=xlsread(e:temporary專業(yè)課數(shù)值代數(shù)cha3_3_4.xls,a2:l29);y=xlsread(e:temporary專業(yè)課數(shù)值代數(shù)cha3_3_4.x

12、ls,m2:m29);q,r=qrfenjie(a);q1=q(:,1:12);c=q1*y;x=huidaifa(r,c);x=x計算結(jié)果為（1）第一章上機習題中的三個線性方程組結(jié)果對比圖依次為以第二個線性方程組為例，比較各方法的運行速度。依次為qr分解，不選主元的gauss消去法，列主元gauss消去法，平方根法，改進的平方根法。elapsed time is 0.034588 seconds.elapsed time is 0.006237 seconds.elapsed time is 0.009689 seconds.elapsed time is 0.030862 seconds.elapsed time is 0.007622 seconds.（2） 二次多項式的系數(shù)為x = 1.0000 1.0000 1.0000（3） 房產(chǎn)估價的線性模型的系數(shù)為x = columns 1 through 6 2.0775 0.7189 9.6802 0.1535 1