线性方程组求解Matlab程序.docx
- 文档编号:2256455
- 上传时间:2022-10-28
- 格式:DOCX
- 页数:31
- 大小:21.61KB
线性方程组求解Matlab程序.docx
《线性方程组求解Matlab程序.docx》由会员分享,可在线阅读,更多相关《线性方程组求解Matlab程序.docx(31页珍藏版)》请在冰豆网上搜索。
线性方程组求解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=dia
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 线性方程组 求解 Matlab 程序