本文详细介绍了使用Matlab求解线性方程组的各种方法，包括直接法（Gauss消元法、列主元消元法、Gauss-Jordan消元法、LU分解法、Cholesky分解法、LDL'分解法）和迭代法（Richardson迭代法、Jacobi迭代法、Gauss-Seidel迭代法、超松驰迭代法、SSOR迭代法、两步迭代法、最速下降法、共轭梯度法）。通过实例展示了每种方法的实现代码和求解过程。
线性方程组求解
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 amax<1e-10
return;
end
if r>k % 交换两行
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 amax<1e-10
return;
end
if r>k % 交换两行
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<1e-10
return;
end
L(k,k)=sqrt(delta);
for i=k+1:n
L(i,k)=A(i,k);
for j=1:k-1
L(i,k)=L(i,k)-L(i,j)*L(k,j);
end
L(i,k)=L(i,k)/L(k,k);
end
end
function x=Chol_Solve(A,b)
% 利用对称正定矩阵之 Cholesky 分解求解线性方程组 Ax=b
n=length(b);
l=Cholesky(A);
x=ones(1,n);
y=ones(1,n);
for i=1:n
z=0;
for k=1:i-1
z=z+l(i,k)*y(k);
end
y(i)=(b(i)-z)/l(i,i);
end
for i=n:-1:1
z=0;
for k=i+1:n
z=z+l(k,i)*x(k);
end
x(i)=(y(i)-z)/l(i,i);
end
Example:
>> 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))<1e-10
return;
end
for i=k+1:n
T(i,k)=A(i,k);
for j=1:k-1
T(i,k)=T(i,k)-T(i,j)*L(k,j);
end
L(i,k)=T(i,k)/d(k);
end
end
D=diag(d);
function x=LDL_Solve(A,b)
% 利用对对称正定矩阵 A 进行 LDL' 分解法求解线性方程组 Ax=b
n=length(b);
[l,d]=LDL_Factor(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)/d(n,n);
for i=n-1:-1:1
z=0;
for k=i+1:n
z=z+l(k,i)*x(k);
end
x(i)=y(i)/d(i,i)-z;
end
Example:
>> 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<3
error
return
elseif nargin ==5
M = varargin{1};
end
D=diag(diag(A)); % 求 A 的对角矩阵
L=-tril(A,-1); % 求 A 的下三角阵
U=-triu(A,1); % 求 A 的上三角阵
B=D/(L+U);
f=D/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]=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<3
error
return;
end
D=diag(diag(A)); % 求 A 的对角矩阵
L=-tril(A,-1); % 求 A 的下三角阵
U=-triu(A,1); % 求 A 的上三角阵
G=(D-L)/U;
f=(D-L)/b;
x=G*x0+f;
n=1; % 迭代次数
while norm(x-x0)>=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<=0 || 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<=0 || 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<3
error
return
elseif nargin ==5
M = varargin{1};
end
D=diag(diag(A)); % 求 A 的对角矩阵
L=-tril(A,-1); % 求 A 的下三角阵
U=-triu(A,1); % 求 A 的上三角阵
B1=(D-L)/U;
B2=(D-U)/L;
f1=(D-L)/b;
f2=(D-U)/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]=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 |
一般矩阵 |
|
TRANSA |
指出求解的方程是 Ax=b 还是 A’x=b ， A’ 为之共轭转置 |
Example:
>> A=[4 -1 2;-1 5 0;2 0 6];b=[1 -1 2]';
>> optt.SYM=true;optt.POSDEF=true;optt.TRANSA=false;
>> x=linsolve(A,b,optt)
4. 超定方程组的解法
利用伪逆求解：
X=pinv(A)*b
左除求解：
x=A/b
最小二乘法求解：
X=lsqnonneg(A,b)
最小二乘法求解：
A’Ax=A’b è x=A’*A/A’*b
Example:
>> A=[1 -2 3;4 1 0;7 1 6;9 5 8];b=[1 0 1 0]';
>> x=pinv(A)*b
x =
0.0903
-0.3248
0.1048
>> x=A/b
x =
0.0903
-0.3248
0.1048
>> x=lsqnonneg(A,b)
x =
0
0
0.0826
>> x=A'*A/A'*b
x =
0.0903
-0.3248
0.1048
5. 有无穷组解的线性方程组的解法
齐次线性方程组的通解：
特解与通解：
1． 求 Ax=b 的一个特解
2． 求 Ax=0 的通解
3． 将特解与通解组合成最终解
Example:
>> A=[4 8 -6 2;1 3 9 5;5 11 3 7;3 5 -15 -3];
>> b=[1 2 3 -1]';
>> d=[A b];
>> s=rref(d) % 采用增广矩阵阖求解
s =
1.0000 0 -22.5000 -8.5000 -3.2500
0 1.0000 10.5000 4.5000 1.7500
0 0 0 0 0
0 0 0 0 0
è 特解：（ -3.25 1.75 0 0 ） ’
基础解系有两个基向量： (-22.5 10.5 1 0)’ (-8.5 4.5 0 1)’
扫一扫
2694
到【灌水乐园】发言
