用矩阵分解法求解线性方程组

1 1 1 2 L = P2 L−1 P2 L−1 = 1 2 2 1 0 7 − 1 0 1 −1 7 2
数值分析
数值分析
P为 排 列 阵 ， 在 计 算 机 中 用 向 量 表 示 例 P = (1 2 3 4)T , P2 = (3 4 1 2)T , Ax = b, PA = LU , LUx = Pb = f f ( i ) = b ( P ( i )) f (1) = b( P (1)) = b(3), f (2) = b( P (2)) = b(4) f (3) = b( P (3)) = b(1), f (4) = b( P (4)) = b( 2 ) P1 = (3 2 1 4)T , P = (3 4 1 2)T PAx = Pb,
2 解 Ux = −1 7 2 0 0 1 -1 x1 1 x1 13 1 0 x 2 2 x2 2 ,得 2 = 13 x = − 13 ,得 x 0 −1 3 7 3 7 x x4 −1 0 11 4 − 11 7 7 0
数值分析
数值分析
lupdsv.m %功能：调用列主元三角分解函数 [LU,p]=lupd(A) % 求解线性方程组Ax=b。 。 求解线性方程组
%解法：PA=LU, Ax=b←→PAx=Pb 解法： 解法 % % LUx=Pb, Ly=f=Pb, y=Ux f(i)=b(p(i))
%输入：方阵A，右端项 （行或列向量均可） 输入：方阵 ，右端项b（行或列向量均可） 输入 %输出：解x（行向量） 输出： 输出 （行向量）
数值分析
数值分析
L2 P2 L1 P1 A = U , 令P = P2 P1
A = P1−1 L−1 P2−1 L−1U = P1 L−1 P2 L−1U 1 2 1 2
PA = P2 P1 P1 L−1 P2 L−1U = P2 L−1 P2 L−1U ， L = P2 L−1 P2 L−1 1 2 1 2 1 2 PA = LU
j =1 k −1
1 l 21 l 31 ⋮ ⋮ l n1
1 l 32 ⋮ ⋮ 1 ⋱ 1
l n 2 ⋯ ⋯ l n , n −1
y1 b1 y2 b2 = ⋮ ⋮ 1 yn bn
0 0 P1 = 1 0
−1 0 0 −1 0 3 0 ↔1 1
2 1 1 1
0 −2 b= −1 6
1 1 2 1
数值分析
0 1 0 1 0 0 0 0 0 0 0 1
2 −1 0 −1 0 0 P1 A = 0 1 −1 1 3 0
2 L1 P1 A =
2 P2 L1 P1 A =
−1 −1 1 7
−1 7 2 1
0 0 −1 0
0
2
2
1 3 2 2 1 2
1 2 −1 3 0 − 1 2 2 0 1 2 −1 7 1 0 2 2 =U L 2 P2 L1 P1 A = 0 − 1 13 7 11 0 0 7
这种方法没有选主元问题，计算稳定， 这种方法没有选主元问题，计算稳定，但计算量略 1 3 n ，正交分解法的计算量 大，消元法的计算量大约为 3 2 3 大约为 n 。 3
数值分析
数值分析
(1) A = Q R q 11 q 12 q 21 q 22 Q = ⋯ ⋯ q n1 q n 2 ( 2 )计 算 b = Q T b bk =
六、利用矩阵奇异值分解A = U ΣV 求解Ax = b
T
七、 三对角方程组的解法
数值分析
数值分析
一、利用三角分解A = LU 求解Ax = b 设已有 A = LU 代入原方程 Ax = b 得 LY = b LUx = b ⇔ Ux = Y
第一步 : 求解下三角方程组LY = b,向前回代求出 Y = ( y1 , y2 ,⋯, yn )T y1 = b1 yk = bk − ∑ lkj y j (k = 2, 3,⋯, n)
数值分析
数值分析
第二步 : 求解上三角方程组 Ux = Y ,向后回代求出 x x n = y n unn xk = ( yk −
j = k +1
∑u
n
kj
x j ) ukk y1 y 2 ⋮ yn
( k = n − 1, n − 2, ⋯ ,1) u11 u12 ⋯ u1 n x1 u22 ⋯ u2 n x 2 = ⋱ ⋮ ⋮ unn x n
数值分析
数值分析
function x=lupdsv(A,b) n=length(b); [LU,p]=lupd(A); y(1)=b(p(1)); for i=2:n y(i)=b(p(i))-LU(i,1:i-1)*y(1:i-1)'; end x(n)=y(n)/LU(n,n); for i=(n-1):-1:1 x(i)=(y(i)-LU(i,i+1:n)*x(i+1:n)')/LU(i,i); end
数值分析
数值分析
function x=lupqdsv(A,b) n=length(b); [LU,p,q]=lupqd(A); y(1)=b(p(1)); for i=2:n y(i)=b(p(i))-LU(i,1:i-1)*y(1:i-1)'; end z(n)=y(n)/LU(n,n);x(q(n))=z(n); for i=(n-1):-1:1 z(i)=(y(i)-LU(i,i+1:n)*z(i+1:n)')/LU(i,i); x(q(i))=z(i); end
数值分析
数值分析
−1 1 y1 1 y1 −1 13 1 y 2 2 y 6 2 = = Pb = f , 得 2 = 解LY = 2 y3 − 13 1 y3 0 0 7 7 y −2 y4 − 11 − 1 − 1 0 1 4 7 2 7
数值分析
数值分析
二、 用列主元的三角分解PA = LU 求解Ax = b LY = Pb Ax = b ⇔ PAx = Pb ⇒ LUx = Pb ⇔ Ux = Y 用列主元三角分解求解Ax Ax= 例：用列主元三角分解求解Ax=b
0 −1 其中A = 2 1 i1 解： i1 = 3,
数值分析
数值分析
三、 用全主元的三角分解PAQT = LU 求解Ax = b Ax = b ⇔ PAQ T (Qx ) = Pb ⇒ LU (Qx ) = Pb
lupqdsv.m %功能：调用全主元三角分解函数[LU,p,q]=lupqd(A) 功能：调用全主元三角分解函数 功能 % 求解线性方程组Ax=b。 。 求解线性方程组
数值分析
第三节 用矩阵分解法求解线性方程组
一、利用三角分解A = LU 求解Ax = b
二、 用列主元的三角分解PA = LU 求解Ax = b
三、 用全主元的三角分解PAQT = LU 求解Ax = b
四、利用Cholesky分解A = LLT 求解Ax = b
五、利用正交分解A = QR求解 Ax = b
%解法：PAQ-1=LU, Ax=b←→(PAQ-1)(Qx)=Pb 解法： 解法 % % % LU(Qx)=Pb, z=Qx, y=Uz Ly=f=Pb, f(i)=b(p(i)) Uz=y, z=Qx , x(q(i))=z(i).
%输入：方阵A，右端项 （行或列向量均可） 输入：方阵 ，右端项b（行或列向量均可） 输入 %输出：解x（行向量） 输出： 输出 （行向量）
⋯ ⋯ ⋯ ⋯
q1 n r11 q2n ,R = ⋯ q nn
r12 r22
⋯ ⋯
r1 n ⋯ ⋯ rnn
∑q
i =1
n
ik
b i , k = 1, 2, ⋯ , n
( 3 )回 代 求 解 R x = b x n = b n rnn n , ( k = n − 1, ⋯ , 2, 1) x k = ( b n − ∑ rkj x j ) rkk j= k +1
数值分析
数值分析
四、利用Cholesky分解A = LLT 求解Ax = b A ∈ R n×n是对称正定的矩阵，有Cholesky分解式 A = LLT 代入原方程 LLT x = b, 可分成两步
(1)求 解 下 三 角 方 程 组 (1) 求 LY = b y i = ( bi − ∑ l ik y k ) / l ii , ( i = 1, 2, ⋯ , n )
a11 a 21 A = 0 ⋮ 0 a12 a 22 ⋱ ⋱ ⋯ 0 a 23 ⋱ a n − 1, n − 2 0 ⋯ ⋱ ⋱ a n −1, n − 1 a n , n −1 0 ⋮ 0 a n −1, n an,n A称为三对 称为三对 角矩阵。 角矩阵。
k =1 i −1
(2)求 (2)求解上三角方程组 LT x = Y xi = ( yi −
k = i +1
∑l
n
ki
xk ) / lii , ( i = n, n − 1,⋯ ,1)
数值分析
数值分析
五、利用正交分解A = QR求解 Ax = b Ax = b ⇔ QRx = b ⇔ Rx = QT b，于是可分成三步 ）正交分解 （1）正交分解 A = QR （2）矩阵乘法 b = QT b ）矩阵乘法 ）回代求解上三角方程组 （3）回代求解上三角方程组 Rx = b