MATLAB代码 解线性方程组的迭代法
matlab逐次超松弛迭代法

matlab逐次超松弛迭代法
逐次超松弛迭代法(Gauss-Seidel迭代法)是一种用于解线性方程组的迭代方法,通常用于求解大型稀疏线性方程组。
在MATLAB 中,可以使用该方法来解决线性方程组的数值解。
首先,让我们来了解一下逐次超松弛迭代法的基本原理。
该方法是基于迭代的思想,通过不断迭代更新解向量的各个分量,直到满足一定的收敛条件为止。
具体步骤如下:
1. 首先,需要将线性方程组表示为矩阵形式 Ax = b,其中A 是系数矩阵,x是未知向量,b是常数向量。
2. 然后,将系数矩阵A分解为下三角矩阵L、对角矩阵D和上三角矩阵U,即A = L + D + U。
3. 接下来,可以根据逐次超松弛迭代法的迭代公式来更新解向量x的各个分量,直到满足一定的精度要求或者迭代次数达到指定的值为止。
在MATLAB中,可以通过编写相应的代码来实现逐次超松弛迭代
法。
具体步骤如下:
1. 首先,需要编写一个函数来实现逐次超松弛迭代法的迭代过程,可以使用for循环来进行迭代更新解向量的各个分量。
2. 其次,需要编写主程序来调用该函数,并传入系数矩阵A、常数向量b以及迭代的初始解向量作为输入参数。
3. 最后,可以设置迭代的终止条件,例如迭代次数的最大值或者解的精度要求,以及初始解向量的初值。
需要注意的是,在实际应用中,逐次超松弛迭代法的收敛性和稳定性需要进行分析和验证,以确保得到正确的数值解。
此外,还需要注意选择合适的松弛因子来加速收敛速度。
总的来说,逐次超松弛迭代法是一种常用的求解线性方程组的数值方法,在MATLAB中可以通过编写相应的代码来实现该方法,并得到线性方程组的数值解。
三种迭代法matlab程序 数值分析

