用Matlab写的雅各比i和高斯塞德尔以及SOR迭代法
可以参考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
function [v,sN,vChain]=gaussSeidel(A,b,x0,errorBound,maxSp)
%Gauss-Seidel迭代法求解线性方程组
%A-系数矩阵 b-右端向量 x0-初始迭代点 errorBound-近似精度 maxSp-最大迭代次数
%v-近似解 sN-迭代次数 vChain-迭代过程的所有值
step=0;
error=inf;
s=size(A);
D=zeros(s(1));
vChain=zeros(15,3);%最多能记录15次迭代次数
k=1;
fx0=x0;
for i=1:s(1)
D(i,i)=A(i,i);
end;
L=-tril(A,-1);
U=-triu(A,1);
while error>=errorBound & step<maxSp
x0=inv(D)*(L+U)*x0+inv(D)*b;
vChain(k,:)=x0';
k=k+1;
error=norm(x0-fx0);
fx0=x0;
step=step+1;
end
v=x0;
sN=step;
(1) 编程求解，并与用数学软件求解的结果对比。
(2) 考察迭代法的收敛性，若均收敛，对比两种方法的收敛速度。
解：源程序：
①雅克比迭代法：建立函数文件jacobi.m
function [n,x]=jacobi(A,b,X,nm,w)
%用雅克比迭代法求解方程组Ax=b
%输入：A为方程组的系数矩阵，b为方程组右端的列向量，X为迭代初值构成的列向量，nm为最大迭代次数，w为误差精度
%输出：x为求得的方程组的解构成的列向量，n为迭代次数
n=1;
m=length(A);
D=diag(diag(A)); %令A=D-L-U,计算矩阵D
L=tril(-A)+D; %令A=D-L-U,计算矩阵L
U=triu(-A)+D; %令A=D-L-U,计算矩阵U
M=inv(D)*(L+U); %计算迭代矩阵
g=inv(D)*b; %计算迭代格式中的常数项
%下面是迭代过程
while n<=nm
x=M*X+g; %用迭代格式进行迭代
if norm(x-X,2)<w
disp('迭代次数为');n
disp('方程组的解为');x
return;
%上面：达到精度要求就结束程序，输出迭代次数和方程组的解
end
X=x;n=n+1;
end
%下面：如果达到最大迭代次数仍不收敛，输出警告语句及迭代的最终结果（并不是方程组的解）
disp('在最大迭代次数内不收敛!');
disp('最大迭代次数后的结果为');x
②高斯赛德尔迭代法：建立函数文件gaussseidel.m
function [n,x]=gaussseidel(A,b,X,nm,w)
%用高斯-赛德尔迭代法求解方程组Ax=b
%输入：A为方程组的系数矩阵，b为方程组右端的列向量，X为迭代初值构成的列向量，nm为最大迭代次数，w为误差精度
%输出：x为求得的方程组的解构成的列向量，n为迭代次数
n=1;
m=length(A);
I=eye(m); %生成m*m阶的单位矩阵
D=diag(diag(A)); %令A=D-L-U,计算矩阵D
L=tril(-A)+D; %令A=D-L-U,计算矩阵L
U=triu(-A)+D; %令A=D-L-U,计算矩阵U
M=inv(D-L)*U; %计算迭代矩阵
g=inv(I-inv(D)*L)*(inv(D)*b); %计算迭代格式中的常数项
%下面是迭代过程
while n<=nm
x=M*X+g; %用迭代格式进行迭代
if norm(x-X,2)<w
disp('迭代次数为');n
disp('方程组的解为');x
return;
%上面：达到精度要求就结束程序，输出迭代次数和方程组的解
end
X=x;n=n+1;
end
%下面：如果达到最大迭代次数仍不收敛，输出警告语句及迭代的最终结果（并不是方程组的解）
disp('在最大迭代次数内不收敛!');
disp('最大迭代次数后的结果为');x
运行过程及结果：
①雅克比迭代法的运行过程及结果为：
>> A=[10,3,1;2,-11,3;1,3,12];
b=[2;-5;4];
X=[0;0;0];
nm=50;
w=10^-6;
>> jacobi(A,b,X,nm,w)
迭代次数为
n =
14
方程组的解为
x =
0.0254
0.5144
0.2026
②高斯赛德尔迭代法的运行过程及结果为：
>> A=[10,3,1;2,-11,3;1,3,12];
b=[2;-5;4];
X=[0;0;0];
nm=50;
w=10^-6;
>> gaussseidel(A,b,X,nm,w)
迭代次数为
n =
6
方程组的解为
x =
0.0254
0.5144
0.2026
③用matlab中的A\b命令的运行过程及结果如下：
>> A=[10,3,1;2,-11,3;1,3,12];
b=[2;-5;4];
>> A\b
ans =
0.0254
0.5144
0.2026
(1)由运行结果可知，编程得到的方程组的解为[0.0254;0.5144;0.2026]。而用matlab中的A\b命令求出的方程组的解为[0.0254;0.5144;0.2026]，完全一致。
(2)由运行过程可知，两种方法均收敛，雅克比迭代次数为14，高斯赛德尔迭代次数为6，说明后者的迭代速度比前者快。
2.用超松弛迭代法求解方程组，分别取松弛因子 ，取迭代初值[0;0;0]，迭代30次或满足 时停止计算。
(1) 编程求解，并与用数学软件求解的结果对比。
(2) 考察迭代是否收敛，若收敛，松弛因子取何值时收敛最快，与有关定理的结论对照，看结果是否一致。
解：源程序：
①超松弛迭代法：建立函数文件sor22.m
function [n,x]=sor22(A,b,X,nm,w,ww)
%用超松弛迭代法求解方程组Ax=b
%输入：A为方程组的系数矩阵，b为方程组右端的列向量，X为迭代初值构成的列向量，nm为最大迭代次数，w为误差精度，ww为松弛因子
%输出：x为求得的方程组的解构成的列向量，n为迭代次数
n=1;
m=length(A);
D=diag(diag(A)); %令A=D-L-U,计算矩阵D
L=tril(-A)+D; %令A=D-L-U,计算矩阵L
U=triu(-A)+D; %令A=D-L-U,计算矩阵U
M=inv(D-ww*L)*((1-ww)*D+ww*U); %计算迭代矩阵
g=ww*inv(D-ww*L)*b; %计算迭代格式中的常数项
%下面是迭代过程
while n<=nm
x=M*X+g; %用迭代格式进行迭代
if norm(x-X,'inf')<w
disp('迭代次数为');n
disp('方程组的解为');x
return;
%上面：达到精度要求就结束程序，输出迭代次数和方程组的解
end
X=x;n=n+1;
end
%下面：如果达到最大迭代次数仍不收敛，输出警告语句及迭代的最终结果（并不是方程组的解）
disp('在最大迭代次数内不收敛!');
disp('最大迭代次数后的结果为');x
②向后超松弛迭代法：建立函数文件sorxh22.m
function [n,x]=sorxh22(A,b,X,nm,w,ww)
%用向后超松弛迭代法求解方程组Ax=b
%输入：A为方程组的系数矩阵，b为方程组右端的列向量，X为迭代初值构成的列向量，nm为最大迭代次数，w为误差精度，ww为松弛因子
%输出：x为求得的方程组的解构成的列向量，n为迭代次数
n=1;
m=length(A);
D=diag(diag(A)); %令A=D-L-U,计算矩阵D
L=tril(-A)+D; %令A=D-L-U,计算矩阵L
U=triu(-A)+D; %令A=D-L-U,计算矩阵U
M=inv(D-ww*U)*((1-ww)*D+ww*L); %计算矩阵迭代矩阵
g=ww*inv(D-ww*U)*b; %计算迭代格式中的常数项
%下面是迭代过程
while n<=nm
x=M*X+g; %用迭代格式进行迭代
if norm(x-X,'inf')<w
disp('迭代次数为');n
disp('方程组的解为');x
return;
%上面：达到精度要求就结束程序，输出迭代次数和方程组的解
end
X=x;n=n+1;
end
%下面：如果达到最大迭代次数仍不收敛，输出警告语句及迭代的最终结果（并不是方程组的解）
disp('在最大迭代次数内不收敛!');
disp('最大迭代次数后的结果为');x
③对称超松弛迭代法：建立函数文件ssor22.m
function [n,x]=ssor22(A,b,X,nm,w,ww)
%用对称超松弛迭代法求解方程组Ax=b
%输入：A为方程组的系数矩阵，b为方程组右端的列向量，X为迭代初值构成的列向量，nm为最大迭代次数，w为误差精度，ww为松弛因子
%输出：x为求得的方程组的解构成的列向量，n为迭代次数
n=1;
m=length(A);
I=eye(m); %生成m*m阶的单位矩阵
D=diag(diag(A)); %令A=D-L-U,计算矩阵D
L=tril(-A)+D; %令A=D-L-U,计算矩阵L
U=triu(-A)+D; %令A=D-L-U,计算矩阵U
%下面计算迭代矩阵和迭代格式中的常数项
M=inv(D-ww*U)*((1-ww)*D+ww*L)*inv(D-ww*L)*((1-ww)*D+ww*U); g=ww*inv(D-ww*U)*(I+((1-ww)*D+ww*L)*inv(D-ww*L))*b;
%下面是迭代过程
while n<=nm
x=M*X+g; %用迭代格式进行迭代
if norm(x-X,'inf')<w
disp('迭代次数为');n
disp('方程组的解为');x
return;
%上面：达到精度要求就结束程序，输出迭代次数和方程组的解
end
X=x;n=n+1;
end
%下面：如果达到最大迭代次数仍不收敛，输出警告语句及迭代的最终结果（并不是方程组的解）
disp('在最大迭代次数内不收敛!');
disp('最大迭代次数后的结果为');x
可以参考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
诶呀,我这学期正好在学~可是都没去上课~
...远场噪声呢?能否用单极子,偶极子举例,用matlab能实现吗
工程实际中作用又相反，因此在近似的处理上可相消，近似地认识其对系统的作用相互抵消了。有了主导极点和偶极子对的概念后，对于高阶系统的分析，在误差精度允许的情况下，可将高阶系统的主导极点分析出来，利用主导极点来分析系统，相当于降低了系统的阶数，给分析带来方便。雅岚？？？
我们老师预习五年级《难忘的一课》这篇课文,要求查找台湾的资料你可以给...
雅美族的甩发舞与雕船本领,阿美族的对位唱法与多姿的舞蹈,排湾族的刺绣与雕刻,布农族的‘打耳祭‘与‘成年祭‘大典,赛夏族的‘矮灵祭‘与佩铃叮咚的舞蹈,泰雅族的绣衣与播种节,卑南族的‘刹猴祭‘与‘狩猎祭‘等,是那样五彩缤纷。为发展旅游业,采取了一些措施,如在桃园县兴建‘小人国‘,更丰富了旅游资源。
带有歌词 嘿我真的好想你,是什么歌曲?求助
苦等2小时: 最火的社交软件是什么 回答 苦等3小时: 你以前喜欢的女生不知道你已经结婚离婚了。聊了几次也... 回答 苦等5小时: 谁知道闲鱼,里面一个叫永不言弃的用户为什么违规 20 回答 苦等5小时: 如何在matlab中对一组变量使用滑动窗口 10 回答 更多等待求助问题 > 登录...
滑模控制的基本原理
滑模变结构控制的原理，是根据系统所期望的动态特性来设计系统的切换超平面，通过滑动模态控制器使系统状态从超平面之外向切换超平面收束。系统一旦到达切换超平面，控制作用将保证系统沿切换超平面到达系统原点，这一沿切换超平面向原点滑动的过程称为滑模控制。由于系统的特性和参数只取决于设计的切换超平面而与...
偶蒲盘得： 1. 用雅克比迭代法和高斯--赛德尔迭代法求解下列方程组,取迭代初值[0;0;0]. (1) 编程求解,并与用数学软件求解的结果对比. (2) 考察迭代法的收敛性,若均收敛,对比两种方法的收敛速度.解:源程序: ①雅克比迭代法:建立函数文件...
怀集县17217876582： Matlab实现雅各比矩阵 - ？
偶蒲盘得： function X=jacobi(A,b,P,delta,max1) % A是n维非奇异阵 % B是n维向量 % P是初值 % delta是误差界 % X为所求的方程组AX=B的近似解 N=length(b); for k=1:max1 for j=1:N X(j)=(b(j)-A(j,[1:j-1,...
怀集县17217876582： 用matlab程序,编写出高斯塞德尔迭代法 - ？
偶蒲盘得： function [v,sN,vChain]=gaussSeidel(A,b,x0,errorBound,maxSp) %Gauss-Seidel迭代法求解线性方程组 %A-系数矩阵 b-右端向量 x0-初始迭代点 errorBound-近似精度 maxSp-最大迭代次数 %v-近似解 sN-迭代次数 vChain-迭代过程的所有值 step=0;...
怀集县17217876582： 求用matlab程序,编写出高斯塞德尔迭代法 - ？
偶蒲盘得： function [x,k,index]=GS(A,b,ep,it_max) n=length(A);k=0; x=zeros(n,1); y=zeros(n,1); index=1; for k=1:it_max y=x; for i=1:n z=b(i); for j=1:n; if j~=i z=z-A(i,j)*x(j); end end if abs(A(i,i))<1e-10 index=0; return; end z=z/A(i:i); x(i)=z; end if norm(y-x,inf)<ep break; endend
怀集县17217876582： 雅可比迭代法与高斯塞德尔迭代法的区别与特征 - ？
偶蒲盘得： Jacobi与Gauss-Seidel迭代法Jacobi(雅可比)迭代法 我们从形式入手学习J和GS迭代法.先使用雅可比方法: 例:解方程组 解: 初值 迭代两步则有下表:k 0 0.1250 0.4000 -0.6000 1 0.2500 0.3150 -0.4950 2 0.2263 0.3005 -0.4870 由解的...
怀集县17217876582： 求!!在MATLAB中用求解高斯—塞德尔迭代法线性方程组程序 - ？
偶蒲盘得： 在Matlab命令行下输入:edit Gauss_seidel.m 然后将下面输入,并保存 function x=Gauss_Seidel(A,b,x0,tol) if (nargin==2) x0=ones(size(b)); tol=1e-6; elseif (nargin==3) tol=1e-6; else sprintf('USAGE:Gauss_Seidel(A,b,x0,tol)') end D=diag(diag(A));...
怀集县17217876582： matlab 计算节点温度高斯塞尔曼迭代算法 - ？
偶蒲盘得： 求解思路:1、根据题意,列出25组方程系数矩阵 A;2、根据边界条件,列出边界向量值 B;3、利用高斯——塞尔曼迭代算法,求解 X值;4、生成网格数据点矩阵,绘制温度分布图.
怀集县17217876582： 有谁知道雅可比、高斯 -- 塞德尔迭代法的程序实现?? - ？
偶蒲盘得： ss_seidel.m 然后将下面输入,并保存 function x=Gauss_Seidel(A,b,x0,tol) if (nargin==2) x0=ones(size(b)); tol=1e-6; elseif (nargin==3) tol=1e-6; else sprintf('USAGE:Gauss_Seidel(A,b,x0,tol)') endD=diag(diag(A)); U=triu(A,1); L=tril(A,-1); G=-(D+L...
怀集县17217876582： 如何判断雅各比迭代法、高斯赛德尔迭代法是否收敛 - ？
偶蒲盘得： 计算谱半径,谱半径小于1,则收敛,否则不收敛.其中谱半径就是迭代矩阵J或者G的最大特征值!!望采纳!!不懂再问!也可用列范数或行范数判断,列范数或者行范数小于1,则收敛.但范数大于1时,不能说明其发散,还要通过计算谱半径来确定其收敛性.
怀集县17217876582： 矩阵A要满足什么条件时,高斯赛德尔法比雅格比法收敛的快. - ？
偶蒲盘得： // seidel.h: interface for the cseidel class. // ////////////////////////////////////////////////////////////////////// #if !defined(afx_seidel_h__35754d65_c3b8_4814_b9d7_8de3ba72eff3__included_) #define afx_seidel_h__35754d65_c3b8_4814_b9d7_8de3ba72eff3__...
