2023-06-27《计算方法》- 陈丽娟 - 线性方程组的直接解法.md

合集 - 《计算方法》-陈丽娟(8)

1.2023-06-16 《计算方法》- 陈丽娟 - 绪论.md2023-06-162.2023-06-16《计算方法》- 陈丽娟 - 插值法（一）2023-06-16

3.2023-06-27《计算方法》- 陈丽娟 - 线性方程组的直接解法.md2023-06-27

4.2023-06-19《计算方法》- 陈丽娟 - 方程的近似解法（注解）2023-06-195.2023-06-18《计算方法》- 陈丽娟 - 方程的近似解法.md2023-06-186.2023-06-17《计算方法》- 陈丽娟 - 插值法（二）.md2023-06-177.2023-06-30《计算方法》- 陈丽娟 - 线性方程组的迭代解法.md2023-06-308.2023-06-28《计算方法》- 陈丽娟 - 向量和矩阵基础.md2023-06-28

收起

2023-06-27《计算方法》- 陈丽娟 - 线性方程组的直接解法

Matlab计算方法高斯消元法矩阵分解

线性方程组的解法这一课题我们在高等代数中已经了解过，对于一个非奇异方阵，通过求解或者克莱姆法则均可以直接得到方程的精确解，但是上述方法计算量很大，难以在实际中应用，因此引出了本章的内容。

首先，我们给出问题的形式以方便叙述：

其中

问题中已知，待求变量为.

1. 高斯消元法

高斯消元法是一个经典的算法，其基本思想是先把方程组转化为一个上三角方程组，然后逆次序逐一求出未知量。

针对, 首先记增广矩阵。

Step 1. 从第2行到第n行中消去的系数. 为了做到这一点，我们从高代中知道，一行加上另一行的倍不影响最终的解，因此可以将第一行的加到第行，, 即是
.

Step 2. 依次重复Step 1中的计算，可以得到最终的三角矩阵.

Step 3. 回代得到解.

高斯消元法的时间复杂度为:

1. 在消元的过程中，第步需要计算的除法显然有次，因有行和列，则有乘法和减法各次，因此乘法和除法的时间复杂度为
   
   
   
   减法和回代过程的时间复杂度显然比消元的过程小，因此高斯消元法的时间复杂度为, 另外由于没有需要额外存储的数据，其空间复杂度即矩阵的大小.

=我们介绍了高斯消元法的基本步骤，那么当满足什么条件的时候高斯消元法起作用呢？这里给出一个定理来说明这个问题。

定理1:
高斯消元法求解的充分必要条件是的顺序主子式全不为0.

这是由于在高斯消元的过程中总是要求.

虽然理论上上述高斯消元法可以得到问题的解，但是当使用计算机编程时，由于不可避免的误差和有效数字，在某些情况下上述方法不能得到目标解，书中给出了如下例子：

例子1:
在三位有效数字下，给定方程组

精确解为.
用上述消元法消去第二个式子的得到, .
而若使用第二个式子消去第一个式子中的, 则有.

因此为了减小误差，给出
列主元高斯消元法:
Step 1. 找到第一列中最大元素所在的行, 交换第1行和第行。 然后将第一行的加到第行，.

Step 2. 依次重复Step 1中的计算, 找第列中的最大元素, 可以得到最终的三角矩阵.

Step 3. 回代得到解.

附上代码：

1. %% 顺序、列主元高斯消元法 
   

   ---

2. % 输入系数矩阵A, 向量b, 是否使用列主元方法iscol 
   

   ---

3. % 返回求解向量x, 最终的上三角矩阵A1, 以及对应的向量b1 
   

   ---

4. function [x, A1, b1] = GaussLE(A, b, iscol) 
   

   ---

5. x = b; 
   

   ---

6. if iscol == 1 
   

   ---

7. % 主元高斯消元法 
   

   ---

8. AA = [A b]; 
   

   ---

9. n = length(b); 
   

   ---