• for k=1:max1
• for j=1:N
•
if j==1
•
X(1)=(b(1)-A(1,2:N)*P(2:N))/A(1,1);
•
elseif j==N
•
X(N)=(b(N)-A(N,1:N-1)*(X(1:N-1))')/A(N,N);
•
else
•
X(j)=(b(j)-A(j,1:j-1)*X(1:j-1)-A(j,j+1:N)*P(j+1:N))/A(j,j);
•
end
• end
• err=abs(norm(X'-P));
• P=X';
• if(err<delta)
•
break
• end
• end
• X=X';
• err,k
雅可比迭代法的Matlab程序
给 定 初 始 值 X P0 , 用 雅 克 比 迭 代 法 求 解 线 性 方 程 组
AX b,并生成序列Pk ,求不超过误差界的近似解。
• for k=1:max1
• for j=1:N
•
if j==1
•
X(1)=(b(1)-A(1,2:N)*P(2:N))/A(1,1);
•
elseif j==N
•
X(N)=(b(N)-A(N,1:N-1)*(X(1:N-1))')/A(N,N);
•
else
•
X(j)=(b(j)-A(j,1:j-1)*X(1:j-1)-A(j,j+1:N)*P(j+1:N))/A(j,j);
function X=jacobi(A,b,P,delta,max1) %A是n维非奇异阵。%b是n维向量。%P是初值。%delta是误差界。 %max1是给定的迭代最高次数。%X为所求的方程组AX=b的近似解。 N=length(b); for k=1:max1 for j=1:N
matlab 解线性方程组的迭代法

小结
➢ 线性方程组求根方法的几何意义
➢ 线性方程组求根函数的理解与应用
设线性代数方程组为
展开为
若对角元素 逐一变量分离得方程组
即
此即为迭代公式
简单迭代解法的过程如下:
1 设定一组初值 2 第一次迭代:
得到
第k次迭代 第i个变量
3 第二次迭代: 得到
4 同样做法,得到第k+1次迭代:
迭代次数k的取值与精度要求有关,按下式判断:
若满足则停止迭代 为了便于编程,迭代公式可改写为:
matlab 解线性方程组的 迭代法
2020年4月22日星期三
第十讲 解线性方程组的迭代解法
内容提要
引言 简单迭代法 赛得尔迭代法 迭代解法的收敛性 MATLAB的线性方程组求解函数2 小结
1、引言
迭代解法的基本思想
根据给定方程组,设计出一个迭代公式,构造一 数组的序列 ,代入迭代公式,计算出 ,再代 入迭代公式,经过k次迭代运算后得到 ,若 收敛于某一极限数组xi,则xi就是方程组的近似解。
while(norm(x-x1)>eps) x1=x; x=(I-A)*x1+b; n = n + 1; if(n>=M) disp('Warning: 迭代次数太多,现
在退出!'); return;
end end
例:求解方程组
clear all; A =[ 1.0170 -0.0092 0.0095;
matlabjacobi迭代法

matlabjacobi迭代法Jacobi迭代法是一种求解线性方程组的迭代法,其基本思想是将原方程组的系数矩阵分解为对角部分和非对角部分,对于对角矩阵使用前、后代替法求解,对于非对角部分使用迭代更新法求解。
Jacobi迭代法的基本形式如下:$\begin{cases}a_{11}x_1+a_{12}x_2+...+a_{1n}x_n=b_1 \\a_{21}x_1+a_{22}x_2+...+a_{2n}x_n=b_2 \\... \\a_{n1}x_1+a_{n2}x_2+...+a_{nn}x_n=b_n \\\end{cases}$其中,$a_{ij}$表示系数矩阵的第$i$行第$j$列的元素,$b_i$表示方程组的第$i$个方程的解。
设向量$x^{(k)}=(x_1^{(k)},x_2^{(k)},...,x_n^{(k)})$表示Jacobi迭代法的第$k$次迭代结果,则迭代公式为:$x_i^{(k+1)}=\frac{1}{a_{ii}}(b_i-\sum_{j=1,j\ne i}^n a_{ij}x_j^{(k)}),i=1,2,...,n$迭代公式的意义是,将第$i$个变量的系数$a_{ii}$看成系数矩阵的一个主对角元,将剩下的系数$a_{ij}(i\ne j)$看成非对角元,同时将当前未知量向量$x^{(k)}$看成已知量,利用这些参数求解第$i$个方程中未知量$x_i$。
Jacobi迭代法的收敛条件为原矩阵的对角线元素不为零,且矩阵的任意一行中非对角线元素绝对值之和小于对角线元素绝对值。
在Matlab中,可通过编写函数的方式实现Jacobi迭代法。
函数jacobi实现了迭代公式,并以向量形式返回迭代结果,如下所示:```function xnew = jacobi(A, b, xold)% Jacobi迭代法求解线性方程组Ax=b% A为系数矩阵,b为常数向量,xold为迭代初值% 输出迭代后的解向量xnew% 初始化迭代初值n = length(b);xnew = zeros(n,1);% 迭代更新for i = 1:nxnew(i) = (b(i) - A(i,:)*xold + A(i,i)*xold(i)) / A(i,i);endend```在主程序中可按以下步骤使用函数jacobi求解线性方程组:1.构造系数矩阵A和常数向量b;2.设定迭代初值xold;3.利用jacobi函数求解迭代结果,并对迭代过程进行循环。
matlab 解线性方程组的迭代法

迭代过程本质上就是计算极限的过程,一般不能 得到精确解。
迭代法的优点是程序简单,适合于大型方程组求 解,但缺点是要判断迭代是否收敛和收敛速度问题 。 1. 雅可比(Jacobi(1804-1851))迭代法(简单迭代法) 2. 赛得尔 (Seidel (1821 - 1896))迭代法
2、简单迭代法
while(norm(x-x1)>eps) x1=x; x=(I-A)*x1+b; n = n + 1; if(n>=M) disp('Warning: 迭代次数太多,现
在退出!'); return;
end end
例:求解方程组
clear all; A =[ 1.0170 -0.0092 0.0095;
遗传算法是一种基于自然选择的用于求解有约束和无约束 最优问题的方法。遗传算法反复修改包含若干个体的种群 。遗传算法在每一步中,随机从当前种群中选择若干个个 体作为父辈,并用它们产生下一代子辈。在若干代之后, 种群就朝着最优解“进化”。我们可以利用遗传算法去解决 各种最优化问题,包括目标函数是不连续、不可微、随机 或者高度非线性的问题。
若不满足收敛条件,适当调整方程次序或作一 定的线性组合,就可能满足收敛条件。
5、MATLAB的线性方程组求解函数 2
格式
solve('eqn1','eqn2',...,'eqnN','var1,var2,...,varN')
matlab jacobi迭代法代码

matlab jacobi迭代法代码Matlab是一种常用的数学软件,它具有强大的矩阵计算和绘图功能。
在数值计算中,迭代法是一种重要的求解方法。
本文将介绍如何使用Matlab实现Jacobi迭代法,并运用实例来说明其应用。
Jacobi迭代法是一种经典的迭代法,用于解线性方程组。
它的基本思想是通过迭代逐步逼近方程组的解。
具体而言,对于线性方程组Ax=b,Jacobi迭代法通过以下步骤进行计算:1. 将方程组表示为x=D^(-1)(L+U)x+b的形式,其中D为A的对角矩阵,L为A的严格下三角矩阵,U为A的严格上三角矩阵。
2. 初始化解向量x^(0)为一个初始猜测值,通常取零向量。
3. 根据迭代公式x^(k+1)=D^(-1)(b-(L+U)x^(k)),计算下一迭代解x^(k+1)。
4. 重复步骤3,直到解向量收敛于方程组的解。
下面是一个使用Matlab实现Jacobi迭代法的示例代码:```matlabfunction x = Jacobi(A, b, maxIter, tolerance)n = size(A, 1);x = zeros(n, 1);xPrev = x;iter = 0;while iter < maxIterfor i = 1:nsigma = A(i, 1:i-1) * xPrev(1:i-1) + A(i, i+1:n) * xPrev(i+1:n);x(i) = (b(i) - sigma) / A(i, i);endif norm(x - xPrev) < tolerancebreak;endxPrev = x;iter = iter + 1;endend```在上面的代码中,函数Jacobi接受四个参数:系数矩阵A,右侧常数向量b,最大迭代次数maxIter和收敛容限tolerance。
函数返回解向量x。
在迭代过程中,我们使用了一个for循环来更新解向量x的每个分量。
用matlab解线性方程组

用matlab解线性方程组2008-04-12 17:00一。
高斯消去法1.顺序高斯消去法直接编写命令文件a=[]d=[]'[n,n]=size(a);c=n+1a(:,c)=d;for k=1:n-1a(k+1:n, k:c)=a(k+1:n, k:c)-(a(k+1:n,k)/ a(k,k))*a(k, k:c); %消去endx=[0 0 0 0]' %回带x(n)=a(n,c)/a(n,n);for g=n-1:-1:1x(g)=(a(g,c)-a(g,g+1:n)*x(g+1:n))/a(g,g)end2.列主高斯消去法*由于“[r,m]=max(abs(a(k:n,k)))”返回的行是“k:n,k”内的第几行,所以要通过修正来把m 改成真正的行的值。
该程序只是演示程序,真正机器计算不需要算主元素所在列以下各行应为零的值。
直接编写命令文件a=[]d=[] '[n,n]=size(a);c=n+1a(:,c)=d; %(增广)for k=1:n-1[r,m]=max(abs(a(k:n,k))); %选主m=m+k-1; %(修正操作行的值)if(a(m,k)~=0)if(m~=k)a([k m],:)=a([m k],:); %换行enda(k+1:n, k:c)=a(k+1:n, k:c)-(a(k+1:n,k)/ a(k,k))*a(k, k:c); %消去endendx=[0 0 0 0]' %回带x(n)=a(n,c)/a(n,n);for g=n-1:-1:1x(g)=(a(g,c)-a(g,g+1:n)*x(g+1:n))/a(g,g)end3.分别用顺序高斯消去法和列主高斯消去法解方程组a*x=d,并比较结果a=[0 1 2 3;9 11 23 34;62.5 23.4 15.5 17.2;192.01 124 25.1 59.3] d=[1;1;1;1]顺序高斯消去法:提示“Warning: Divide by zero.” x =NaN NaN NaN NaN 列主高斯消去法:x =-1.2460 2.8796 5.5206 -4.3069由此可见列主高斯消去法可以解决顺序高斯消去法所不能解决的问题。
基于Matlab的解线性方程组的几种迭代法的实现及比较

基于Matlab的解线性方程组的几种迭代法的实现及比较线性方程组的解法有很多种,其中一类常用的方法是迭代法。
迭代法根据一个初值逐步逼近方程组的解,在每一次迭代中利用现有的信息产生新的近似值,并不断地修正。
下面介绍基于Matlab的三种迭代法:雅可比迭代法、高斯-赛德尔迭代法和超松弛迭代法,并进行比较。
1. 雅可比迭代法雅可比迭代法是迭代法中最简单的一种方法。
对于线性方程组Ax=b,雅可比迭代法的迭代公式为:x_{i+1}(j)=1/a_{jj}(b_j-\\sum_{k=1,k\eq j}^n a_{jk}x_i(k))其中,i表示迭代次数,j表示未知数的下标,x_i表示第i次迭代的近似解,a_{jk}表示系数矩阵A的第j行第k列元素,b_j 表示方程组的常数项第j项。
在Matlab中,可以使用以下代码实现雅可比迭代:function [x,flag]=jacobi(A,b,X0,tol,kmax)n=length(b);x=X0;for k=1:kmaxfor i=1:nx(i)=(b(i)-A(i,:)*x+A(i,i)*x(i))/A(i,i);endif norm(A*x-b)<tolflag=1;returnendendflag=0;return其中,参数A为系数矩阵,b为常数项列向量,X0为初值列向量,tol为迭代误差容许值(默认为1e-6),kmax为最大迭代次数(默认为1000)。
函数返回值x为近似解列向量,flag表示是否满足容许误差要求。
2. 高斯-赛德尔迭代法高斯-赛德尔迭代法是雅可比迭代法的改进。
其基本思想是,每次迭代时,利用已经求出的新解中的信息来更新其他未知数的值。
迭代公式为:x_{i+1}(j)=(1/a_{jj})(b_j-\\sum_{k=1}^{j-1}a_{jk}x_{i+1}(k)-\\sum_{k=j+1}^n a_{jk}x_i(k))与雅可比迭代法相比,高斯-赛德尔迭代法的每一次迭代都利用了前面已求得的近似解,因此可以更快地收敛。
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
解线性方程组的迭代法1.rs里查森迭代法求线性方程组Ax=b的解function[x,n]=rs(A,b,x0,eps,M)if(nargin==3)eps=1.0e-6;%eps表示迭代精度M=10000;%M表示迭代步数的限制值elseif(nargin==4)M=10000;endI=eye(size(A));n=0;x=x0;tol=1;%迭代过程while(tol>eps)x=(I-A)*x0+b;n=n+1;%n为最终求出解时的迭代步数tol=norm(x-x0);x0=x;if(n>=M)disp('Warning:迭代次数太多,可能不收敛!');return;endend2.crs里查森参数迭代法求线性方程组Ax=b的解function[x,n]=crs(A,b,x0,w,eps,M)if(nargin==4)eps=1.0e-6;%eps表示迭代精度M=10000;%M表示迭代步数的限制值elseif(nargin==5)M=10000;endI=eye(size(A));n=0;x=x0;tol=1;%迭代过程while(tol>eps)x=(I-w*A)*x0+w*b;n=n+1;%n为最终求出解时的迭代步数tol=norm(x-x0);x0=x;if(n>=M)disp('Warning:迭代次数太多,可能不收敛!');return;endend3.grs里查森迭代法求线性方程组Ax=b的解function[x,n]=grs(A,b,x0,W,eps,M)if(nargin==4)eps=1.0e-6;%eps表示迭代精度M=10000;%M表示迭代步数的限制值elseif(nargin==5)M=10000;endI=eye(size(A));n=0;x=x0;tol=1;%前后两次迭代结果误差%迭代过程while(tol>eps)x=(I-W*A)*x0+W*b;%迭代公式n=n+1;%n为最终求出解时的迭代步数tol=norm(x-x0);x0=x;if(n>=M)disp('Warning:迭代次数太多,可能不收敛!');return;endend4.jacobi雅可比迭代法求线性方程组Ax=b的解function[x,n]=jacobi(A,b,x0,eps,varargin)if nargin==3eps=1.0e-6;M=200;elseif nargin<3errorreturnelseif nargin==5M=varargin{1};endD=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;%迭代次数while norm(x-x0)>=epsx0=x;x=B*x0+f;n=n+1;if(n>=M)disp('Warning:迭代次数太多,可能不收敛!');return;endend5.gauseidel高斯-赛德尔迭代法求线性方程组Ax=b的解function[x,n]=gauseidel(A,b,x0,eps,M)if nargin==3eps=1.0e-6;M=200;elseif nargin==4M=200;elseif nargin<3errorreturn;endD=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;%迭代次数while norm(x-x0)>=epsx0=x;x=G*x0+f;n=n+1;if(n>=M)disp('Warning:迭代次数太多,可能不收敛!');return;endend6.SOR超松弛迭代法求线性方程组Ax=b的解function[x,n]=SOR(A,b,x0,w,eps,M)if nargin==4eps=1.0e-6;M=200;elseif nargin<4errorreturnelseif nargin==5M=200;endif(w<=0||w>=2)error;return;endD=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)>=epsx0=x;x=B*x0+f;n=n+1;if(n>=M)disp('Warning:迭代次数太多,可能不收敛!');return;endend7.SSOR对称逐次超松弛迭代法求线性方程组Ax=b的解function[x,n]=SSOR(A,b,x0,w,eps,M)if nargin==4eps=1.0e-6;M=200;elseif nargin<4errorreturnelseif nargin==5M=200;endif(w<=0||w>=2)error;return;endD=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;%迭代次数while norm(x-x0)>=epsx0=x;x12=B1*x0+f1;x=B2*x12+f2;n=n+1;if(n>=M)disp('Warning:迭代次数太多,可能不收敛!');return;endend8.JOR雅可比超松弛迭代法求线性方程组Ax=b的解function[x,n]=JOR(A,b,x0,w,eps,M)if nargin==4eps=1.0e-6;M=10000;elseif nargin==5M=10000;endif(w<=0||w>=2)%收敛条件要求error;return;endD=diag(diag(A));%求A的对角矩阵B=w*inv(D);%迭代过程x=x0;n=0;%迭代次数tol=1;%迭代过程while tol>=epsx=x0-B*(A*x0-b);n=n+1;tol=norm(x-x0);x0=x;if(n>=M)disp('Warning:迭代次数太多,可能不收敛!');return;endend9.twostep两步迭代法求线性方程组Ax=b的解function[x,n]=twostep(A,b,x0,eps,varargin)if nargin==3eps=1.0e-6;M=200;elseif nargin<3errorreturnelseif nargin==5M=varargin{1};endD=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;%迭代次数while norm(x-x0)>=epsx0=x;x12=B1*x0+f1;x=B2*x12+f2;n=n+1;if(n>=M)disp('Warning:迭代次数太多,可能不收敛!');return;endend10.fastdown最速下降法求线性方程组Ax=b的解function[x,n]=fastdown(A,b,x0,eps)if(nargin==3)eps=1.0e-6;endx=x0;n=0;tol=1;while(tol>eps)%以下过程可参考算法流程r=b-A*x0;d=dot(r,r)/dot(A*r,r);x=x0+d*r;tol=norm(x-x0);x0=x;n=n+1;end11.conjgrad共轭梯度法求线性方程组Ax=b的解function[x,n]=conjgrad(A,b,x0)r1=b-A*x0;p=r1;n=0;for i=1:rank(A)%以下过程可参考算法流程if(dot(p,A*p)< 1.0e-50)%循环结束条件break;endalpha=dot(r1,r1)/dot(p,A*p);x=x0+alpha*p;r2=r1-alpha*A*p;if(r2< 1.0e-50)%循环结束条件break;endbelta=dot(r2,r2)/dot(r1,r1);p=r2+belta*p;n=n+1;end12.preconjgrad预处理共轭梯度法求线性方程组Ax=b的解function[x,n]=preconjgrad(A,b,x0,M,eps)if nargin==4eps=1.0e-6;endr1=b-A*x0;iM=inv(M);z1=iM*r1;p=z1;n=0;tol=1;while tol>=epsalpha=dot(r1,z1)/dot(p,A*p);x=x0+alpha*p;r2=r1-alpha*A*p;z2=iM*r2;belta=dot(r2,z2)/dot(r1,z1);p=z2+belta*p;n=n+1;tol=norm(x-x0);x0=x;%更新迭代值r1=r2;z1=z2;end13.BJ块雅克比迭代法求线性方程组Ax=b的解function[x,N]=BJ(A,b,x0,d,eps,M)if nargin==4eps=1.0e-6;M=10000;elseif nargin<4errorreturnelseif nargin==5M=10000;%参数的默认值endNS=size(A);n=NS(1,1);if(sum(d)~=n)disp('分块错误!');return;endbnum=length(d);bs=ones(bnum,1);for i=1:(bnum-1)bs(i+1,1)=sum(d(1:i))+1;%获得对角线上每个分块矩阵元素索引的起始值endDB=zeros(n,n);for i=1:bnumendb=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的对角分块矩阵endfor i=1:bnumendb=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的对角分块矩阵的逆矩阵endN=0;tol=1;while tol>=epsx=invDB*(DB-A)*x0+invDB*b;%由于LB+DB=DB-AN=N+1;%迭代步数tol=norm(x-x0);%前后两步迭代结果的误差x0=x;if(N>=M)disp('Warning:迭代次数太多,可能不收敛!');return;endend14.BGS块高斯-赛德尔迭代法求线性方程组Ax=b的解function[x,N]=BGS(A,b,x0,d,eps,M)if nargin==4eps=1.0e-6;M=10000;elseif nargin<4errorreturnelseif nargin==5M=10000;endNS=size(A);n=NS(1,1);bnum=length(d);bs=ones(bnum,1);for i=1:(bnum-1)bs(i+1,1)=sum(d(1:i))+1;%获得对角线上每个分块矩阵元素索引的起始值endDB=zeros(n,n);for i=1:bnumendb=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的对角分块矩阵endLB=-tril(A-DB);%求A的下三角分块阵UB=-triu(A-DB);%求A的上三角分块阵N=0;tol=1;while tol>=epsinvDL=inv(DB-LB);x=invDL*UB*x0+invDL*b;%块迭代公式N=N+1;tol=norm(x-x0);x0=x;if(N>=M)disp('Warning:迭代次数太多,可能不收敛!');return;endend15.BSOR块逐次超松弛迭代法求线性方程组Ax=b的解function[x,N]=BSOR(A,b,x0,d,w,eps,M)if nargin==5eps=1.0e-6;M=10000;elseif nargin<5errorreturnelseif nargin==6M=10000;%参数默认值endNS=size(A);n=NS(1,1);bnum=length(d);bs=ones(bnum,1);for i=1:(bnum-1)bs(i+1,1)=sum(d(1:i))+1;%获得对角线上每个分块矩阵元素索引的起始值endDB=zeros(n,n);for i=1:bnumendb=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的对角矩阵endLB=-tril(A-DB);%求A的下三角阵UB=-triu(A-DB);%求A的上三角阵N=0;tol=1;iw=1-w;while tol>=epsinvDL=inv(DB-w*LB);x=invDL*(iw*DB+w*UB)*x0+w*invDL*b;%块迭代公式N=N+1;tol=norm(x-x0);x0=x;if(N>=M)disp('Warning:迭代次数太多,可能不收敛!');return;endend。