百度文库-
百度文库-让每个人平等地提升自我
PAGE
PAGE #
上机习题
1?先用你所熟悉的的计算机语言将不选主元和列主元 Gauss消去法编写成通用的子程序；
后用你编写的程序求解 84阶方程组；最后将你的计算结果与方程的精确解进行比较，并就 此谈谈你对Gauss消去法的看法。
Sol:
先用matlab将不选主元和列主元 Gauss消去法编写成通用的子程序，得到 L,U,P
不选主元Gauss消去法：L,U GaussLA(A)得到L,U满足A LU
列主元 Gauss消去法：L,U ,P GaussCol(A)得到 L,U , P 满足 PA LU
用前代法解Ly b or Pb，得y
用回代法解Ux y，得x
求解程序为x Gauss代b,L,U,P ( P可缺省，缺省时默认为单位矩阵)
计算脚本为ex1_1
代码
%!法(计算三角分解：Gauss消去法)
fun ctio n [L,U]=GaussLA(A)
n=len gth(A);
for k=1:n-1
A(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);
end
百度文库-
百度文库-让每个人平等地提升自我
PAGE
PAGE #
U=triu(A);
L=tril(A);
L=L-diag(diag(L))+diag( on es(1, n));
end
%!法计算列主元三角分解：列主元 Gauss消去法)
fun ctio n [L,U,P]=GaussCol(A)
n=len gth(A);
for k=1:n-1
[s,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)~=0
A(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 ;
end
end
L=tril(A);U=triu(A);L=L-diag(diag(L))+diag(o nes(1, n));
百度文库-
百度文库-让每个人平等地提升自我
PAGE
PAGE #
P=eye( n);
for i=1:n-1
temp=P(i,:);
P(i,:)=P(u(i),:);
P(u(i),:)=temp;
end
end
%高斯消去法解线性方程组
fun ctio n x=Gauss(A,b,L,U,P)
if nargin<5
P=eye(le ngth(A));
end
n=len gth(A);
b=P*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);
end
b(n )=b( n)/L( n,n);
y=b;
for j=n:-1:2
y(j)=y(j)/U(j,j);
百度文库-
百度文库-让每个人平等地提升自我
PAGE
PAGE #
y(1:j-1)=y(1:j-1)-y(j)*U(1:j-1,j);
end
y(1)=y(1)/U(1,1);
x=y;
end ex1_1
clc;clear;
%第一题
A=6*eye(84)+diag(8*o nes(1,83),-1)+diag(o nes(1,83),1);
b=[7;15*o nes(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);
'o-' );title( 'Gauss');'.-' );title( 'PGauss');'*-' );title('
'o-' );title( 'Gauss');
'.-' );title( 'PGauss');
'*-' );title('精确解');
subplot(1,3,1);plot(1:84,x1_1,
subplot(1,3,2);plot(1:84,x1_2,
subplot(1,3,3);plot(1:84,o nes(1,84),
结果为（其中Gauss表示不选主元的Gauss消去法，PGauss表示列主元Gauss 消去法，精确解为1, ,1184）： \ /
百度文库-
百度文库-让每个人平等地提升自我
PAGE
PAGE #
由图，显然列主元消去法与精确解更为接近, 大，且不如列主元消去法稳定。不选主元的Gauss消去法误差比列主元消去法Gauss消去法重点在于
由图，显然列主元消去法与精确解更为接近, 大，且不如列主元消去法稳定。
不选主元的Gauss消去法误差比列主元消去法
Gauss消去法重点在于 A的分解过程，无论
A如何分解，后面两步的运算过程不变。
2?先用你所熟悉的的计算机语言将平方根法和改进的平方根法编写成通用的子程序； 然后用
你编写的程序求解对称正定方程组 Ax=b。
Sol:
先用matlab将平方根法和改进的平方根法编写成通用的子程序，得到 L, (D)：
平方根法：L=Cholesky(A)
改进的平方根法：[L,D]=LDLt(A)
求解得Ly b
求解得 LTx y or DLTx y
求解程序为x=Gauss(A,b,L,U,P)( U LT or U DLT，P此时缺省，缺省时默认为单
位矩阵)
百度文库-
百度文库-让每个人平等地提升自我
PAGE
PAGE #
(3)计算脚本为ex1_2
代码
%|法(计算 Cholesky分解：平方根法)
fun ctio n L=Cholesky(A)
n=len gth(A);
for k=1:n
A(k,k)=sqrt(A(k,k));
A(k+1: n,k)=A(k+1: n,k)/A(k,k);
for j=k+1:n
A(j: n,j)=A(j: n,j)-A(j: n,k)*A(j,k);
end
end
L=tril(A);
end
%十算LDL '分解：改进的平方根法
fun ctio n [L,D]=LDLt(A)
n=len gth(A);
for j=1:n
for i=1:n
v(i,1)=A(j,i)*A(i,i);
end
A(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 PAGE
end
PAGE #
end
L=tril(A);
D=diag(diag(A));
L=L-diag(diag(L))+diag( on es(1, n)); end
%高斯消去法解线性方程组
fun ctio n x=Gauss(A,b,L,U,P)
if nargin<5
P=eye(le ngth(A));
end
n=len gth(A);
b=P*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);
end
b(n )=b( n)/L( n,n);
y=b;
for j=n:-1:2
y(j)=y(j)/U(j,j);
y(1:j-1)=y(1:j-1)-y(j)*U(1:j-1,j);
百度文库-
百度文库-让每个人平等地提升自我
PAGE
PAGE #
y(1)=y(1)/U(1,1);
x=y;
end /
ex1_2
嘟二题
嘟一问
A=10*eye(100)+diag(o nes(1,99),-1)+diag(o nes(1,99),1);
b=ro un d(100*ra nd(100,1));
那方根法
L=Cholesky(A);
x1_2_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=Gauss(A,b,L,D*L');
百度文库-
百度文库-让每个人平等地提升自我
PAGE
PAGE #
结果分别为
x1_2_1_1 =
百度文库-
百度文库-让每个人平等地提升自我
\z PAGE
\z
PAGE #
百度文库-
百度文库-让每个人平等地提升自我
\z PAGE
\z
PAGE #
百度文库-
百度文库-让每个人平等地提升自我
\z PAGE
\z
PAGE #
百度文库-
百度文库-让每个人平等地提升自我
PAGE
PAGE #
x1_2_1_2 =
百度文库-
百度文库-让每个人平等地提升自我
\z PAGE
\z
PAGE #
百度文库-
百度文库-让每个人平等地提升自我
\z PAGE
\z
PAGE #
百度文库-
百度文库-让每个人平等地提升自我
PAGE
PAGE #
x1_2_2_1 =
百度文库-
百度文库-让每个人平等地提升自我
PAGE
PAGE #
+07 *
百度文库-
百度文库-让每个人平等地提升自我
1818
1818
x1_2_2_2 =
百度文库-
百度文库-让每个人平等地提升自我
\z1919
\z
1919
百度文库-
百度文库-让每个人平等地提升自我
2020
2020
然后评价各个方法的优Gauss消去法3?用第1题的程序求解第2
然后评价各个方法的优
Gauss消去法
Sol:
Gauss表示不选主元的Gauss消去法，PGauss表示列主元
计算脚本为：
嘟三题
嘟一问
A=10*eye(100)+diag(o nes(1,99),-1)+diag(o nes(1,99),1);
b=ro un d(100*ra nd(100,1));
%F选主元Gauss消去法
[L,U]=GaussLA(A);
x1_3_1_ 仁Gauss(A,b,L,U);
%^主元Gauss消去法
[L,U,P]=GaussCol(A);
x1_3_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);plot(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');
2121
平方根法1
10
20
30
40
50
60
90
10
5
0
-5
0
10
5
0
-5
0 10
改进的平方根法1
20 30 40 50 60 70 80
90
Gaussl
10
5
0
-5
0
5
0
-5
0
200
100
PGauss1
100
10
20
30 40 50 60 70 80 90 100
10
/平方根法2
7
10
15
20
Gauss2
25
30 35 40
10
-5
5
0
0 10 20 \ 30 40 50 60 70 80
90
改进的平方根法2
100
0 5
10
15
20
PGauss2
25
30 35
40
-200
0
1
1
1
1
1 I J \ J
1
0
L
I
I
I ■<
■ r \
1J
| X
jf ~
V
1 i
r
r
[
r
-500
r
r
[
[
5
10
15 20 25 30
35
40
0
5
10 15
20 25 30 35
0
500
40
subplot(4,2,8);plot(y2,x1_3_2_2);title( 'PGauss2');
平方根法和改进的平法根法计算量更小，计算过程稳定，但使用范围窄；
不选主元和列主元的 Gauss消去法计算量较大，但适用范围广。
例题考虑对称正定线性方程组 Ax=b,其中向量b是随机生成的，其元素是服从区间 ［0,1］上
均匀分布的随机数，矩阵 A LLT,这里L是随机生成的一个下三角矩阵，其元素是服从区
间［1,2］上均匀分布的随机数。
对n=10, 20，...,500分别应用 Gauss消去法、列主元Gauss消去法和Cholesky分解法求解该
方程组，画出它们所用的 CPU时间，其中"Gauss"表示Gauss消去法、“ PGausg表示列主
元Gauss消去法，“ Cholesky"表示Cholesky分解法。
Sol:
经试验知，对应课本上图所示的 Cholesky分解法应为改进后的 Cholesky分解法即LDLT分
解。
2222
百度文库-
百度文库-让每个人平等地提升自我
end2323
end
2323
此处所用的CPU时间利用cputime测量。
计算脚本为eg1_3_2
clc;clear; /
for i=1:50;
n=i*10;
b=ra nd(n ,1);
L=tril(u nifrnd(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 ,L 3');
t3(i)=cputime-t3(i);
百度文库-
百度文库-让每个人平等地提升自我
2424
2424
N=10:10:500;
plot(N,t1, 'o-' ,N,t2, '.-' ,N,t3, '*-');
legend( 'Gauss' , 'PGauss' , 'Cholesky');
结果为
1.5
0.5
50
350
150 200
250 300
————Gauss ,,
—*— PGauss
Cholesky
-
fl _
?
*
H0
?
* *
#
*
a
肩机战就!卄卩;V；.t
100
500
400 450