概述

用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向量。并绘出收敛曲线。

计算线性方程组Ax=b，利用A=D-L-U，将线性方程组改写为x=Bx+f,用迭代的方法进行计算
以6阶希尔伯特矩阵为例，如果需计算其他矩阵，改一下A，要计算其他阶数，改一下n的值就好

clc
clear
format rat%也可以用format long
n=6
%定义希尔伯特矩阵
H=zeros(n,n);
for i=1:n
    for j=1:n
        H(i,j)=1/(i+j-1);
    end
end
H
%此处可定义不同的矩阵A，向量b
A=H;
b=A*ones(n,1);
x0=zeros(n,1);%初始向量
D=zeros(n,n);
L=zeros(n,n);
U=zeros(n,n);
for i=1:n
    D(i,i)=A(i,i);%定义对角阵
end
D
for i=2:n
    for j=1:i-1
        L(i,j)=-A(i,j);%定义单位下三角阵
    end
end
L
for i=1:n
    for j=i+1:n
        U(i,j)=-A(i,j);
    end
end
U
m=1000;%可根据情况改变m的值
t1=0;
t2=0;
%在计算之前先判断两种迭代方法是否收敛
if max(abs(eig(D^-1*(L+U))))<1
    disp('雅可比迭代法收敛，接下来使用此方法进行求解');
    x1=x0;
    for i=1:m
        x1=D^-1*(L+U)*x1+D^-1*b;
        if(abs(t1-x1)<10^-5)
            sprintf('在经过%d次迭代后，得到收敛解',i)
            break
        end
        t1=x1;
    end
    disp('雅可比迭代法得到的x:');
    disp(x1)
else
    disp('雅可比迭代法不收敛，计算停止')
end
if max(abs(eig((D-L)^-1*U)))<1
    disp('高斯-塞德尔迭代法收敛，接下来使用此方法进行求解');
    x2=x0;
    for i=1:m
        x2=(D-L)^-1*U*x2+(D-L)^-1*b;
        if(abs(t2-x2)<10^-5)
            sprintf('在经过%d次迭代后，得到收敛解',i)
            break
        end
        t2=x2;
    end
    disp('高斯-塞德尔迭代法得到的x:');
    disp(x2)
else
    disp('高斯-塞德尔迭代法不收敛，计算停止')
end

matlab数值分析 线性方程组的迭代解法 高斯-赛德尔迭代法

Function [x,iter]=gs(A,b,tol)
D=diag(diag(A));
L=D-tril(A);
U=D-triu(A);
x=zeros(size(b));      %从x=[0;0…]T开始
for iter=1:500
    x=(D-L)\(b+U*x);       %此句换为x=(D)\(b+L*x+U*x);即为Jacobi迭代
    error=norm(b-A*x)/norm(b);
    if(error<tol)
        break;
    end
end

主函数(调用程序)

A=[2,-1,0;-1,3,-1;0,-1,2];b=[1;8;-5];tol=1e-4;[x,iter]=gs(A,b,tol)