用matlab实现高斯列主元消去法解线性方程及LU分解

Posion

1. <font face="Times New Roman">
   
2. 用matlab实现高斯列主元消去法解线性方程及LU分解
   
3. 
   
4. function x=gaussLinearEquation(A,b)
   
5. 
   
6. %高斯法解线性方程Ax=b
   
7. 
   
8. disp('原方程为AX=b：')
   
9. 
   
10. A
    
11. 
    
12. b
    
13. 
    
14. disp('------------------------')
    
15. 
    
16. n=length(b);
    
17. 
    
18. eps=10^-2;
    
19. 
    
20. for k=1:n-1
    
21. 
    
22.     %找列主元
    
23. 
    
24.     [mainElement,index]=max(abs(A(k:n,k)));
    
25. 
    
26.     index=index+k-1;%index在A(k:n,k)中的行号转换为在A中的行号
    
27. 
    
28.     if abs(mainElement)<eps
    
29. 
    
30.         disp('列元素太小！！');
    
31. 
    
32.         break;
    
33. 
    
34.     elseif index>k
    
35. 
    
36.         %列主元所在行不是当前行，将当前行与列主元所在行交换
    
37. 
    
38.         temp=A(k,:);
    
39. 
    
40.         A(k,:)=A(index,:);
    
41. 
    
42.         A(index,:)=temp;
    
43. 
    
44.     end
    
45. 
    
46.     %消元
    
47. 
    
48.     for i=k+1:n
    
49. 
    
50.         m(i,k)=A(i,k)/A(k,k);%A(k,k)将A(i,k)消为0所乘系数
    
51. 
    
52.         A(i,k:n)=A(i,k:n)-m(i,k)*A(k,k:n);%第i行消元处理
    
53. 
    
54.         b(i)=b(i)-m(i,k)*b(k);%还有b也需要处理！！      
    
55. 
    
56.     end
    
57. 
    
58. end
    
59. 
    
60. disp('消元后所得到的上三角阵是')
    
61. 
    
62. A
    
63. 
    
64. %回代
    
65. 
    
66. b(n)=b(n)/A(n,n);
    
67. 
    
68. for i=n-1:-1:1
    
69. 
    
70.     %sum(A(i,i+1:n).*b(i+1:n)')表示已知
    
71. 
    
72.     b(i)=(b(i)-sum(A(i,i+1:n).*b(i+1:n)'))/A(i,i);
    
73. 
    
74. end
    
75. 
    
76. clear x;
    
77. 
    
78. x=b;
    
79. 
    
80. disp('AX=b的解x是')
    
81. 
    
82. x
    
83. 
    
84. 用法：
    
85. 
    
86. 在控制台输入：
    
87. 
    
88. A=[1.003 0.333 1.504 -0.333;
    
89. 
    
90.    -2.011 1.455 0.506 2.956;
    
91. 
    
92.    4.329 -1.952 0.006 2.087;
    
93. 
    
94.    5.113 -4.004 3.332 -1.112];
    
95. 
    
96. b=[ 3.005,5.407,0.136,3.772 ]';
    
97. 
    
98. 执行gaussLinearEquation(A,b);即可得到结果。
    
99. 
    
100. 使用matlab实现矩阵的LU分解
     
101. function [L,U]=myLU(A)
     
102. %实现对矩阵A的LU分解，L为下三角矩阵
     
103. A
     
104. [n,n]=size(A);
     
105. L=zeros(n,n);
     
106. U=zeros(n,n);
     
107. for i=1:n
     
108.     L(i,i)=1;
     
109. end
     
110. for k=1:n
     
111.     for j=k:n
     
112.         U(k,j)=A(k,j)-sum(L(k,1:k-1).*U(1:k-1,j)');
     
113.       
     
114.     end
     
115.     for i=k+1:n
     
116.         L(i,k)=(A(i,k)-sum(L(i,1:k-1).*U(1:k-1,k)'))/U(k,k);      
     
117.     end
     
118.    
     
119. end
     
120. 用法，在控制台输入
     
121. A=[1 2 3 -4;-3 -4 -12 13;2 10 0 -3;4 14 9 -13];
     
122. 然后执行[l,u]=myLU(A);
     
123. 这样得到l和u，可以通过l*u与A比较来验证LU分解的正确性。
     
124. 
     
125. 注意：
     
126. [L1,U1]=lu(x)
     
127. [L2,U2,P]=lu(x)
     
128. [L3,U3,P,Q] = lu(X)
     
129. MATLAB中[L1,U1]=lu(x)的结果：L是下三角的置换矩阵即L1=p*L2，U是上三角阵。Matlab中LU分解采用高斯消元法，结果是不唯一的，只要[L1,U1]=lu(x)满足L1*U1=x, [L2,U2,P]=lu(x)满足L2*U2=p*x，[L3,U3,P,Q] = lu(X)满足 L3*U3= P*X*Q就行了。</font>

复制代码

转自：http://blog.sina.com.cn/s/blog_49c02a8c0100yt1x.html