§ 2 线性方程组的误差分析
/* Error Analysis for Linear system of Equations */
求解 时,A 和 的误差对解 有何影响?bxA b? x?
设 A 精确,有误差,得到的解为,即b? b xx
bbxxA )(
bAx 1 |||||||||||| 1 bAx
绝对误差放大因子
|||||||||||||||| xAxAb又
||||
||||
||||
1
b
A
x
||||
||||||||||||
||||
|||| 1
b
bAA
x
x
相对误差放大因子
§ 2 Error Analysis for,bxA
设 精确,A有误差,得到的解为,即b? A? xx
bxxAA ))((
bxxAxxA )()(
)(1 xxAAx
||||
||||
||||||||
||||||||
||||
||||
1
1
A
A
AA
AA
xx
x
bxAAxAA )()(
xAxAA )(
xAxAAIA )( 1
xAAAAIx 111 )(
Wait a minute …
Who said that ( I + A?1?A )
is invertible?
(只要?A充分小,使得
)1|||||||||||| 11 AAAA
||||
||||
||||||||1
||||
||||
||||||||
||||||||1
||||||||
||||
||||
1
1
1
1
A
A
AA
A
A
AA
AA
AA
x
x
是关键的误差放大因子,称为
A的 条件数,记为 cond (A),
越 则 A 越病态,
难得准确解。
|||||||| 1 AA
大
§ 2 Error Analysis for,bxA
注,?cond (A) 的具体大小与 || · || 的取法有关,但相对大小一致。
cond (A) 取决于 A,与解题方法无关。
||||
||||
||||
||||
||||||||)(1
)(
||||
||||
b
b
A
A
AAAc o nd
Ac o nd
x
x
常用条件数有:
cond (A)1 cond (A)? cond (A)2 )(/)(
m inm a x AAAA TT
特别地,若 A 对称,则
||m i n
||m a x)(
2?
Ac o nd
条件数的性质:
A可逆,则 cond (A)p? 1;
A可逆, R 则 cond (? A) = cond (A) ;
A正交,则 cond (A)2=1;
A可逆,R正交,则 cond (RA)2 = cond (AR)2 = cond (A)2 。
§ 2 Error Analysis for,bxA
精确解 为,
1
1
x?例,
97.1
99.1,
98.099.0
99.01 bA?
计算 cond (A)2 。?
1 0 0 0 09 9 0 0
9 9 0 09 8 0 0A?1 =
解,考察 A 的特征根
0)d e t ( AI?
0 0 0 0 5 0 5 0 4.0
9 8 0 0 5 0 5 0 4.1
2
1
2
1
2)(?
Ac o n d 39206 >> 1
测试病态程度:
给 一个扰动b?
3
4
10106.0
1097.0b,其相对误差为
%01.0105 1 3.0|||| |||| 4
2
2
b
b 此时 精确解 为
0203.1
3*x?
0 2 0 3.2
2* xxx
2
2
||||
||||
x
x 2.0102 > 200%
§ 2 Error Analysis for,bxA
例,Hilbert 阵
12
1
1
11
3
1
2
1
1
2
11
nnn
n
nH
cond (H2)?= 27 cond (H3) 748
cond (H6)?= 2.9? 106 cond (Hn) as n
注,一般判断矩阵是否病态,并不计算 A?1,而由经验得出。
行列式很大或很小(如某些行、列近似相关);
元素间相差大数量级,且无规则;
主元消去过程中出现小主元;
特征值相差大数量级。
§ 2 Error Analysis for,bxA
近似解的误差估计及改善:
设 的近似解为,则一般有bxA *x? 0* xAbr
b
r
x
xx ||||
||||
||*|| cond (A)
误差上限
改善方法:
Step 1,近似解 bxA ;1x?
Step 2,;11 xAbr
Step 3,;
111 drdA
Step 4,;
112 dxx
若 可被精确解出,则有就是精确解了。
1d
bAxAbAxx 11112 )(
2x?
经验表明,若 A 不是非常病态(例如,),
则如此迭代可达到机器精度;但若 A 病态,则此算法也不能改进。
1)(Ac o n d?
HW,p.66
#2,#4,#5
§ 3 Jacobi 法和 Gauss - Seidel 法
/* Jacobi & Gauss-Seidel Iterative Methods */
Jacobi Iterative Method
nnnnnn
nn
nn
bxaxaxa
bxaxaxa
bxaxaxa
...
............
...
...
2211
22222121
11212111
nnnnn
nn
n
nn
nn
bxaxa
a
x
bxaxa
a
x
bxaxa
a
x
1111
22121
22
2
11212
11
1
...
1
............
...
1
...
1
0?iia
写成 矩阵形式,
A =
L
UD bxULxD bxULDbxA )( )(
bDxULDx 11 )(
B f?
Jacobi 迭代阵 bDxULDx kk 1)(1)1( )(
§ 3 Jacobi & Gauss-Seidel Iterative Methods
Algorithm,Jacobi Iterative Method
Solve given an initial approximation,
Input,the number of equations and unknowns n; the matrix entries a[ ][ ];
the entries b[ ]; the initial approximation X0[ ]; tolerance TOL;
maximum number of iterations Nmax.
Output,approximate solution X[ ] or a message of failure.
Step 1 Set k = 1;
Step 2 While ( k? Nmax) do steps 3-6
Step 3 For i = 1,…,n
Set ; /* compute xk */
Step 4 If then Output (X[ ]);
STOP; /* successful */
Step 5 For i = 1,…,n Set X 0[ ] = X [ ]; /* update X0 */
Step 6 Set k ++;
Step 7 Output (Maximum number of iterations exceeded);
STOP,/* unsuccessful */
bxA )0(x?
ii
n
j
ij
jiji
i a
Xab
X
1
)0(
T O LXXXX iini |0|m a x||0|| 1
What if aii = 0?
迭代过程中,A的元素不改变,故可以 事先调整 好 A 使得
aii? 0,否则 A不可逆 。
必须等 X(k)完全计算好了才能计算 X(k+1),因此需要 两组向量 存储。
A bit wasteful,
isn’t it?
§ 3 Jacobi & Gauss-Seidel Iterative Methods
Gauss - Seidel Iterative Method
)(1 1)(1)(414)(313)(212
11
)1(
1 bxaxaxaxaax
k
nn
kkkk
)(1 2)(2)(424)(323)1(121
22
)1(
2 bxaxaxaxaax
k
nn
kkkk
)(1 3)(3)(434)1(232)1(131
33
)1(
3 bxaxaxaxaax
k
nn
kkkk
)(1 )1( 11)1(33)1(22)1(11)1( nknnnknknkn
nn
k
n bxaxaxaxaax
… … … …
只存一组向量即可。
写成 矩阵形式,bDxUxLDx kkk 1)()1(1)1( )(
bxUxLD kk )()1()(
bLDxULDx kk 1)(1)1( )()(
B f?Gauss-Seidel
迭代阵
A mathematician about his colleague,
" He made a lot of mistakes,but he made them
in a good direction,I tried to copy this,but I found
out that it is very difficult to make good mistakes,"
§ 3 Jacobi & Gauss-Seidel Iterative Methods
注,二种方法都存在 收敛性问题 。
有例子表明,Gauss-Seidel法收敛时,Jacobi法可能不收敛;而 Jacobi法收敛时,Gauss-Seidel法也可能不收敛。
p.76 #2 给出了例子。
收敛性分析将在下节课讨论。
§ 3 Jacobi & Gauss-Seidel Iterative Methods
Lab 07,Gauss-Seidel Method
Use the Gauss-Seidel method to solve a given n× n linear system
with an initial approximation and a given tolerance TOL.
Input
There are several sets of inputs,For each set:
The 1st line contains an integer 100? n? 0 which is the size of a
matrix,n =?1 signals the end of file.
The following n lines contain the
augmented matrix in the following format:
The numbers are separated by spaces and new lines.
The last line of each test case contains a real number TOL,which
is the tolerance for || ·||? norm,and an integer N? 0 which is the
maximum number of iteration.
bxA 0)0(x
nnnn
n
n
baa
baa
baa
.,,
.,,.,,.,,.,,
.,,
.,,
1
2221
1111
§ 3 Jacobi & Gauss-Seidel Iterative Methods
Output /*?represents a space */
Each entry of the solution is to be printed as in the C fprintf:
fprintf(outfile,"%12.8f\n",x );
If the matrix has a zero column,print the message,Matrix?has?a?zero?
column.No?unique?solution?exists.\n”,
If the method fails to give a solution after N iterations,print the message
,Maximum?number?of?iterations?exceeded.\n”,
If there is an entry of that is out of the range,print the message
,No?convergence.\n”,
The outputs of two test cases must be seperated by a blank line.
Sample Input
(?represents a space)
3
10?–1?0?9
–1?10?–2?7
0?–4?10?6
0.000000001?1000
3
3?–1?0?1
3?6?0?0
3?3?0?4
0.000000001?1000
–1
)(kx? ]2,2[ 127127?
Sample Output
(?represents a space)
1.00000000
1.00000000
1.00000000
Matrix?has?a?zero?column.No?unique
solution?exists.
/* Error Analysis for Linear system of Equations */
求解 时,A 和 的误差对解 有何影响?bxA b? x?
设 A 精确,有误差,得到的解为,即b? b xx
bbxxA )(
bAx 1 |||||||||||| 1 bAx
绝对误差放大因子
|||||||||||||||| xAxAb又
||||
||||
||||
1
b
A
x
||||
||||||||||||
||||
|||| 1
b
bAA
x
x
相对误差放大因子
§ 2 Error Analysis for,bxA
设 精确,A有误差,得到的解为,即b? A? xx
bxxAA ))((
bxxAxxA )()(
)(1 xxAAx
||||
||||
||||||||
||||||||
||||
||||
1
1
A
A
AA
AA
xx
x
bxAAxAA )()(
xAxAA )(
xAxAAIA )( 1
xAAAAIx 111 )(
Wait a minute …
Who said that ( I + A?1?A )
is invertible?
(只要?A充分小,使得
)1|||||||||||| 11 AAAA
||||
||||
||||||||1
||||
||||
||||||||
||||||||1
||||||||
||||
||||
1
1
1
1
A
A
AA
A
A
AA
AA
AA
x
x
是关键的误差放大因子,称为
A的 条件数,记为 cond (A),
越 则 A 越病态,
难得准确解。
|||||||| 1 AA
大
§ 2 Error Analysis for,bxA
注,?cond (A) 的具体大小与 || · || 的取法有关,但相对大小一致。
cond (A) 取决于 A,与解题方法无关。
||||
||||
||||
||||
||||||||)(1
)(
||||
||||
b
b
A
A
AAAc o nd
Ac o nd
x
x
常用条件数有:
cond (A)1 cond (A)? cond (A)2 )(/)(
m inm a x AAAA TT
特别地,若 A 对称,则
||m i n
||m a x)(
2?
Ac o nd
条件数的性质:
A可逆,则 cond (A)p? 1;
A可逆, R 则 cond (? A) = cond (A) ;
A正交,则 cond (A)2=1;
A可逆,R正交,则 cond (RA)2 = cond (AR)2 = cond (A)2 。
§ 2 Error Analysis for,bxA
精确解 为,
1
1
x?例,
97.1
99.1,
98.099.0
99.01 bA?
计算 cond (A)2 。?
1 0 0 0 09 9 0 0
9 9 0 09 8 0 0A?1 =
解,考察 A 的特征根
0)d e t ( AI?
0 0 0 0 5 0 5 0 4.0
9 8 0 0 5 0 5 0 4.1
2
1
2
1
2)(?
Ac o n d 39206 >> 1
测试病态程度:
给 一个扰动b?
3
4
10106.0
1097.0b,其相对误差为
%01.0105 1 3.0|||| |||| 4
2
2
b
b 此时 精确解 为
0203.1
3*x?
0 2 0 3.2
2* xxx
2
2
||||
||||
x
x 2.0102 > 200%
§ 2 Error Analysis for,bxA
例,Hilbert 阵
12
1
1
11
3
1
2
1
1
2
11
nnn
n
nH
cond (H2)?= 27 cond (H3) 748
cond (H6)?= 2.9? 106 cond (Hn) as n
注,一般判断矩阵是否病态,并不计算 A?1,而由经验得出。
行列式很大或很小(如某些行、列近似相关);
元素间相差大数量级,且无规则;
主元消去过程中出现小主元;
特征值相差大数量级。
§ 2 Error Analysis for,bxA
近似解的误差估计及改善:
设 的近似解为,则一般有bxA *x? 0* xAbr
b
r
x
xx ||||
||||
||*|| cond (A)
误差上限
改善方法:
Step 1,近似解 bxA ;1x?
Step 2,;11 xAbr
Step 3,;
111 drdA
Step 4,;
112 dxx
若 可被精确解出,则有就是精确解了。
1d
bAxAbAxx 11112 )(
2x?
经验表明,若 A 不是非常病态(例如,),
则如此迭代可达到机器精度;但若 A 病态,则此算法也不能改进。
1)(Ac o n d?
HW,p.66
#2,#4,#5
§ 3 Jacobi 法和 Gauss - Seidel 法
/* Jacobi & Gauss-Seidel Iterative Methods */
Jacobi Iterative Method
nnnnnn
nn
nn
bxaxaxa
bxaxaxa
bxaxaxa
...
............
...
...
2211
22222121
11212111
nnnnn
nn
n
nn
nn
bxaxa
a
x
bxaxa
a
x
bxaxa
a
x
1111
22121
22
2
11212
11
1
...
1
............
...
1
...
1
0?iia
写成 矩阵形式,
A =
L
UD bxULxD bxULDbxA )( )(
bDxULDx 11 )(
B f?
Jacobi 迭代阵 bDxULDx kk 1)(1)1( )(
§ 3 Jacobi & Gauss-Seidel Iterative Methods
Algorithm,Jacobi Iterative Method
Solve given an initial approximation,
Input,the number of equations and unknowns n; the matrix entries a[ ][ ];
the entries b[ ]; the initial approximation X0[ ]; tolerance TOL;
maximum number of iterations Nmax.
Output,approximate solution X[ ] or a message of failure.
Step 1 Set k = 1;
Step 2 While ( k? Nmax) do steps 3-6
Step 3 For i = 1,…,n
Set ; /* compute xk */
Step 4 If then Output (X[ ]);
STOP; /* successful */
Step 5 For i = 1,…,n Set X 0[ ] = X [ ]; /* update X0 */
Step 6 Set k ++;
Step 7 Output (Maximum number of iterations exceeded);
STOP,/* unsuccessful */
bxA )0(x?
ii
n
j
ij
jiji
i a
Xab
X
1
)0(
T O LXXXX iini |0|m a x||0|| 1
What if aii = 0?
迭代过程中,A的元素不改变,故可以 事先调整 好 A 使得
aii? 0,否则 A不可逆 。
必须等 X(k)完全计算好了才能计算 X(k+1),因此需要 两组向量 存储。
A bit wasteful,
isn’t it?
§ 3 Jacobi & Gauss-Seidel Iterative Methods
Gauss - Seidel Iterative Method
)(1 1)(1)(414)(313)(212
11
)1(
1 bxaxaxaxaax
k
nn
kkkk
)(1 2)(2)(424)(323)1(121
22
)1(
2 bxaxaxaxaax
k
nn
kkkk
)(1 3)(3)(434)1(232)1(131
33
)1(
3 bxaxaxaxaax
k
nn
kkkk
)(1 )1( 11)1(33)1(22)1(11)1( nknnnknknkn
nn
k
n bxaxaxaxaax
… … … …
只存一组向量即可。
写成 矩阵形式,bDxUxLDx kkk 1)()1(1)1( )(
bxUxLD kk )()1()(
bLDxULDx kk 1)(1)1( )()(
B f?Gauss-Seidel
迭代阵
A mathematician about his colleague,
" He made a lot of mistakes,but he made them
in a good direction,I tried to copy this,but I found
out that it is very difficult to make good mistakes,"
§ 3 Jacobi & Gauss-Seidel Iterative Methods
注,二种方法都存在 收敛性问题 。
有例子表明,Gauss-Seidel法收敛时,Jacobi法可能不收敛;而 Jacobi法收敛时,Gauss-Seidel法也可能不收敛。
p.76 #2 给出了例子。
收敛性分析将在下节课讨论。
§ 3 Jacobi & Gauss-Seidel Iterative Methods
Lab 07,Gauss-Seidel Method
Use the Gauss-Seidel method to solve a given n× n linear system
with an initial approximation and a given tolerance TOL.
Input
There are several sets of inputs,For each set:
The 1st line contains an integer 100? n? 0 which is the size of a
matrix,n =?1 signals the end of file.
The following n lines contain the
augmented matrix in the following format:
The numbers are separated by spaces and new lines.
The last line of each test case contains a real number TOL,which
is the tolerance for || ·||? norm,and an integer N? 0 which is the
maximum number of iteration.
bxA 0)0(x
nnnn
n
n
baa
baa
baa
.,,
.,,.,,.,,.,,
.,,
.,,
1
2221
1111
§ 3 Jacobi & Gauss-Seidel Iterative Methods
Output /*?represents a space */
Each entry of the solution is to be printed as in the C fprintf:
fprintf(outfile,"%12.8f\n",x );
If the matrix has a zero column,print the message,Matrix?has?a?zero?
column.No?unique?solution?exists.\n”,
If the method fails to give a solution after N iterations,print the message
,Maximum?number?of?iterations?exceeded.\n”,
If there is an entry of that is out of the range,print the message
,No?convergence.\n”,
The outputs of two test cases must be seperated by a blank line.
Sample Input
(?represents a space)
3
10?–1?0?9
–1?10?–2?7
0?–4?10?6
0.000000001?1000
3
3?–1?0?1
3?6?0?0
3?3?0?4
0.000000001?1000
–1
)(kx? ]2,2[ 127127?
Sample Output
(?represents a space)
1.00000000
1.00000000
1.00000000
Matrix?has?a?zero?column.No?unique
solution?exists.