1、上机习题1.先用你所熟悉的的计算机语言将不选主元和列主元 Gauss 消去法编写成通用的子程序;然后用你编写的程序求解 84 阶方程组;最后将你的计算结果与方程的精确解进行比较,并就此谈谈你对 Gauss 消去法的看法。Sol:(1 )先用 matlab 将不选主元和列主元 Gauss 消去法编写成通用的子程序,得到 :PUL,不选主元 Gauss 消去法: 得到 满足)(,AGausLUU,A列主元 Gauss 消去法: 得到 满足,ColPLP,L(2 ) 用前代法解 ,得bory y用回代法解 ,得Ux求解程序为 ( 可缺省,缺省时默认为单位矩阵)PLbAGaus,(3 )计算脚本为 e
2、x1_1代码%算法 1.1.3(计算三角分解:Gauss 消去法)function L,U=GaussLA(A)n=length(A);for k=1:n-1A(k+1:n,k)=A(k+1:n,k)/A(k,k);A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-A(k+1:n,k)*A(k,k+1:n);endU=triu(A);L=tril(A);L=L-diag(diag(L)+diag(ones(1,n);end%算法 1.2.2(计算列主元三角分解:列主元 Gauss 消去法)function L,U,P=GaussCol(A)n=length(A);for k=1:n-
3、1s,t=max(abs(A(k:n,k);p=t+k-1;temp=A(k,1:n);A(k,1:n)=A(p,1:n);A(p,1:n)=temp;u(k)=p;if A(k,k)=0A(k+1:n,k)=A(k+1:n,k)/A(k,k);A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-A(k+1:n,k)*A(k,k+1:n);else break;endendL=tril(A);U=triu(A);L=L-diag(diag(L)+diag(ones(1,n);P=eye(n);for i=1:n-1temp=P(i,:);P(i,:)=P(u(i),:);P(u(i),
4、:)=temp;endend%高斯消去法解线性方程组function x=Gauss(A,b,L,U,P)if nargin5P=eye(length(A);endn=length(A);b=P*b;for j=1:n-1b(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);y=b;for j=n:-1:2y(j)=y(j)/U(j,j);y(1:j-1)=y(1:j-1)-y(j)*U(1:j-1,j);endy(1)=y(1)/U(1,1);x=y;endex1_1clc;clear;%第一题A=6*e
5、ye(84)+diag(8*ones(1,83),-1)+diag(ones(1,83),1);b=7;15*ones(82,1);14;%不选主元 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);%解的比较subplot(1,3,1);plot(1:84,x1_1,o-);title(Gauss);subplot(1,3,2);plot(1:84,x1_2,.-);title(PGauss);subplot(1,3,3);plot(1:84,on
6、es(1,84),*-);title(精确解);结果为(其中 Gauss 表示不选主元的 Gauss 消去法, PGauss 表示列主元 Gauss消去法,精确解为 ):841,0 50 100-6-4-20246x 108 Gauss0 50 10011111111 PGauss0 50 10000.20.40.60.811.21.41.61.82 平 平 平由图,显然列主元消去法与精确解更为接近,不选主元的 Gauss 消去法误差比列主元消去法大,且不如列主元消去法稳定。Gauss 消去法重点在于 A 的分解过程,无论 A 如何分解,后面两步的运算过程不变。2.先用你所熟悉的的计算机语言将
7、平方根法和改进的平方根法编写成通用的子程序;然后用你编写的程序求解对称正定方程组 Ax=b。Sol:(1 )先用 matlab 将平方根法和改进的平方根法编写成通用的子程序,得到 L,(D):平方根法:L=Cholesky(A)改进的平方根法:L,D=LDLt(A)(2 ) 求解得 bLy求解得 yxDorxTT 求解程序为 x=Gauss(A,b,L,U,P)( , 此时缺省,缺省时默认为TTDLUor P单位矩阵)(3 ) 计算脚本为 ex1_2代码%算法 1.3.1(计算 Cholesky 分解:平方根法)function L=Cholesky(A)n=length(A);for k=1
8、:nA(k,k)=sqrt(A(k,k);A(k+1:n,k)=A(k+1:n,k)/A(k,k);for j=k+1:nA(j:n,j)=A(j:n,j)-A(j:n,k)*A(j,k);endendL=tril(A);end%计算 LDL分解:改进的平方根法function L,D=LDLt(A)n=length(A);for j=1:nfor i=1:nv(i,1)=A(j,i)*A(i,i);endA(j,j)=A(j,j)-A(j,1:j-1)*v(1:j-1,1);A(j+1:n,j)=(A(j+1:n,j)-A(j+1:n,1:j-1)*v(1:j-1,1)/A(j,j);end
9、L=tril(A);D=diag(diag(A);L=L-diag(diag(L)+diag(ones(1,n);end%高斯消去法解线性方程组function x=Gauss(A,b,L,U,P)if nargin5P=eye(length(A);endn=length(A);b=P*b;for j=1:n-1b(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);y=b;for j=n:-1:2y(j)=y(j)/U(j,j);y(1:j-1)=y(1:j-1)-y(j)*U(1:j-1,j);endy(
10、1)=y(1)/U(1,1);x=y;endex1_2%第二题%第一问A=10*eye(100)+diag(ones(1,99),-1)+diag(ones(1,99),1);b=round(100*rand(100,1);%平方根法L=Cholesky(A);x1_2_1_1=Gauss(A,b,L,L);%改进的平方根法L,D=LDLt(A);x1_2_1_2=Gauss(A,b,L,D*L);%第二问A=hilb(40);b=sum(A);b=b;%平方根法L=Cholesky(A);x1_2_2_1=Gauss(A,b,L,L);%改进的平方根法L,D=LDLt(A);x1_2_2_2
11、=Gauss(A,b,L,D*L);结果分别为x1_2_1_1 =7.25868.4143-0.40138.59845.41770.22492.33364.43898.27728.7890-0.16678.87848.38243.29787.64010.30143.34578.24186.23688.39065.8575-0.96567.79817.98425.36016.41436.49662.62016.30240.35637.1342-0.69852.85060.19270.22277.58015.97621.65839.4409-1.06774.23622.70606.70377.25
12、700.72654.47783.49595.56375.86756.76141.51806.05825.90020.93970.70314.02949.00321.93825.61500.91207.26521.43604.37495.81467.47918.39424.57890.81691.25231.66038.14480.89157.94010.70758.98492.44371.57771.77905.63193.90182.35067.59254.72454.16278.64831.35436.80876.55892.60275.41400.25770.00904.65226.46
13、858.6626-0.09485.28564.2385-0.67063.4671x1_2_1_2 =7.25868.4143-0.40138.59845.41770.22492.33364.43898.27728.7890-0.16678.87848.38243.29787.64010.30143.34578.24186.23688.39065.8575-0.96567.79817.98425.36016.41436.49662.62016.30240.35637.1342-0.69852.85060.19270.22277.58015.97621.65839.4409-1.06774.236
14、22.70606.70377.25700.72654.47783.49595.56375.86756.76141.51806.05825.90020.93970.70314.02949.00321.93825.61500.91207.26521.43604.37495.81467.47918.39424.57890.81691.25231.66038.14480.89157.94010.70758.98492.44371.57771.77905.63193.90182.35067.59254.72454.16278.64831.35436.80876.55892.60275.41400.257
15、70.00904.65226.46858.6626-0.09485.28564.2385-0.67063.4671x1_2_2_1 =1.0e+07 *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.0036-0.00190.00090.0002-0.0002-0.00060.00040.0001-0.00020.0
16、0010.0000-0.00000.0000-0.0000-0.0000x1_2_2_2 =1.00001.00000.99981.00111.00640.86811.8034-1.56935.5763-2.5315-1.769310.4883-6.28070.5882-4.715722.8299-19.91348.703210.3265-25.214010.028212.3882-1.942514.1891-12.0552-0.5803-12.47918.56529.8724-10.550216.3871-5.813213.421611.1767-64.315446.383712.6957-
17、21.755612.1204-1.93423.用第 1 题的程序求解第 2 题的两个方程组并比较所有的计算结果,然后评价各个方法的优劣。Sol:Gauss 表示不选主元的 Gauss 消去法,PGauss 表示列主元 Gauss 消去法。计算脚本为:%第三题%第一问A=10*eye(100)+diag(ones(1,99),-1)+diag(ones(1,99),1);b=round(100*rand(100,1);%不选主元 Gauss 消去法L,U=GaussLA(A);x1_3_1_1=Gauss(A,b,L,U);%列主元 Gauss 消去法L,U,P=GaussCol(A);x1_3
18、_1_2=Gauss(A,b,L,U,P);%第二问A=hilb(40);b=sum(A);b=b;%不选主元 Gauss 消去法L,U=GaussLA(A);x1_3_2_1=Gauss(A,b,L,U);%列主元 Gauss 消去法L,U,P=GaussCol(A);x1_3_2_2=Gauss(A,b,L,U,P);ex1_2;y1=1:100;y2=1:40;subplot(4,2,1);plot(y1,x1_2_1_1);title(平方根法 1);subplot(4,2,2);plot(y1,x1_2_1_2);title(改进的平方根法 1);subplot(4,2,3);plo
19、t(y1,x1_3_1_1);title(Gauss1);subplot(4,2,4);plot(y1,x1_3_1_2);title(PGauss1);subplot(4,2,5);plot(y2,x1_2_2_1);title(平方根法 2);subplot(4,2,6);plot(y2,x1_2_2_2);title(改进的平方根法 2);subplot(4,2,7);plot(y2,x1_3_2_1);title(Gauss2);subplot(4,2,8);plot(y2,x1_3_2_2);title(PGauss2);平方根法和改进的平法根法计算量更小,计算过程稳定,但使用范围窄
20、;不选主元和列主元的 Gauss 消去法计算量较大,但适用范围广。例题 1.3.2考虑对称正定线性方程组 Ax=b,其中向量 b 是随机生成的,其元素是服从区间0,1上均匀分布的随机数,矩阵 ,这里 L 是随机生成的一个下三角矩阵,其元素是服从区间TA1,2上均匀分布的随机数。对 n=10,20,.,500 分别应用 Gauss 消去法、列主元 Gauss 消去法和 Cholesky 分解法求解该方程组,画出它们所用的 CPU 时间,其中“Gauss”表示 Gauss 消去法、 “PGauss”表示列主元 Gauss 消去法, “Cholesky”表示 Cholesky 分解法。Sol:0 1
21、0 20 30 40 50 60 70 80 90 100-50510 平平平平10 10 20 30 40 50 60 70 80 90 100-50510 平平平平平平平10 10 20 30 40 50 60 70 80 90 100-50510 Gauss10 10 20 30 40 50 60 70 80 90 100-50510 PGauss10 5 10 15 20 25 30 35 40-505x 107 平平平平20 5 10 15 20 25 30 35 40-100-50050 平平平平平平平20 5 10 15 20 25 30 35 40-2000200 Gauss2
22、0 5 10 15 20 25 30 35 40-5000500 PGauss2经试验知,对应课本上图 1.1 所示的 Cholesky 分解法应为改进后的 Cholesky 分解法即分解。TLD此处所用的 CPU 时间利用 cputime 测量。计算脚本为 eg1_3_2clc;clear; for i=1:50;n=i*10;b=rand(n,1);L=tril(unifrnd(1,2,n,n);A=L*L;t1(i)=cputime;L1,U1=GaussLA(A);x1=Gauss(A,b,L1,U1);t1(i)=cputime-t1(i);t2(i)=cputime;L2,U2,P2=GaussCol(A);x2=Gauss(A,b,L2,U2,P2);t2(i)=cputime-t2(i);t3(i)=cputime;L3=LDLt(A);x3=Gauss(A,b,L3,L3);t3(i)=cputime-t3(i);endN=10:10:500;plot(N,t1,o-,N,t2,.-,N,t3,*-);legend(Gauss,PGauss,Cholesky);结果为0 50 100 150 200 250 300 350 400 450 50000.511.5GaussPGaussCholesky