10. for i = 1:(n-1) 
    

    ---

11. if AA(i,i) == 0 
    

    ---

12. % 主元素为0，不能进行高斯消元返回并报错 
    

    ---

13. x = 0; 
    

    ---

14. A1 = 0; 
    

    ---

15. b1 = 0; 
    

    ---

16. print("The element is 0."); 
    

    ---

17. break 
    

    ---

18. else 
    

    ---

19. % 找到最大的主元 
    

    ---

20. [M I] = max(AA(i:n,i)); 
    

    ---

21. % 交换位置 
    

    ---

22. I = I + i - 1; 
    

    ---

23. temp = AA(I,:); 
    

    ---

24. AA(I,:) = AA(i,:); 
    

    ---

25. AA(i,:) = temp; 
    

    ---

26. for j = (i+1):n 
    

    ---

27. mij = - AA(j,i) / AA(i,i); 
    

    ---

28. AA(j,:) = AA(j,:) + mij .* AA(i,:); 
    

    ---

29. end 
    

    ---

30. end 
    

    ---

31. end 
    

    ---

32. % 回代过程 
    

    ---

33. x(n) = AA(n, n+1) / AA(n, n); 
    

    ---

34. for i = flip(1:(n-1)) 
    

    ---

35. x(i) = (AA(i, n+1) - sum(AA(i, (i+1):n) .* x((i+1):n)')) / AA(i, i); 
    

    ---

36. end 
    

    ---

37. A1 = AA(:,1:n); 
    

    ---

38. b1 = AA(:,n+1); 
    

    ---

39. else 
    

    ---

40. % 顺序高斯消元法 
    

    ---

41. AA = [A b]; 
    

    ---

42. n = length(b); 
    

    ---

43. for i = 1:(n-1) 
    

    ---

44. if AA(i,i) == 0 
    

    ---

45. % 主元素为0，不能进行高斯消元返回并报错 
    

    ---

46. x = 0; 
    

    ---

47. A1 = 0; 
    

    ---

48. b1 = 0; 
    

    ---

49. print("The element is 0."); 
    

    ---

50. break 
    

    ---

51. else 
    

    ---

52. for j = (i+1):n 
    

    ---

53. mij = - AA(j,i) / AA(i,i); 
    

    ---

54. AA(j,:) = AA(j,:) + mij .* AA(i,:); 
    

    ---

55. end 
    

    ---

56. end 
    

    ---

57. end 
    

    ---

58. % 回代过程 
    

    ---

59. x(n) = AA(n, n+1) / AA(n, n); 
    

    ---

60. for i = flip(1:(n-1)) 
    

    ---

61. x(i) = (AA(i, n+1) - sum(AA(i, (i+1):n) .* x((i+1):n)')) / AA(i, i); 
    

    ---

62. end 
    

    ---

63. A1 = AA(:,1:n); 
    

    ---

64. b1 = AA(:,n+1); 
    

    ---

65. end 
    

    ---

66. end 
    

    ---

虽然可以将判断是否列主元消元置于循环当中以减少代码的重复度，但是这样需要增加n次判断，因此还是分开写着两类消元法。
然后我们给出一个例子：

1. format long 
   

   ---

2.  
   

   ---

3. A = rand(100,100)*1000; 
   

   ---

4. b = rand(100,1)*1000; 
   

   ---

5. % 顺序消元 
   

   ---

6. [x1, A1, b1] = GaussLE(A, b, 0); 
   

   ---

7. % 列主元消元 
   

   ---

8. [x2, A2, b2] = GaussLE(A, b, 1); 
   

   ---

9. % 内置求解器 
   

   ---

10. x3 = linsolve(A, b); 
    

    ---

11. plot(abs(x1-x3)) 
    

    ---

12. hold on 
    

    ---

13. plot(abs(x2-x3)) 
    

    ---

14. legend(["顺序消元" "列主元消元"]) 
    

    ---

和matlab内置求解器的误差为

这个不严谨的例子可以部分说明列主元是要比顺序消元更好。

2. 高斯消元法的矩阵形式和LU分解

由于在矩阵的某一行加上另一行的某倍实际上等价于矩阵左乘, 为基坐标, 为倍数, 为单位阵。 因此前述的一系列变换相当于矩阵左乘了一系列矩阵：, 同时由于均为下三角矩阵，因此也是下三角矩阵，则原问题可表述为

其中是消元后的系数（显然是上三角矩阵）。
利用三角矩阵的优良形式可以快速求得原问题的解：令, 先求解, 解得, 再求解.

那么如何通过分解得到呢？首先假设的对角线元素均为1, 则立即可得.

我们有

依次看，即对的第2行有

对的第3行有

对的第r行有

对的第r列有

上面我们说明了对于顺序主子式都大于0的矩阵，存在分解，使得是单位下三角矩阵，是上三角矩阵，下面我们说明这个分解是唯一的。设存在两个分解, 可得, 由于是上三角矩阵，是单位下三角矩阵，因此, 即, .

下面给出三角分解算法求解线性方程组的程序：

1. %% 求解线性方程组的LU分解 
   

   ---

2. % 输入系数矩阵A，向量b， 如果b为空，则只返回分解 
   

   ---

3. function [x, L, U] = LULE(A, b) 
   

   ---

4. %% 先对A进行LU分解 
   

   ---

5. x = 0; 
   

   ---

6. n = size(A,1); 
   

   ---

7. L = eye(n); 
   

   ---

8. U = zeros(n, n); 
   

   ---

9. U(1,:) = A(1,:); 
   

   ---

10. L(:,1) = A(:,1) / U(1,1); 
    

    ---

11. for i = 2:n 
    

    ---

12. for j = i:n 
    

    ---

13. U(i,j) = A(i,j) - L(i,1:(i-1)) * U(1:(i-1),j); 
    

    ---

14. L(j,i) = (A(j,i) - L(j,1:(i-1)) * U(1:(i-1),i)) / U(i,i); 
    

    ---

15. end 
    

    ---

16. end 
    

    ---

17. if length(b) == n 
    

    ---

18. %% 然后计算 y 和 x 
    

    ---

19. x = zeros(n,1); 
    

    ---

20. y = x; 
    

    ---

21. y(1) = b(1); 
    

    ---

22. for i = 2:n 
    

    ---

23. y(i) = b(i) - L(i,1:(i-1)) * y(1:(i-1)); 
    

    ---

24. end 
    

    ---

25. x(n) = y(n) / U(n,n); 
    

    ---

26. for i = flip(1:(n-1)) 
    

    ---

27. x(i) = (y(i) - U(i,(i+1):n) * x((i+1):n))/U(i,i); 
    

    ---

28. end 
    

    ---

29. end 
    

    ---

30. end 
    

    ---

测试程序

1. format long 
   

   ---

2.  
   

   ---

3. A = rand(100,100)*1000; 
   

   ---

4. b = rand(100,1)*1000; 
   

   ---

5. % 顺序消元 
   

   ---

6. [x1, A1, b1] = GaussLE(A, b, 0); 
   

   ---

7. % 列主元消元 
   

   ---

8. [x2, A2, b2] = GaussLE(A, b, 1); 
   

   ---

9. % 内置求解器 
   

   ---

10. x3 = linsolve(A, b); 
    

    ---

11.  
    

    ---

12. [x4, L, U] = LULE(A, b); 
    

    ---

13. plot(abs(x1-x3)) 
    

    ---

14. hold on 
    

    ---

15. plot(abs(x2-x3)) 
    

    ---

16. hold on 
    

    ---

17. plot(abs(x4-x3)) 
    

    ---

18. legend(["顺序消元" "列主元消元" "LU"]) 
    

    ---

测试结果：

enter description here

在上矩阵理论的时候还有接触过其他分解形式，但是距离久远以及在家没有教材因此不便叙述，等会学校再补充相关知识。