文章目录
 
 
前言
一、Jacobi迭代法是什么？
二、对应的编程思想以及公式推导 
   
1.Jacobi迭代法 公式推导
2.Jacobi迭代法求解线性方程组 例子
3.Jacobi迭代法 编程实现
总结
 
 


 前言
 
 
 
                                       雅克比（Jacobi）迭代法求解线性方程组
 
 

一、Jacobi迭代法是什么？ 
 
 
 
 
简单的讲其实就是我们平时求解的方法（最常用的方法）
 
 
以下是Jacobi的迭代过程：
 
 
 
 
 
 
 
 
 
 
 

二、对应的编程思想以及公式推导
 
 

1.Jacobi迭代法 公式推导
 
 
 
线性方程组为：                         
 
 
将A分裂：                                   
 
 
从而得到迭代公式为：     
 
 
由于D为对角元素，从而有：          
 
 
即编程语言可以写为：       
 
 
推到过程为：
 
 
                                            
 
 

2.Jacobi迭代法求解线性方程组 例子
 
 
 
 
 
 
 
 
 
 
 
 
3.Jacobi迭代法 编程实现
 
 
 
完整程序代码为：
 
 
function [x,error,iter]=GJJacobi_solve(A,b,epsilon,max_iter)
% x , error , iter 为输出变量     A ，b epsilon,max_iter为输入变量
n=length(b);
x=zeros(n,1);
% 首先编写第一步代码目的  x =-(L+U)*x/D+b/D
%第二部明确符号的意义 D表示对角线的元素 L表示下三角的元素 U表示上三角的元素 b表示列向量的长度
%确定误差error = norm(x-y)

%输入矩阵A 以及 列向量b


%表示L  U  D
L=tril(A); %tril（A）表示取矩阵A的下三角元素,  包含了对角线元素
U=triu(A); %triu(A) 表示取矩阵A的上三角元素，  包含了对角线元素
D=diag(diag(A)); %表示取矩阵A的主对角元素



%想办法使得L，U对角线元素都为0
                             % L(logical(eye(size(L))))=0; %由于维度不对，从而只是等号右边为零不对
                             %  U(logical(eye(size(U))))=0;
%改正
                              %B =  zeros(1,n);              %[0,0,0]
L(logical(eye(size(L))))=0;   %[8;11;12]
U(logical(eye(size(U))))=0;

                             

%需要看看是否输出的是对角线元素为0
disp(L)
disp(U)







% 先写一个循环(发现运用一个变量时需要先初始化error=1)
% 设置最小误差为10e-6
error = 1 ;%初始化误差变量
iter = 0;  %初始化迭代步数变量
while error>epsilon && iter<max_iter % 给出循环条件，以及函数输入变量epsilon,max_iter
    for i =1:max_iter
        y=x;%迭代前的矩阵，用来计算迭代误差
        x=-D\(L+U)*x+D\b ;  %从这里知道需要把L,U,D,b表示出来 x=-(L+U)*x/D+b/D 这个维度错误
        error=norm(x-y);  %知道迭代前的数据y，以及迭代后的数据x，设定误差为（x-y)的范数
        
        if error>epsilon %给出判断误差是否减小
            iter =iter+1;
            break
        end     %if 需要 end 结束
    end         %for 需要 end 结束
end             %while 需要 end 结束
 
 
测试代码为：
 
 
 A=[8 -3 2
    4 11 -1
    6 3 12];
b=[20;33;36];
epsilon=10e-6;
max_iter=15;
[x,error,iter]=GJJacobi_solve(A,b,epsilon,max_iter);
disp('程序计算的精确解为：');
disp(x);
disp('最大迭代次数下的误差：');
disp(error);
disp('最小迭代次数：');
disp(iter+1);
 
 
运行结果为：
 
 
     0     0     0
     4     0     0
     6     3     0

     0    -3     2
     0     0    -1
     0     0     0

程序计算的精确解为：
    3.0000
    2.0000
    1.0000

最大迭代次数下的误差：
   4.1060e-12

最小迭代次数：
    14
 
 
 
 


 总结
 
以上就是今天要讲的内容，本文仅仅简单介绍了雅克比迭代求解线性方程组

function x=myjacob(A,b,x0)
myeps=0.0001;
n=length(b);
err0=1;
while err0>myeps
    x1=zeros(n,1);
    for i=1:n
        s1=A(i,i:i-1)*x0(1:i-1);
        s2=A(i,i+1:n)*x0(i+1:n);
        x1(i)=(b(i)-s1-s2)/A(i,i);
    end
    err0=norm(x1-x0);x0=x1;
end
x=x1;
end
 

计算线性方程组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

说明推导见此博客：https://blog.csdn.net/zengxyuyu/article/details/53054880
 
源码见下面：
 
main.m
 
clear
clc
A = [8 -3 2;4 11 -1;6 3 12];
b = [20;33;36];
[x, n] = jacobi(A,b,[0,0,0]',1.0e-7,30)
 
jacobi.m 
 
function [x,n] = jacobi(A,b,x0,eps,it_max)
%  求线性方程组的Jacobi迭代法，调用格式为

%  [x, k] = jacobi(A,b,x0,eps,it_max)

%  其中, A 为线性方程组的系数矩阵，b 为常数项，eps 为精度要求，默认为1e-6,

%  it_max 为最大迭代次数，默认为200

%  x 为线性方程组的解，k迭代次数
if nargin ==3
    eps = 1.0e-6;
    M = 200;
elseif nargin<3
    disp('输入参数数目不足3个');
    return
elseif nargin ==5
    M = it_max;
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