用Matlab写的雅各比i和高斯塞德尔以及SOR迭代法

2025-03-26 23:20:28
推荐回答(1个)
回答1:

可以参考MATLAB语言常用算法程序集
SOR迭代法AX=b
function
[x,n]=SOR(A,b,x0,w,eps,M)
if
nargin==4
eps=
1.0e-6;
M
=
200;
elseif
nargin<4
error
return
elseif
nargin
==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;
%迭代次数
while
norm(x-x0)>=eps
x0=x;
x
=B*x0+f;
n=n+1;
if(n>=M)
disp('Warning:
迭代次数太多,可能不收敛!');
return;
end
end
块雅克比迭代法求线性方程组Ax=b的解
function
[x,N]=
BJ(A,b,x0,d,eps,M)
if
nargin==4
eps=
1.0e-6;
M
=
10000;
elseif
nargin<4
error
return
elseif
nargin
==5
M
=
10000;
%参数的默认值
end
NS
=
size(A);
n
=
NS(1,1);
if(sum(d)
~=
n)
disp('分块错误!');
return;
end
bnum
=
length(d);
bs
=
ones(bnum,1);
for
i=1:(bnum-1)
bs(i+1,1)=sum(d(1:i))+1;
%获得对角线上每个分块矩阵元素索引的起始值
end
DB
=
zeros(n,n);
for
i=1:bnum
endb
=
bs(i,1)+d(i,1)-1;
DB(bs(i,1):endb,bs(i,1):endb)=A(bs(i,1):endb,bs(i,1):endb);
%求A的对角分块矩阵
end
for
i=1:bnum
endb
=
bs(i,1)+d(i,1)-1;
invDB(bs(i,1):endb,bs(i,1):endb)=inv(DB(bs(i,1):endb,bs(i,1):endb));
%求A的对角分块矩阵的逆矩阵
end
N
=
0;
tol
=
1;
while
tol>=eps
x
=
invDB*(DB-A)*x0+invDB*b;
%由于LB+DB=DB-A
N
=
N+1;
%迭代步数
tol
=
norm(x-x0);
%前后两步迭代结果的误差
x0
=
x;
if(N>=M)
disp('Warning:
迭代次数太多,可能不收敛!');
return;
end
end
高斯-赛德尔迭代法求线性方程组Ax=b的解
function
[x,N]=
BJ(A,b,x0,d,eps,M)
if
nargin==4
eps=
1.0e-6;
M
=
10000;
elseif
nargin<4
error
return
elseif
nargin
==5
M
=
10000;
%参数的默认值
end
NS
=
size(A);
n
=
NS(1,1);
if(sum(d)
~=
n)
disp('分块错误!');
return;
end
bnum
=
length(d);
bs
=
ones(bnum,1);
for
i=1:(bnum-1)
bs(i+1,1)=sum(d(1:i))+1;
%获得对角线上每个分块矩阵元素索引的起始值
end
DB
=
zeros(n,n);
for
i=1:bnum
endb
=
bs(i,1)+d(i,1)-1;
DB(bs(i,1):endb,bs(i,1):endb)=A(bs(i,1):endb,bs(i,1):endb);
%求A的对角分块矩阵
end
for
i=1:bnum
endb
=
bs(i,1)+d(i,1)-1;
invDB(bs(i,1):endb,bs(i,1):endb)=inv(DB(bs(i,1):endb,bs(i,1):endb));
%求A的对角分块矩阵的逆矩阵
end
N
=
0;
tol
=
1;
while
tol>=eps
x
=
invDB*(DB-A)*x0+invDB*b;
%由于LB+DB=DB-A
N
=
N+1;
%迭代步数
tol
=
norm(x-x0);
%前后两步迭代结果的误差
x0
=
x;
if(N>=M)
disp('Warning:
迭代次数太多,可能不收敛!');
return;
end
end