科学计算方法6(高斯消元法)
合集下载
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
注释: 8是主元(pivot), -0.5是乘子(multiplier)。 20/29
更加简洁的表示:
8 6 2 28 8 6 2 28 8 6 2 28
4
11
7 40 0
8
6 26 0
8
6 26
4 7 6 33 0 4 5 19 0 0 2 6
8 6 2 1 0 0 8 6 2 A= 4 11 7 = 0.5 1 0 0 8 6 =LU
28/29
My own faith has been strong for years. Back in 1985, I made a $100 bet with Peter Alfeld that a fast matrix inverse would be found by 1995. It wasn’t, so I paid him $100 and renewed the bet for another decade. It still hadn’t been found by 2005, so I paid him a further $100 and we renewed once more. One day I will win!—or my heirs?
我们经常使用记号O表示算法是“是多少阶的”。
因此高斯消元法为O(n3/3)的算法。复杂度(complexity)
是衡量算法性能优劣的主要指标。
26/29
例1. 在一台特定的计算机上, 估计用高斯消元法 求解5000*5000的线性方程组所需时间。
高斯消元法求解50*50的线性方程组所需时间为t1。
22/29
代码:
%elimination phase
for k = 1:n-1 %pivot row
for i= k+1:n %rows to be transformed
if A(i,k) ˜= 0 %multiplier
lambda = A(i,k)/A(k,k);
A(i,k+1:n) = A(i,k+1:n) - lambda*A(k,k+1:n);
回代阶段 0.1242 0.0147 0.0015 0.0001
高斯消元法由不相等的两部分组成:工作量相对多
的消元阶段和工作量相对少的回代阶段。对于大的n, 幂次较低的n相比而言可以忽略不计,最高次幂决定了 当n趋近于无穷时的极限形态。换而言之, 对大的n, 低 阶项对算法的执行时间的估计没有太大的影响, 当仅需 近似估计执行时间时可以忽略不计。
11/29
(Cramer法则) 设n阶矩阵A可逆,则线性方程组Ax b
有唯一解x ( x1, x2 , , xn )T ,其中
xj
det Aj det A
( j 1,2,
, n)
计算一个n阶行列式需要的乘法次数为(n -1)n。
以n 30为例,总的乘法次数大约为7.7e+33。
7.7e+33/365(天) / 24(小时) / 3600(秒) / 109 2.4e+18(年) 7.7e+33/365(天) / 24(小时) / 3600(秒)（/ 4.7 1015） 5.2e+10(年)
36
51
13
52
34
74
0
7
1.1
减肥所要求 每日营养量
33 45 3
1/29
2/29
3/29
4/29
5/29
6/29
大规模代数系统的求解是科学与工程计算 的基础和关键,其计算量也占有非常大的比重, 有的甚至占整个计算的80%以上！
Gene Golub (1932—2007)
美国科学院,工程院、艺术与科学学院院士 以及瑞典皇家工程科学院院士
503 / 50003
3 /3
t503 / t 50003
3 /3
t1 t2
t2 1003 t1
27/29
回顾:
高斯消元法是时间复杂度为O(n3)的算法, 最 终在有限步之内给出一个解。因此高斯消元法 叫做求解线性方程组的直接法。
理论上直接法在有限步之内给出精确解(当 然在一有限精度的计算机上执行时得到的结果 将只是近似的)。
x3 y3 / 2 3 x2 ( y2 6x3 ) / 8 1 x1 ( y1 6 x2 2 x3 ) / 8 2
16/29
Erdős 相信上帝手中有一本包含世间所有精妙 证明的天书。上帝相信这本书在 Gauss 手上。
17/29
高斯消元法
高斯消元法是求解线性方程组的经典方法之一。高 斯消元法由两个部分构成: 消元阶段和回代阶段 。消元 阶段的目的是将一般线性方程组转化成容易求解的上三 角线性方程组。
a11 a12 a1n x1 b1
a21
a22
a2n
x2
b2
an1
an2
ann
xn
bn
Ax=b
9/29
存在性
秩是定量刻画矩阵性质的指标, 相对于矩阵行列式更加精细。
增广矩阵:
a11 a12
[
A,
b]
a21
a22
an1 an2
a1n b1
a2n
b2
ann
y1
28
0.5 y1 + y2 40
0.5y1 0.5 y2 y3 33
y1 28 y2 (40 0.5 y1 )2 26 y3 33 0.5 y1 0.5 y2 6
然后求解上三角方程组Ux = y :
8 x1 6 x2 2 x3 y1 8x2 6x3 y2 2x3 y3
14/29
例 求解方程组Ax = b, 其中
8 6 2 28 A 4 11 7 , b 40
4 7 6 33
假设系数矩阵A存在如下LU分解:
1 0 0 8 6 2 A LU 0.5 1 0 0 8 6
0.5 0.5 1 0 0 2
15/29
解: 首先求解下三角线性方阵组Ly = b:
设3种食物每100克中蛋白质、碳水化合物和脂肪的含量 如下表; 表中还给出了美国流行的剑桥大学医学院的简捷营养处方。 现在的问题是:如果用这3种食物作为每天的主要食物,那么它们的 用量应各取多少才能全面准确地实现这个营养。
营养
蛋白质 碳水化合物 脂肪
每100g食物所含营养(g)
脱脂牛奶 大豆面粉
乳清
b(i)= b(i) - lambda*b(k);
end
end end
xk [bk (ak,k1 xk1 ak ,n xn )] / ak ,k (k n 1, 1)
%back-substitution phase
for k = n:-1:1
b(k) = (b(k) - A(k,k+1:n)*x(k+1:n))/A(k,k);
7/29
求解
a11 x1 a12 x2 a1n xn b1 a21x1a22 x2 a2n xn b2 an1 x1 an2 x2 ann xn bn
8/29
➢线性方程组的矩阵形式
a11 x1 a12 x2 a1n xn b1 a21x1a22 x2 a2n xn b2 an1 x1 an2 x2 ann xn bn
12/29
求解如下下三角线性方程组Lx=c:
l11 x1 c1 l21 x1 l22 x2 c2 l31 x1 l32 x2 l33 x3 c3
从第一个方程开始, 从前往后每次求解一个方程, 则下三 角方程组的解可以表示如下
x1 c1 / l11 x2 (c2 l21 x1 ) / l22 x3 (c3 l31 x1 l32 x2 ) / l33
例如
8x1 6x2 2x3 28 (a) 4x1 11x2 7 x3 40 (b)
4x1 7 x2 6x3 33 (c)
18/29
如何将线性方程组化为更容易求解的线性方程组且 不改变线性方程组的解？
➢ 交换两个方程; ➢ 某一方程乘以非零常数c; ➢ 把某一方程的c倍加于另一方程。
O((n-1)n!) (Cramer) O(n3) (Gauss) O(n2.81) (Strassen) O(n2.373) O(n2+eps)?
Ref: The smart money’s on numerical analysis, SIAM News, 2012
29/29
19/29
消元阶段:
经过第一轮消元:
8x1 6x2 2x3 28
(a)
8 x2 6 x3 26 (b)
4 x2 5 x3 19
(c)
注释: 8是主元(pivot), -0.5和0.5是乘子(multiplier)。
经过第二轮消元:
8x1 6x2 2x3 28 (a) 8x2 6x3 26 (b) 2x3 6 (c)
共计： 2n3/3+n2/2-7n/6
回代阶段:
除法次数 n次, 乘法次数 n(n-1)/2次,
加(减)法次数 n(n-1)/2次
共计： n2
25/29
n=10 n=100 n=1000 n=10000
消元阶段 0.8282 0.0621 0.0145 0.9782 0.0073 0.0001 0.9978 0.0007 0.0000 0.9998 0.0000 0.0000
4 7 6 0.5 0.5 1 0 0 2
Guass消元法和LU分解的关系？ 21/29
回代阶段:
从最后一个方程开始, 从后往前每次求解一个方程, 则 三角方程组的解可以表示如下
x3 6 / 2 3 x2 (26 6x3 ) / 8 1 x1 (28 6 x2 2 x3 ) / 8 2
13/29
类似的, 求解如下上三角线性方程组Ux=c:
u11 x1 u12 x2 u13 x3 c1 u22 x2 u23 x3 c2 u33 x3 c3
从最后一个方程开始, 从后往前每次求解一个方程, 则上 三角方程组的解可以表示如下
x3 c3 / u33 x2 (c2 u23 x3 ) / u22 x1 (c1 u13 x3 u12 x2 ) / u11
end
23/29
24/29
12 22 32
消元阶段:
n 2 n(n1)(2n1) 6
除法次数(n-1)+ (n-2)+···+1= n(n-1)/2,
乘法次数(n-1)*n+ (n-2)* (n-1) +···+1*2= (n3-n)/3,
加(减)法次数(n-1)*n+ (n-2)* (n-1) +···+1*2= (n3-n)/3。
bn
如果rank([A,b])=rank(A), 则线性方程组的有解。
10/29
唯一性 如果系数矩阵非奇异即|A| ≠ 0, 则线性方程组的解唯一。 如果系数矩阵奇异即|A| = 0,则线性方程组的有无穷多 解或无解。例如, 2x + y = 3和4x + 2y = 6有无穷多解; 而 2x + y = 3和4x + 2y = 0则无解。
更加简洁的表示:
8 6 2 28 8 6 2 28 8 6 2 28
4
11
7 40 0
8
6 26 0
8
6 26
4 7 6 33 0 4 5 19 0 0 2 6
8 6 2 1 0 0 8 6 2 A= 4 11 7 = 0.5 1 0 0 8 6 =LU
28/29
My own faith has been strong for years. Back in 1985, I made a $100 bet with Peter Alfeld that a fast matrix inverse would be found by 1995. It wasn’t, so I paid him $100 and renewed the bet for another decade. It still hadn’t been found by 2005, so I paid him a further $100 and we renewed once more. One day I will win!—or my heirs?
我们经常使用记号O表示算法是“是多少阶的”。
因此高斯消元法为O(n3/3)的算法。复杂度(complexity)
是衡量算法性能优劣的主要指标。
26/29
例1. 在一台特定的计算机上, 估计用高斯消元法 求解5000*5000的线性方程组所需时间。
高斯消元法求解50*50的线性方程组所需时间为t1。
22/29
代码:
%elimination phase
for k = 1:n-1 %pivot row
for i= k+1:n %rows to be transformed
if A(i,k) ˜= 0 %multiplier
lambda = A(i,k)/A(k,k);
A(i,k+1:n) = A(i,k+1:n) - lambda*A(k,k+1:n);
回代阶段 0.1242 0.0147 0.0015 0.0001
高斯消元法由不相等的两部分组成:工作量相对多
的消元阶段和工作量相对少的回代阶段。对于大的n, 幂次较低的n相比而言可以忽略不计,最高次幂决定了 当n趋近于无穷时的极限形态。换而言之, 对大的n, 低 阶项对算法的执行时间的估计没有太大的影响, 当仅需 近似估计执行时间时可以忽略不计。
11/29
(Cramer法则) 设n阶矩阵A可逆,则线性方程组Ax b
有唯一解x ( x1, x2 , , xn )T ,其中
xj
det Aj det A
( j 1,2,
, n)
计算一个n阶行列式需要的乘法次数为(n -1)n。
以n 30为例,总的乘法次数大约为7.7e+33。
7.7e+33/365(天) / 24(小时) / 3600(秒) / 109 2.4e+18(年) 7.7e+33/365(天) / 24(小时) / 3600(秒)（/ 4.7 1015） 5.2e+10(年)
36
51
13
52
34
74
0
7
1.1
减肥所要求 每日营养量
33 45 3
1/29
2/29
3/29
4/29
5/29
6/29
大规模代数系统的求解是科学与工程计算 的基础和关键,其计算量也占有非常大的比重, 有的甚至占整个计算的80%以上！
Gene Golub (1932—2007)
美国科学院,工程院、艺术与科学学院院士 以及瑞典皇家工程科学院院士
503 / 50003
3 /3
t503 / t 50003
3 /3
t1 t2
t2 1003 t1
27/29
回顾:
高斯消元法是时间复杂度为O(n3)的算法, 最 终在有限步之内给出一个解。因此高斯消元法 叫做求解线性方程组的直接法。
理论上直接法在有限步之内给出精确解(当 然在一有限精度的计算机上执行时得到的结果 将只是近似的)。
x3 y3 / 2 3 x2 ( y2 6x3 ) / 8 1 x1 ( y1 6 x2 2 x3 ) / 8 2
16/29
Erdős 相信上帝手中有一本包含世间所有精妙 证明的天书。上帝相信这本书在 Gauss 手上。
17/29
高斯消元法
高斯消元法是求解线性方程组的经典方法之一。高 斯消元法由两个部分构成: 消元阶段和回代阶段 。消元 阶段的目的是将一般线性方程组转化成容易求解的上三 角线性方程组。
a11 a12 a1n x1 b1
a21
a22
a2n
x2
b2
an1
an2
ann
xn
bn
Ax=b
9/29
存在性
秩是定量刻画矩阵性质的指标, 相对于矩阵行列式更加精细。
增广矩阵:
a11 a12
[
A,
b]
a21
a22
an1 an2
a1n b1
a2n
b2
ann
y1
28
0.5 y1 + y2 40
0.5y1 0.5 y2 y3 33
y1 28 y2 (40 0.5 y1 )2 26 y3 33 0.5 y1 0.5 y2 6
然后求解上三角方程组Ux = y :
8 x1 6 x2 2 x3 y1 8x2 6x3 y2 2x3 y3
14/29
例 求解方程组Ax = b, 其中
8 6 2 28 A 4 11 7 , b 40
4 7 6 33
假设系数矩阵A存在如下LU分解:
1 0 0 8 6 2 A LU 0.5 1 0 0 8 6
0.5 0.5 1 0 0 2
15/29
解: 首先求解下三角线性方阵组Ly = b:
设3种食物每100克中蛋白质、碳水化合物和脂肪的含量 如下表; 表中还给出了美国流行的剑桥大学医学院的简捷营养处方。 现在的问题是:如果用这3种食物作为每天的主要食物,那么它们的 用量应各取多少才能全面准确地实现这个营养。
营养
蛋白质 碳水化合物 脂肪
每100g食物所含营养(g)
脱脂牛奶 大豆面粉
乳清
b(i)= b(i) - lambda*b(k);
end
end end
xk [bk (ak,k1 xk1 ak ,n xn )] / ak ,k (k n 1, 1)
%back-substitution phase
for k = n:-1:1
b(k) = (b(k) - A(k,k+1:n)*x(k+1:n))/A(k,k);
7/29
求解
a11 x1 a12 x2 a1n xn b1 a21x1a22 x2 a2n xn b2 an1 x1 an2 x2 ann xn bn
8/29
➢线性方程组的矩阵形式
a11 x1 a12 x2 a1n xn b1 a21x1a22 x2 a2n xn b2 an1 x1 an2 x2 ann xn bn
12/29
求解如下下三角线性方程组Lx=c:
l11 x1 c1 l21 x1 l22 x2 c2 l31 x1 l32 x2 l33 x3 c3
从第一个方程开始, 从前往后每次求解一个方程, 则下三 角方程组的解可以表示如下
x1 c1 / l11 x2 (c2 l21 x1 ) / l22 x3 (c3 l31 x1 l32 x2 ) / l33
例如
8x1 6x2 2x3 28 (a) 4x1 11x2 7 x3 40 (b)
4x1 7 x2 6x3 33 (c)
18/29
如何将线性方程组化为更容易求解的线性方程组且 不改变线性方程组的解？
➢ 交换两个方程; ➢ 某一方程乘以非零常数c; ➢ 把某一方程的c倍加于另一方程。
O((n-1)n!) (Cramer) O(n3) (Gauss) O(n2.81) (Strassen) O(n2.373) O(n2+eps)?
Ref: The smart money’s on numerical analysis, SIAM News, 2012
29/29
19/29
消元阶段:
经过第一轮消元:
8x1 6x2 2x3 28
(a)
8 x2 6 x3 26 (b)
4 x2 5 x3 19
(c)
注释: 8是主元(pivot), -0.5和0.5是乘子(multiplier)。
经过第二轮消元:
8x1 6x2 2x3 28 (a) 8x2 6x3 26 (b) 2x3 6 (c)
共计： 2n3/3+n2/2-7n/6
回代阶段:
除法次数 n次, 乘法次数 n(n-1)/2次,
加(减)法次数 n(n-1)/2次
共计： n2
25/29
n=10 n=100 n=1000 n=10000
消元阶段 0.8282 0.0621 0.0145 0.9782 0.0073 0.0001 0.9978 0.0007 0.0000 0.9998 0.0000 0.0000
4 7 6 0.5 0.5 1 0 0 2
Guass消元法和LU分解的关系？ 21/29
回代阶段:
从最后一个方程开始, 从后往前每次求解一个方程, 则 三角方程组的解可以表示如下
x3 6 / 2 3 x2 (26 6x3 ) / 8 1 x1 (28 6 x2 2 x3 ) / 8 2
13/29
类似的, 求解如下上三角线性方程组Ux=c:
u11 x1 u12 x2 u13 x3 c1 u22 x2 u23 x3 c2 u33 x3 c3
从最后一个方程开始, 从后往前每次求解一个方程, 则上 三角方程组的解可以表示如下
x3 c3 / u33 x2 (c2 u23 x3 ) / u22 x1 (c1 u13 x3 u12 x2 ) / u11
end
23/29
24/29
12 22 32
消元阶段:
n 2 n(n1)(2n1) 6
除法次数(n-1)+ (n-2)+···+1= n(n-1)/2,
乘法次数(n-1)*n+ (n-2)* (n-1) +···+1*2= (n3-n)/3,
加(减)法次数(n-1)*n+ (n-2)* (n-1) +···+1*2= (n3-n)/3。
bn
如果rank([A,b])=rank(A), 则线性方程组的有解。
10/29
唯一性 如果系数矩阵非奇异即|A| ≠ 0, 则线性方程组的解唯一。 如果系数矩阵奇异即|A| = 0,则线性方程组的有无穷多 解或无解。例如, 2x + y = 3和4x + 2y = 6有无穷多解; 而 2x + y = 3和4x + 2y = 0则无解。