2023-06-27《计算方法》- 陈丽娟 - 线性方程组的直接解法
Matlab计算方法高斯消元法矩阵分解线性方程组的解法这一课题我们在高等代数中已经了解过，对于一个非奇异方阵，通过求解或者克莱姆法则均可以直接得到方程的精确解，但是上述方法计算量很大，难以在实际中应用，因此引出了本章的内容。
首先，我们给出问题的形式以方便叙述：
|
|
其中
问题中
已知，待求变量为
.
1. 高斯消元法
高斯消元法是一个经典的算法，其基本思想是先把方程组转化为一个上三角方程组，然后逆次序逐一求出未知量。
针对
, 首先记增广矩阵
。
Step 1. 从第2行到第n行中消去
的系数. 为了做到这一点，我们从高代中知道，一行加上另一行的
倍不影响最终的解，因此可以将第一行的
加到第
行，
, 即是
.
Step 2. 依次重复Step 1中的计算，可以得到最终的三角矩阵.
Step 3. 回代得到解.
高斯消元法的时间复杂度为
:
- 在消元的过程中，第
步需要计算的除法
显然有
次，因有
行和
列，则有乘法和减法各
次，因此乘法和除法的时间复杂度为
减法和回代过程的时间复杂度显然比消元的过程小，因此高斯消元法的时间复杂度为
, 另外由于没有需要额外存储的数据，其空间复杂度即矩阵的大小
.
=我们介绍了高斯消元法的基本步骤，那么当
满足什么条件的时候高斯消元法起作用呢？这里给出一个定理来说明这个问题。
定理1:
高斯消元法求解
的充分必要条件是
的顺序主子式全不为0.
这是由于在高斯消元的过程中总是要求
.
虽然理论上上述高斯消元法可以得到问题
的解，但是当使用计算机编程时，由于不可避免的误差和有效数字，在某些情况下上述方法不能得到目标解，书中给出了如下例子：
例子1:
在三位有效数字下，给定方程组
精确解为
.用上述消元法消去第二个式子的
得到
, 
.而若使用第二个式子消去第一个式子中的
, 则有
.
因此为了减小误差，给出
列主元高斯消元法:
Step 1. 找到第一列中最大元素所在的行
, 交换第1行和第
行。 然后将第一行的
加到第
行，
.
Step 2. 依次重复Step 1中的计算, 找第
列中的最大元素, 可以得到最终的三角矩阵.
Step 3. 回代得到解.
附上代码：
- %% 顺序、列主元高斯消元法
- % 输入系数矩阵A, 向量b, 是否使用列主元方法iscol
- % 返回求解向量x, 最终的上三角矩阵A1, 以及对应的向量b1
- function [x, A1, b1] = GaussLE(A, b, iscol)
- x = b;
- if iscol == 1
- % 主元高斯消元法
- AA = [A b];
- n = length(b);
- for i = 1:(n-1)
- if AA(i,i) == 0
- % 主元素为0，不能进行高斯消元返回并报错
- x = 0;
- A1 = 0;
- b1 = 0;
- print("The element is 0.");
- break
- else
- % 找到最大的主元
- [M I] = max(AA(i:n,i));
- % 交换位置
- I = I + i - 1;
- temp = AA(I,:);
- AA(I,:) = AA(i,:);
- AA(i,:) = temp;
- for j = (i+1):n
- mij = - AA(j,i) / AA(i,i);
- AA(j,:) = AA(j,:) + mij .* AA(i,:);
- end
- end
- end
- % 回代过程
- x(n) = AA(n, n+1) / AA(n, n);
- for i = flip(1:(n-1))
- x(i) = (AA(i, n+1) - sum(AA(i, (i+1):n) .* x((i+1):n)')) / AA(i, i);
- end
- A1 = AA(:,1:n);
- b1 = AA(:,n+1);
- else
- % 顺序高斯消元法
- AA = [A b];
- n = length(b);
- for i = 1:(n-1)
- if AA(i,i) == 0
- % 主元素为0，不能进行高斯消元返回并报错
- x = 0;
- A1 = 0;
- b1 = 0;
- print("The element is 0.");
- break
- else
- for j = (i+1):n
- mij = - AA(j,i) / AA(i,i);
- AA(j,:) = AA(j,:) + mij .* AA(i,:);
- end
- end
- end
- % 回代过程
- x(n) = AA(n, n+1) / AA(n, n);
- for i = flip(1:(n-1))
- x(i) = (AA(i, n+1) - sum(AA(i, (i+1):n) .* x((i+1):n)')) / AA(i, i);
- end
- A1 = AA(:,1:n);
- b1 = AA(:,n+1);
- end
- end
虽然可以将判断是否列主元消元置于循环当中以减少代码的重复度，但是这样需要增加n次判断，因此还是分开写着两类消元法。
然后我们给出一个例子：
- format long
-
- A = rand(100,100)*1000;
- b = rand(100,1)*1000;
- % 顺序消元
- [x1, A1, b1] = GaussLE(A, b, 0);
- % 列主元消元
- [x2, A2, b2] = GaussLE(A, b, 1);
- % 内置求解器
- x3 = linsolve(A, b);
- plot(abs(x1-x3))
- hold on
- plot(abs(x2-x3))
- legend(["顺序消元" "列主元消元"])
和matlab内置求解器的误差为
这个不严谨的例子可以部分说明列主元是要比顺序消元更好。
2. 高斯消元法的矩阵形式和LU分解
由于在矩阵的某一行加上另一行的某倍实际上等价于矩阵左乘
, 
为基坐标, 
为倍数, 
为单位阵。 因此前述的一系列变换相当于矩阵
左乘了一系列矩阵：
, 同时由于
均为下三角矩阵，因此
也是下三角矩阵，则原问题
可表述为
其中
是消元后的系数（显然是上三角矩阵）。利用三角矩阵的优良形式可以快速求得原问题的解：令
, 先求解
, 解得
, 再求解
.
那么如何通过
分解得到
呢？首先假设
的对角线元素均为1, 则立即可得
.
我们有
依次看，即对
的第2行有
对
的第3行有
对
的第r行有
对
的第r列有
上面我们说明了对于顺序主子式都大于0的矩阵
，存在
分解，使得
是单位下三角矩阵，
是上三角矩阵，下面我们说明这个分解是唯一的。设存在两个分解
, 可得
, 由于
是上三角矩阵，
是单位下三角矩阵，因此
, 即
, 
.
下面给出三角分解算法求解线性方程组的程序：
- %% 求解线性方程组的LU分解
- % 输入系数矩阵A，向量b， 如果b为空，则只返回分解
- function [x, L, U] = LULE(A, b)
- %% 先对A进行LU分解
- x = 0;
- n = size(A,1);
- L = eye(n);
- U = zeros(n, n);
- U(1,:) = A(1,:);
- L(:,1) = A(:,1) / U(1,1);
- for i = 2:n
- for j = i:n
- U(i,j) = A(i,j) - L(i,1:(i-1)) * U(1:(i-1),j);
- L(j,i) = (A(j,i) - L(j,1:(i-1)) * U(1:(i-1),i)) / U(i,i);
- end
- end
- if length(b) == n
- %% 然后计算 y 和 x
- x = zeros(n,1);
- y = x;
- y(1) = b(1);
- for i = 2:n
- y(i) = b(i) - L(i,1:(i-1)) * y(1:(i-1));
- end
- x(n) = y(n) / U(n,n);
- for i = flip(1:(n-1))
- x(i) = (y(i) - U(i,(i+1):n) * x((i+1):n))/U(i,i);
- end
- end
- end
测试程序
- format long
-
- A = rand(100,100)*1000;
- b = rand(100,1)*1000;
- % 顺序消元
- [x1, A1, b1] = GaussLE(A, b, 0);
- % 列主元消元
- [x2, A2, b2] = GaussLE(A, b, 1);
- % 内置求解器
- x3 = linsolve(A, b);
-
- [x4, L, U] = LULE(A, b);
- plot(abs(x1-x3))
- hold on
- plot(abs(x2-x3))
- hold on
- plot(abs(x4-x3))
- legend(["顺序消元" "列主元消元" "LU"])
测试结果：
enter description here
在上矩阵理论的时候还有接触过其他分解形式，但是距离久远以及在家没有教材因此不便叙述，等会学校再补充相关知识。
标签：AA,md,27,06,矩阵,元法,消元,end,高斯消 From： https://www.cnblogs.com/NEFPHYS/p/17509669.html