用matlab解线性方程组
2008-04-12 17:00
一。高斯消去法
1.顺序高斯消去法
直接编写命令文件
a=[]
d=[]'
[n,n]=size(a);
c=n+1
a(:,c)=d;
for k=1:n-1
a(k+1:n, k:c)=a(k+1:n, k:c)-(a(k+1:n,k)/ a(k,k))*a(k, k:c); %消去
end
x=[0 0 0 0]' %回带
x(n)=a(n,c)/a(n,n);
for g=n-1:-1:1
x(g)=(a(g,c)-a(g,g+1:n)*x(g+1:n))/a(g,g)
end
2.列主高斯消去法
*由于“[r,m]=max(abs(a(k:n,k)))”返回的行是“k:n,k”内的第几行，所以要通过修正来把m 改成真正的行的值。该程序只是演示程序，真正机器计算不需要算主元素所在列以下各行应为零的值。
直接编写命令文件
a=[]
d=[] '
[n,n]=size(a);
c=n+1
a(:,c)=d; %(增广)
for k=1:n-1
[r,m]=max(abs(a(k:n,k))); %选主
m=m+k-1; %(修正操作行的值) if(a(m,k)~=0)
if(m~=k)
a([k m],:)=a([m k],:); %换行
end
a(k+1:n, k:c)=a(k+1:n, k:c)-(a(k+1:n,k)/ a(k,k))*a(k, k:c); %消去
end
end
x=[0 0 0 0]' %回带
x(n)=a(n,c)/a(n,n);
for g=n-1:-1:1
x(g)=(a(g,c)-a(g,g+1:n)*x(g+1:n))/a(g,g)
end
3．分别用顺序高斯消去法和列主高斯消去法解方程组a*x=d，并比较结果
a=[0 1 2 3;9 11 23 34;62.5 23.4 15.5 17.2;192.01 124 25.1 59.3] d=[1;1;1;1]
顺序高斯消去法：提示“Warning: Divide by zero.” x =NaN NaN NaN NaN
列主高斯消去法：x =-1.2460 2.8796 5.5206 -4.3069
由此可见列主高斯消去法可以解决顺序高斯消去法所不能解决的问题。
4. 将上述矩阵中的“2”改为2.05，“34”改为“34.6”，“15.5”改为“15.57”，“124”改为“124.7”再用列主高斯消去法计算，与上述结果比较。
x =-0.8081 1.8226 3.5568 -2.7047
很明显虽然系数矩阵只有很小的变化但结果的变化却很大，所以系数矩阵是病态的。 这是因为系数矩阵的条件数很大：cond(a)=2.0695e+003
二。迭代法
J迭代公式
a=[5 2 1;-1 4 2;2 -3 10]
d=[-12;20;3]
x=[0;0;0]; %初始向量
stop=1.0e-4 %迭代的精度
L=-tril(a,-1)
U=-triu(a,1)
D=inv(diag(diag(a)))
X=D*(L+U)*x+D*d; % J迭代公式
n=1;
while norm(X-x,inf)>=stop % 时迭代中止否则继续
x=X;
X=D*(L+U)*x+D*d;
n=n+1;
end
X
n
G-S迭代公式
a=[5 2 1;-1 4 2;2 -3 10]
x=[0;0;0]; %初始向量
d=[-12;20;3]
stop=1.0e-4
L=-tril(a,-1)
U=-triu(a,1)
D=(diag(diag(a)))
X=inv(D-L)*U*x+inv(D-L)*d; % G-S迭代公式
n=1;
while norm(X-x,inf)>= stop % 时迭代中止否则继续
x=X;
X=inv(D-L)*U*x+inv(D-L)*d;
n=n+1;
end
X
n
SOR迭代公式
a=[5 2 1;-1 4 2;2 -3 10]
x=[0;0;0]; %初始向量
d=[-12;20;3]
stop=1.0e-4
w=1.44; %松弛因子
L=-tril(a,-1)
U=-triu(a,1)
D=(diag(diag(a)))
X=inv(D-w*L)*((1-w)*D+w*U)*x+w*inv(D-w*L)*d % SOR迭代公式
n=1;
while norm(X-x,inf)>=stop % 时迭代中止否则继续
x=X;
X=inv(D-w*L)*((1-w)*D+w*U)*x+w*inv(D-w*L)*d;
n=n+1;
end
X
n
3. a*x=d a=[5 2 1;-1 4 2;2 -3 10] d=[-12;20;3]
(1)考察用J、G-S及SOR解此方程组的收敛性.
(2)分别用雅可比迭代法、高斯-赛德尔迭代法及逐次超松弛迭代法解此方程组。要求时 迭代中止，观察收敛速度。
(1) J迭代矩阵的谱半径：max(abs(eig(D*(L+U))))= 0.506
G-S迭代矩阵的谱半径：max(abs(eig(inv(D-L)*U)))= 0.2000
SOR迭代矩阵的谱半径：max(abs(eig(inv(D-w*L)*((1-w)*D+w*U))))=0.9756
所以用J、G-S及SOR均收敛（且有 ）。
(2)
J： X =-4.0000 3.0000 2.0000 n =18
G-S： X =-4.0000 3.0000 2.0000 n =8
SOR( )： X =-4.0000 3.0000 2.0000 n =482
可见高斯-赛德尔迭代法要比雅可比迭代公式的收敛速度快。同时，如果超松弛迭代法的 选取不当不但不会加速收敛还会对收敛速度造成很恶劣的结果。
藏
线性方程组求解
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
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
return;
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
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
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))
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;
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
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
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
error
return
elseif nargin ==5
M = 200;
end
if(w=2)
error;
return;
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
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
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;
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
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);
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 =
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.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)’
本文来自CSDN博客，转载请标http://blog.csdn.net/Guassfans/archive/2008/04/26/2330862.aspx
出处 明：
用matlab解线性方程组
2008-04-12 17:00
一。高斯消去法
1.顺序高斯消去法
直接编写命令文件
a=[]
d=[]'
[n,n]=size(a);
c=n+1
a(:,c)=d;
for k=1:n-1
a(k+1:n, k:c)=a(k+1:n, k:c)-(a(k+1:n,k)/ a(k,k))*a(k, k:c); %消去
end
x=[0 0 0 0]' %回带
x(n)=a(n,c)/a(n,n);
for g=n-1:-1:1
x(g)=(a(g,c)-a(g,g+1:n)*x(g+1:n))/a(g,g)
end
2.列主高斯消去法
*由于“[r,m]=max(abs(a(k:n,k)))”返回的行是“k:n,k”内的第几行，所以要通过修正来把m 改成真正的行的值。该程序只是演示程序，真正机器计算不需要算主元素所在列以下各行应为零的值。
直接编写命令文件
a=[]
d=[] '
[n,n]=size(a);
c=n+1
a(:,c)=d; %(增广)
for k=1:n-1
[r,m]=max(abs(a(k:n,k))); %选主
m=m+k-1; %(修正操作行的值) if(a(m,k)~=0)
if(m~=k)
a([k m],:)=a([m k],:); %换行
end
a(k+1:n, k:c)=a(k+1:n, k:c)-(a(k+1:n,k)/ a(k,k))*a(k, k:c); %消去
end
end
x=[0 0 0 0]' %回带
x(n)=a(n,c)/a(n,n);
for g=n-1:-1:1
x(g)=(a(g,c)-a(g,g+1:n)*x(g+1:n))/a(g,g)
end
3．分别用顺序高斯消去法和列主高斯消去法解方程组a*x=d，并比较结果
a=[0 1 2 3;9 11 23 34;62.5 23.4 15.5 17.2;192.01 124 25.1 59.3] d=[1;1;1;1]
顺序高斯消去法：提示“Warning: Divide by zero.” x =NaN NaN NaN NaN
列主高斯消去法：x =-1.2460 2.8796 5.5206 -4.3069
由此可见列主高斯消去法可以解决顺序高斯消去法所不能解决的问题。
4. 将上述矩阵中的“2”改为2.05，“34”改为“34.6”，“15.5”改为“15.57”，“124”改为“124.7”再用列主高斯消去法计算，与上述结果比较。
x =-0.8081 1.8226 3.5568 -2.7047
很明显虽然系数矩阵只有很小的变化但结果的变化却很大，所以系数矩阵是病态的。 这是因为系数矩阵的条件数很大：cond(a)=2.0695e+003
二。迭代法
J迭代公式
a=[5 2 1;-1 4 2;2 -3 10]
d=[-12;20;3]
x=[0;0;0]; %初始向量
stop=1.0e-4 %迭代的精度
L=-tril(a,-1)
U=-triu(a,1)
D=inv(diag(diag(a)))
X=D*(L+U)*x+D*d; % J迭代公式
n=1;
while norm(X-x,inf)>=stop % 时迭代中止否则继续
x=X;
X=D*(L+U)*x+D*d;
n=n+1;
end
X
n
G-S迭代公式
a=[5 2 1;-1 4 2;2 -3 10]
x=[0;0;0]; %初始向量
d=[-12;20;3]
stop=1.0e-4
L=-tril(a,-1)
U=-triu(a,1)
D=(diag(diag(a)))
X=inv(D-L)*U*x+inv(D-L)*d; % G-S迭代公式
n=1;
while norm(X-x,inf)>= stop % 时迭代中止否则继续
x=X;
X=inv(D-L)*U*x+inv(D-L)*d;
n=n+1;
end
X
n
SOR迭代公式
a=[5 2 1;-1 4 2;2 -3 10]
x=[0;0;0]; %初始向量
d=[-12;20;3]
stop=1.0e-4
w=1.44; %松弛因子
L=-tril(a,-1)
U=-triu(a,1)
D=(diag(diag(a)))
X=inv(D-w*L)*((1-w)*D+w*U)*x+w*inv(D-w*L)*d % SOR迭代公式
n=1;
while norm(X-x,inf)>=stop % 时迭代中止否则继续
x=X;
X=inv(D-w*L)*((1-w)*D+w*U)*x+w*inv(D-w*L)*d;
n=n+1;
end
X
n
3. a*x=d a=[5 2 1;-1 4 2;2 -3 10] d=[-12;20;3]
(1)考察用J、G-S及SOR解此方程组的收敛性.
(2)分别用雅可比迭代法、高斯-赛德尔迭代法及逐次超松弛迭代法解此方程组。要求时 迭代中止，观察收敛速度。
(1) J迭代矩阵的谱半径：max(abs(eig(D*(L+U))))= 0.506
G-S迭代矩阵的谱半径：max(abs(eig(inv(D-L)*U)))= 0.2000
SOR迭代矩阵的谱半径：max(abs(eig(inv(D-w*L)*((1-w)*D+w*U))))=0.9756
所以用J、G-S及SOR均收敛（且有 ）。
(2)
J： X =-4.0000 3.0000 2.0000 n =18
G-S： X =-4.0000 3.0000 2.0000 n =8
SOR( )： X =-4.0000 3.0000 2.0000 n =482
可见高斯-赛德尔迭代法要比雅可比迭代公式的收敛速度快。同时，如果超松弛迭代法的 选取不当不但不会加速收敛还会对收敛速度造成很恶劣的结果。
藏
线性方程组求解
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
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
return;
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
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
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))
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;
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
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
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
error
return
elseif nargin ==5
M = 200;
end
if(w=2)
error;
return;
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
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
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;
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
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);
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 =
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.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)’
本文来自CSDN博客，转载请标http://blog.csdn.net/Guassfans/archive/2008/04/26/2330862.aspx
出处 明：