线性方程组求解Matlab程序.docx

线性方程组求解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-180];

>>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;

00001.02.