用matlab解线性方程组.docx
《用matlab解线性方程组.docx》由会员分享，可在线阅读，更多相关《用matlab解线性方程组.docx（35页珍藏版）》请在冰豆网上搜索。
用matlab解线性方程组
用matlab解线性方程组
2008-04-1217:
00
一。
高斯消去法
1.顺序高斯消去法
直接编写命令文件
a=[]
d=[]'
[n,n]=size（a）;
c=n+1
a（:
c）=d;
fork=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=[0000]'%回带
x（n）=a（n,c）/a（n,n）;
forg=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;%（增广）
fork=1:
n-1
[r,m]=max（abs（a（k:
n,k）））;%选主
m=m+k-1;%（修正操作行的值）
if（a（m,k）~=0）
if（m~=k）
a（[km],:
）=a（[mk],:
）;%换行
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=[0000]'%回带
x（n）=a（n,c）/a（n,n）;
forg=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=[0123;9112334;62.523.415.517.2;192.0112425.159.3]d=[1;1;1;1]
顺序高斯消去法：
提示“Warning:
Dividebyzero.”x=NaNNaNNaNNaN
列主高斯消去法：
x=-1.24602.87965.5206-4.3069
由此可见列主高斯消去法可以解决顺序高斯消去法所不能解决的问题。
4.将上述矩阵中的“2”改为2.05，“34”改为“34.6”，“15.5”改为“15.57”，“124”改为“124.7”再用列主高斯消去法计算，与上述结果比较。
x=-0.80811.82263.5568-2.7047
很明显虽然系数矩阵只有很小的变化但结果的变化却很大，所以系数矩阵是病态的。
这是因为系数矩阵的条件数很大：
cond（a）=2.0695e+003
二。
迭代法
J迭代公式
a=[521;-142;2-310]
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;
whilenorm（X-x,inf）>=stop%时迭代中止否则继续
x=X;
X=D*（L+U）*x+D*d;
n=n+1;
end
X
n
G-S迭代公式
a=[521;-142;2-310]
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;
whilenorm（X-x,inf）>=stop%时迭代中止否则继续
x=X;
X=inv（D-L）*U*x+inv（D-L）*d;
n=n+1;
end
X
n
SOR迭代公式
a=[521;-142;2-310]
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;
whilenorm（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=da=[521;-142;2-310]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.00003.00002.0000n=18
G-S：
X=-4.00003.00002.0000n=8
SOR（）：
X=-4.00003.00002.0000n=482
可见高斯-赛德尔迭代法要比雅可比迭代公式的收敛速度快。
同时，如果超松弛迭代法的选取不当不但不会加速收敛还会对收敛速度造成很恶劣的结果。
藏
线性方程组求解
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;%迭代次数
升级会员 