线性方程组求解Matlab程序

线性方程组求解 1 直接法 Gauss 消元法 function x DelGauss a b Gauss 消去法 n m size a nb length b det 1 存储行列式值 x zeros n 1 for k 1 n 1 for i k 1 n if a k k 0 return end m a i k a k k for j k 1 n a i j a i j m a k j end b i b i m b k end det det a k k end det det a n n for k n 1 1 回代 for j k 1 n b k b k a k j x j end x k b k a k k end Example A 1 0170 0 0092 0 0095 0 0092 0 9903 0 0136 0 0095 0 0136 0 9898 b 1 0 1 x DelGauss A b x 0 9739 0 0047 1 0010 列主元 Gauss 消去法 function x detGauss a b Gauss 列主元消去法 n m size a nb length b det 1 存储行列式值 x zeros n 1 for k 1 n 1 amax 0 选主元 for i k n if abs a i k amax amax abs a i k r i end end if amaxk 交换两行 for j k n z a k j a k j a r j a r j z end z b k b k b r b r z det det end for i k 1 n 进行消元 m a i k a k k for j k 1 n a i j a i j m a k j end b i b i m b k end det det a k k end det det a n n for k n 1 1 回代 for j k 1 n b k b k a k j x j end x k b k a k k end Example x detGauss A b x 0 9739 0 0047 1 0010 Gauss Jordan 消去法 function x GaussJacobi a b Gauss Jacobi 消去法 n m size a nb length b x zeros n 1 for k 1 n amax 0 选主元 for i k n if abs a i k amax amax abs a i k r i end end if amaxk 交换两行 for j k n z a k j a k j a r j a r j z end z b k b k b r b r z end 进行消元 b k b k a k k for j k 1 n a k j a k j a k k end for i 1 n if i k for j k 1 n a i j a i j a i k a k j end b i b i a i k b k end end end for i 1 n x i b i end Example x GaussJacobi A b x 0 9739 0 0047 1 0010 LU 分解法 function l u lu a LU 分解 n length a l eye n u zeros n for i 1 n u 1 i a 1 i end for i 2 n l i 1 a i 1 u 1 1 end for r 2 n for i r n uu 0 for k 1 r 1 uu uu l r k u k i end u r i a r i uu end for i r 1 n ll 0 for k 1 r 1 ll ll l i k u k r end l i r a i r ll u r r end End function x lusolv a b LU 分解求解线性方程组 aX b if length a length b error Error in inputing return end n length a l u lu a y 1 b 1 for i 2 n z 0 for k 1 i 1 z z l i k y k end y i b i z end x n y n u n n for i n 1 1 1 z 0 for k i 1 n z z u i k x k end x i y i z u i i end Example x lusolv A b x 0 9739 0 0047 1 0010 对称正定矩阵之 Cholesky 分解法 function L Cholesky A 对对称正定矩阵 A 进行 Cholesky 分解 n length A L zeros n for k 1 n delta A k k for j 1 k 1 delta delta L k j 2 end if delta a 9 36 30 36 192 180 30 180 180 b 1 1 1 x Chol Solve a b x 1 8333 1 0833 0 7833 对称正定矩阵之 LDL 分解法 function L D LDL Factor A 对称正定矩阵 A 进行 LDL 分解 n length A L eye n D zeros n d zeros 1 n T zeros n for k 1 n d k A k k for j 1 k 1 d k d k L k j T k j end if abs d k x LDL Solve a b x 1 8333 1 0833 0 7833 2 迭代法 Richardson 迭代法 function x n richason A b x0 eps M Richardson 法求解线性方程组 Ax b 方程组系数矩阵 A 方程组之常数向量 b 迭代初始向量 X0 e 解的精度控制 eps 迭代步数控制 M 返回值线性方程组的解 x 返回值迭代步数 n if nargin 3 eps 1 0e 6 M 200 elseif nargin 4 M 200 end I eye size A x1 x0 x I A x0 b n 1 while norm x x1 eps x1 x x I A x1 b n n 1 if n M disp Warning 迭代次数太多 现在退出 return end end Example A 1 0170 0 0092 0 0095 0 0092 0 9903 0 0136 0 0095 0 0136 0 9898 b 1 0 1 x0 0 0 0 x n richason A b x0 x 0 9739 0 0047 1 0010 n 5 Jacobi 迭代法 function x n jacobi A b x0 eps varargin if nargin 3 eps 1 0e 6 M 200 elseif nargin eps x0 x x B x0 f n n 1 if n M disp Warning 迭代次数太多 可能不收敛 return end end Example x n Jacobi A b x0 x 0 9739 0 0047 1 0010 n 5 Gauss Seidel 迭代法 function x n gauseidel A b x0 eps M if nargin 3 eps 1 0e 6 M 200 elseif nargin 4 M 200 elseif nargin eps x0 x x G x0 f n n 1 if n M disp Warning 迭代次数太多 可能不收敛 return end end Example x n gauseidel A b x0 x 0 9739 0 0047 1 0010 n 4 超松驰迭代法 function x n SOR A b x0 w eps M if nargin 4 eps 1 0e 6 M 200 elseif nargin 4 error return elseif nargin 5 M 200 end if w 2 error return end D diag diag A 求 A 的对角矩阵 L tril A 1 求 A 的下三角阵 U triu A 1 求 A 的上三角阵 B inv D L w 1 w D w U f w inv D L w b x B x0 f n 1 迭代次数 while norm x x0 eps x0 x x B x0 f n n 1 if n M disp Warning 迭代次数太多 可能不收敛 return end end Example x n SOR A b x0 1 x 0 9739 0 0047 1 0010 n 4 对称逐次超松驰迭代法 function x n SSOR A b x0 w eps M if nargin 4 eps 1 0e 6 M 200 elseif nargin 4 error return elseif nargin 5 M 200 end if w 2 error return end D diag diag A 求 A 的对角矩阵 L tril A 1 求 A 的下三角阵 U triu A 1 求 A 的上三角阵 B1 inv D L w 1 w D w U B2 inv D U w 1 w D w L f1 w inv D L w b f2 w inv D U w b x12 B1 x0 f1 x B2 x12 f2 n 1 迭代次数 while norm x x0 eps x0 x x12 B1 x0 f1 x B2 x12 f2 n n 1 if n M disp Warning 迭代次数太多 可能不收敛 return end end Example x n SSOR A b x0 1 x 0 9739 0 0047 1 0010 n 3 两步迭代法 function x n twostep A b x0 eps varargin if nargin 3 eps 1 0e 6 M 200 elseif nargin eps x0 x x12 B1 x0 f1 x B2 x12 f2 n n 1 if n M disp Warning 迭代次数太多 可能不收敛 return end end Example x n twostep A b x0 x 0 9739 0 0047 1 0010 n 3 最速下降法 function x n fastdown A b x0 eps if nargin 3 eps 1 0e 6 end r b A x0 d dot r r dot A r r x x0 d r n 1 while norm x x0 eps x0 x r b A x0 d dot r r dot A r r x x0 d r n n 1 end Example x n fastdown A b x0 x 0 9739 0 0047 1 0010 n 5 共轭梯度法 function x n conjgrad A b x0 if nargin 3 eps 1 0e 6 end r1 b A x0 p1 r1 d dot r1 r1 dot p1 A p1 x x0 d p1 r2 r1 d A p1 f dot r2 r2 dot r1 r1 p2 r2 f p1 n 1 for i 1 rank A 1 x0 x p1 p2 r1 r2 d dot r1 r1 dot p1 A p1 x x0 d p1 r2 r1 d A p1 f dot r2 r2 dot r1 r1 p2 r2 f p1 n n 1 end d dot r2 r2 dot p2 A p2 x x d p2 n n 1 Example x n conjgrad A b x0 x 0 9739 0 0047 1 0010 n 4 预处理的共轭梯度法 当 AX B 为病态方程组时 共轭梯度法收敛很慢 预处理技术是在用共轭梯度法求解之前 对系数矩阵做一些变换后再求解 Example A 25 300 1050 1400 630 300 4800 18900 26880 12600 1050 18900 79380 117600 56700 1400 26880 117600 179200 88200 630 12600 56700 88200 44100 b 5 3 1 0 2 x0 0 0 0 0 0 M pascal 5 预处理矩阵 x flag re it pcg A b 1 e 8 1000 M M x0 flag 0 表示在指定迭代次数之内按要求精度收敛 re 表示相对误差 it 表示迭代次数 x 5 7667 2 9167 1 9310 1 4333 1 1349 flag 0 re 5 7305e 012 it 10 其他迭代法 函数 说明 x symmlq A b 线性方程组的 LQ 解法 x bicg A b 线性方程组的双共轭梯度法 x bicgstab A b 线性方程组的稳定双共轭梯度法 x lsqr A b 线性方程组的共轭梯度 LSQR 解法 x gmres A b 线性方程组的广义最小残差解法 x minres A b 线性方程组的最小残差解法 x qmr A b 线性方程组的准最小残差解法 3 特殊解法 解三对角线性方程组之追赶法 function x followup A b n rank A for i 1 n if A i i 0 disp Error 对角有元素为 0 return end end d ones n 1 a ones n 1 1 c ones n 1 for i 1 n 1 a i 1 A i 1 i c i 1 A i i 1 d i 1 A i i end d n 1 A n n for i 2 n d i 1 d i 1 a i 1 1 d i 1 1 c i 1 1 b i 1 b i 1 a i 1 1 d i 1 1 b i 1 1 end x n 1 b n 1 d n 1 for i n 1 1 1 x i 1 b i 1 c i 1 x i 1 1 d i 1 end Example A 2 5 1 0 0 0 0 0 1 1 5 1 0 0 0 0 0 1 0 0 5 1 0 0 0 0 0 1 0 0 5 1 0 0 0 0 0 1 0 1 5 1 0 0 0 0 0 1 0 2 5 b ones 6 1 x followup A b x 0 4615 0 1538 0 7692 0 7692 0 1538 0 4615 快速求解法 通用求解线性方程组的函数 x linsolve A b options 其意义为快速求解方程组 Ax b 其中 A 之结构由决定 内容如下表 options 说明 LT 上三角阵 UT 下三角阵 UHESS 上三角 Hessenberg SYM 实对称矩阵 POSDEF 正定矩阵 RECT 一般矩阵 TRAN