MATLAB实现雅可比与高斯塞德尔迭代

概述

用MATLAB编程实现，形成m函数文件。输入A，b矩阵，无返回值，解得x向量直接显示在命令行窗口，同时绘制出x向量的收敛曲线。
$$A=\left[\begin{array}{ccc}2& -1& 1\\ 1& 1& 1\\ 1& 1& -2\end{array}\right]b=\left[\begin{array}{c}1\\ 1\\ 1\end{array}\right]$$

高斯塞德尔方法

代码

function [] = gaussseidel(A,b)
%高斯塞德尔迭代法
G=-(diag(diag(A))+tril(A,-1))\triu(A,1);
R=max(abs(eig(G)));
if(R>=1)
     disp('Gauss-Seidel is not work');
else
     disp('Gauss-Seidel is going well');
     disp('x vector:')
     d1=(diag(diag(A))+tril(A,-1))\b;
     x=zeros(size(b));
     for k=1:100
        x(:,k+1)=G*x(:,k)+d1;
        if(max(abs(x(:,k+1)-x(:,k)))<10^-4)
            break;
        end
     end
     disp(x(:,k+1));
     e=1:k+1;
     plot(e,x,'- o');
     xlabel('generations');
     ylabel('vector_x');
end
end

结果

A=[2 -1 1;1 1 1;1 1 -2];b=[1;1;1];
gaussseidel(A,b);
Gauss-Seidel is going well
x vector:
0.6667
0.3333
0

可见，高斯塞德尔方法收敛且成功求出结果。

雅可比方法

代码

function [] = jacobi(A,b)
%雅可比迭代法
B=-(diag(diag(A)))\(tril(A,-1)+triu(A,1));
R=max(abs(eig(B)));
if(R>=1)
     disp('Jacobi is not work');
else
     disp('Jacobi is going well');
     disp('x vector:')
     d1=diag(diag(A))\b;
     x=zeros(size(b));
     for k=1:100
        x(:,k+1)=B*x(:,k)+d1;
        if(max(abs(x(:,k+1)-x(:,k)))<10^-4)
            break;
        end
     end
     disp(x(:,k+1));
     e=1:k+1;
     plot(e,x,'- o');
     xlabel('generations');
     ylabel('vector_x');
end
end

结果

A=[2 -1 1;1 1 1;1 1 -2];b=[1;1;1];
jacobi(A,b)
Jacobi is not work
可见，雅可比方法不收敛。

代码思路

首先判断收敛性

判断迭代矩阵的谱半径的上界是否小于1。
若大于1则结束程序。

迭代求解

高斯塞德尔迭代矩阵

$$G=-D^{-1}(L+U)$$
$$x^{(k+1)}=G{x}^{(k)}+d_1$$

雅可比迭代矩阵

$$B=-(D+L{)}^{-1}U$$
$$x^{(k+1)}=B{x}^{(k)}+d_1$$
每次迭代后判断
$$|x^{k+1}-x^k|<10^{-4}$$
结果为“真”跳出迭代循环，输出迭代末尾x向量。并绘出收敛曲线。