基于matlab的jacobi(雅可比)迭代法求解线性方程组

1.jacobi迭代法的概念

2.jacobi迭代法的解释

3.jacobi迭代法的推导

4.matlab实现

概念

考虑线性方程组Ax=b

其中A为非奇异矩阵，当A为低阶稠密矩阵是，选主元消去法是有效方法。

但对于A 的阶数n很大，零元素较多的大型稀疏矩阵方程组，利用迭代法求解则更为合适。迭代法通常适用于A中有大量零元素的特点。

解释

简单迭代法是究竟是什么呢？

给下面这个例子你就懂了

例1.已知9x2=sinx+1,在x=0.4附近有根，假定我们已会计算1+sinx−−−−−−−√那么我们就能从x=0.4开始，通过迭代公式xn+1=131+sinx−−−−−−−√逐步求得所要求的根

简单说就是给定一个初始值，代入式子求得一个x，再代入，不断代入，逐步求得根。

下面给出求解方程组的例子。

例2.

求解方程组

⎧⎩⎨⎪⎪⎪⎪8.x1−3x2+2x34x1+11x2−x36x1+3x2+12x3=20=33=36

记为Ax=b,其中

A=⎛⎝⎜846−31132−112⎞⎠⎟,x=⎛⎝⎜x1x2x3⎞⎠⎟,b=⎛⎝⎜203336⎞⎠⎟,

方程组的精确解是

x∗=(3,2,1)T

现在求解这个方程组

方程组可以写为

⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪x1x2x3=18(3x2−2x3+20)=111(−4x1+x3+33)=112(−6x1−3x2+36)1.1式

或写为

x=G0x+g

其中

G0=⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪0−411−612380−312−281110⎫⎭⎬⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪,g=⎛⎝⎜⎜⎜⎜⎜⎜20833113612⎞⎠⎟⎟⎟⎟⎟⎟,

任取初始值，例如取x(0)=(0,0,0)T

将这些值代入1.1式，得到新的值

x(1)=(x(1)1,x(1)2,x(1)3)T=(2.5,3,3)T,

再将x(1)代如1.1式，得到x(2)，

迭代到第十次有

x(10)=(3.000032,1.999838,0.9998813)T

由此可见，由迭代法产生的向量序列x(k)逐步逼近方程组的精确解x∗

推导

由上可知重复步骤得到迭代公式如下：

x(0)=⎛⎝⎜⎜⎜x(0)1x(0)2x(0)3⎞⎠⎟⎟⎟,x(1)=⎛⎝⎜⎜⎜x(1)1x(1)2x(1)3⎞⎠⎟⎟⎟,⋯,x(k)=⎛⎝⎜⎜⎜x(k)1x(k)2x(k)3⎞⎠⎟⎟⎟,⋯,

简写为：

x(k+1)=G0x(k)+g,

雅可比迭代法计算公式简单，每次迭代只需要计算一次矩阵和向量的乘法且计算过程中原始矩阵A始终不变。

迭代法收敛的充要条件是

对任意选取初始向量x^{(0)},G0的谱半径p(G0)<1

所谓“谱半径”，就是最大特征值（对于实数而言），如果最大特征值是复数的话，谱半径就是特征值的最大模。

设aii≠0(i=1,2,...,n)，并将A写成三部分

使得A = D-L-U

D=⎛⎝⎜⎜a110⋱0ann⎞⎠⎟⎟

L=⎛⎝⎜⎜⎜⎜⎜⎜⎜0−a21⋮−an10⋱⋱⋯0ann00⎞⎠⎟⎟⎟⎟⎟⎟⎟

U=⎛⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜0−a120⋯⋱⋱0−a1n⋮−a(n−1)n0⎞⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟

我们可以推导一下G和g,如下：

(D−L−U)x=b

Dx=(L+U)x+b

x=D−1(L+U)x+D−1b

所以

G=D−1(L+U)

g=D−1b

实现

matlab实现代码如下

function [x,n] = jacobi(A,b,x0,eps,varargin)
if nargin ==3
eps = 1.0e-6;
M = 200;
elseif nargin<3
disp('输入参数数目不足3个');
return
elseif nargin ==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;%迭代次数
while norm(x-x0)>=eps
x0 = x;
x = B*x0+f
n = n+1;
if(n>=M)
disp('Warning:迭代次数太多,可能不收敛!')
return;
end
end