线性方程组求解Matlab程序精选.docx
《线性方程组求解Matlab程序精选.docx》由会员分享，可在线阅读，更多相关《线性方程组求解Matlab程序精选.docx（35页珍藏版）》请在冰点文库上搜索。
线性方程组求解Matlab程序精选
线性方程组求解
1.直接法
Gauss消元法：
functionx=DelGauss（a,b）
%Gauss消去法
[n,m]=size（a）;
nb=length（b）;
det=1;%存储行列式值
x=zeros（n,1）;
fork=1:
n-1
fori=k+1:
n
ifa（k,k）==0
return
end
m=a（i,k）/a（k,k）;
forj=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）;
fork=n:
-1:
1%回代
forj=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.00920.0095;-0.00920.99030.0136;0.00950.01360.9898];
>>b=[101]';
>>x=DelGauss（A,b）
x=
0.9739
-0.0047
1.0010
列主元Gauss消去法：
functionx=detGauss（a,b）
%Gauss列主元消去法
[n,m]=size（a）;
nb=length（b）;
det=1;%存储行列式值
x=zeros（n,1）;
fork=1:
n-1
amax=0;%选主元
fori=k:
n
ifabs（a（i,k））>amax
amax=abs（a（i,k））;r=i;
end
end
ifamax<1e-10
return;
end
ifr>k%交换两行
forj=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
fori=k+1:
n%进行消元
m=a（i,k）/a（k,k）;
forj=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）;
fork=n:
-1:
1%回代
forj=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消去法：
functionx=GaussJacobi（a,b）
%Gauss-Jacobi消去法
[n,m]=size（a）;
nb=length（b）;
x=zeros（n,1）;
fork=1:
n
amax=0;%选主元
fori=k:
n
ifabs（a（i,k））>amax
amax=abs（a（i,k））;r=i;
end
end
ifamax<1e-10
return;
end
ifr>k%交换两行
forj=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）;
forj=k+1:
n
a（k,j）=a（k,j）/a（k,k）;
end
fori=1:
n
ifi~=k
forj=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
fori=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）;
fori=1:
n
u（1,i）=a（1,i）;
end
fori=2:
n
l（i,1）=a（i,1）/u（1,1）;
end
forr=2:
n
%%%%
fori=r:
n
uu=0;
fork=1:
r-1
uu=uu+l（r,k）*u（k,i）;
end
u（r,i）=a（r,i）-uu;
end
%%%%
fori=r+1:
n
ll=0;
fork=1:
r-1
ll=ll+l（i,k）*u（k,r）;
end
l（i,r）=（a（i,r）-ll）/u（r,r）;
end
%%%%
End
functionx=lusolv（a,b）
%LU分解求解线性方程组aX=b
iflength（a）~=length（b）
error（'Errorininputing!
'）
return;
end
n=length（a）;
[l,u]=lu（a）;
y
（1）=b
（1）;
fori=2:
n
z=0;
fork=1:
i-1
z=z+l（i,k）*y（k）;
end
y（i）=b（i）-z;
end
x（n）=y（n）/u（n,n）;
fori=n-1:
-1:
1
z=0;
fork=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.00471.0010
对称正定矩阵之Cholesky分解法：
functionL=Cholesky（A）
%对对称正定矩阵A进行Cholesky分解
n=length（A）;
L=zeros（n）;
fork=1:
n
delta=A（k,k）;
forj=1:
k-1
delta=delta-L（k,j）^2;
end
ifdelta<1e-10
return;
end
L（k,k）=sqrt（delta）;
fori=k+1:
n
L（i,k）=A（i,k）;
forj=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
functionx=Chol_Solve（A,b）
%利用对称正定矩阵之Cholesky分解求解线性方程组Ax=b
n=length（b）;
l=Cholesky（A）;
x=ones（1,n）;
y=ones（1,n）;
fori=1:
n
z=0;
fork=1:
i-1
z=z+l（i,k）*y（k）;
end
y（i）=（b（i）-z）/l（i,i）;
end
fori=n:
-1:
1
z=0;
fork=i+1:
n
z=z+l（k,i）*x（k）;
end
x（i）=（y（i）-z）/l（i,i）;
end
Example:
>>a=[9-3630;-36192-180;30-180180];
>>b=[111]';
>>x=Chol_Solve（a,b）
x=
1.83331.08330.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）;
fork=1:
n
d（k）=A（k,k）;
forj=1:
k-1
d（k）=d（k）-L（k,j）*T（k,j）;
end
ifabs（d（k））<1e-10
return;
end
fori=k+1:
n
T（i,k）=A（i,k）;
forj=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）;
functionx=LDL_Solve（A,b）
%利用对对称正定矩阵A进行LDL'分解法求解线性方程组Ax=b
n=length（b）;
[l,d]=LDL_Factor（A）;
y
（1）=b
（1）;
fori=2:
n
z=0;
fork=1:
i-1
z=z+l（i,k）*y（k）;
end
y（i）=b（i）-z;
end
x（n）=y（n）/d（n,n）;
fori=n-1:
-1:
1
z=0;
fork=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.83331.08330.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.00920.0095;-0.00920.99030.0136;0.00950.01360.9898];
>>b=[101]';x0=[000]';
>>[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）
ifnargin==3
eps=1.0e-6;
M=200;
elseifnargin<3
error
return
elseifnargin==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;%迭代次数
whilenorm（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）
ifnargin==3
eps=1.0e-6;
M=200;
elseifnargin==4
M=200;
elseifnargin<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;%迭代次数
whilenorm（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）
ifnargin==4
eps=1.0e-6;
M=200;
elseifnargin<4
error
return
elseifnargin==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;%迭代次数
whilenorm（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）
ifnargin==4
eps=1.0e-6;
M=200;
elseifnargin<4
error
return
elseifnargin==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;%迭代次数
whilenorm（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）
ifnargin==3
eps=1.0e-6;
M=200;
elseifnargin<3
error
return
elseifnargin==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;%迭代次数
whilenorm（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-3001050-1400630;
-3004800-1890026880-12600;
1050-1890079380-11760056700;
-140026880-117600179200-88200;
630-1260056700-8820044100;];
b=[53-10-2]';
x0=[00000]';
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．特殊解法
解三对角线性方程组之追赶法：
functionx=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.51.00000;
11.51.0000;
01.00.51.000;
001.00.51.00;
0001.01.51.0;
0000
升级会员 