标准雅可比迭代算法-实现

发表于 分类于 阅读次数： Changyan：

标准雅可比迭代算法解线性方程组-算法实现

  • 函数实现
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
function [ output_args ] = cqj_Jacobi(init_X,coefficient_mat,b_mat,maxNum,epsilon)
% 雅可比迭代法
% init_X：初始值
% coefficient_mat：系数矩阵
% b_mat：b
% maxNum：最大迭代次数
% epsilon：精度
output_args = [];
flag = 0;
%% 预处理 检查用户调用函数时参数格式是否传错
if size(coefficient_mat,1) ~= size(coefficient_mat,2)
    disp('系数矩阵格式错误(需要方阵)');
    return;
end
if size(b_mat,1) ~= 1 && size(b_mat,2) ~= 1
    disp('AX=b中的b必须为行向量或者列向量');
    return;
end
if size(b_mat,1) == 1
    b_mat = b_mat';
end
if size(init_X,1) ~= 1 && size(init_X,2) ~= 1
    disp('初始值必须为行向量或者列向量');
    return;
end
if size(init_X,1) == 1
    init_X = init_X';
end
if size(init_X,1) ~= size(coefficient_mat,1)
    disp('初始值与系数矩阵维数不一致');
    return;
end
if size(init_X,1) ~= size(b_mat,1)
    disp('AX=b中的b向量维数错误(需与系数矩阵的维数一致');
    return;
end
if ~exist('maxNum', 'var')
    maxNum = 500; % 默认值
end
if ~exist('epsilon', 'var')
    epsilon = 1e-14; % 默认值
end
%% 对角线
D = diag(diag(coefficient_mat));
% 可能matlab的源码inv函数对角函数求逆会像普通函数一样处理，
% 这种特殊函数求逆会比较快
D_inv = diag(1./diag(coefficient_mat));
%% 开始迭代
% 下面代码段也行 但是对内存的消耗比较大 
% x(:,0) = init_X;
% for k = 1:maxNum
%     x(:,k) = D_inv*(b_mat-(coefficient_mat-D)*x(:,k-1));
%     if norm(x(:,k)-x(:,k-1)) < epsilon  % 达到精度
%         flag = 1;
%         break;
%     end
% end
x = init_X;
%%
% 迭代公式：
% $$x_{k+1}=D^{-1}(b-(L+U)x_k),k=0,1,2,\cdots$
% 
for k = 1:maxNum
    temp = D_inv*(b_mat-(coefficient_mat-D)*x);
    if norm(temp-x) < epsilon
        flag = 1;
        x = temp;
        break;
    end
    x = temp;
end
%% 结果
text(0,0.95,'方程：','Color','red','FontSize',14);
str1 = '\left[\begin{array}';
% 正则表达式
str2 = @(X) regexprep(regexprep(mat2str(zeros(1,size(X,2))),...
    '[\s\[\]]',''),'0','c'); 
mat_start = @(X) [str1 '{' str2(X) '}'];
my_mat2str = @(mat) regexprep(regexprep(regexprep(mat2str(mat),...
    ' *','&'),';','\\\'),'[\[\]]','');
mat_end = '\end{array}\right]';
str5 = '\times X=';
% LaTex 输出
text('Interpreter','latex', 'String',['$$'...
    mat_start(coefficient_mat),my_mat2str(coefficient_mat),mat_end,...
    str5,mat_start(b_mat),my_mat2str(b_mat),mat_end '$$'],...
    'Position',[0.1 0.7],'FontSize',10);
text(0,0.5,'标准雅可比迭代结果：','Color','red','FontSize',14);
if flag == 1  %好的结果
    output_args.itrNum = k;
    output_args.result = x;
    fprintf('迭代次数：%d\n',output_args.itrNum);
%     disp(x(:,k));
    disp('解：');
    disp(x);
    text('Interpreter','latex', 'String',['$$X^*=' ...
        mat_start(x) my_mat2str(x) mat_end '$$'],...
        'Position',[0.1 0.3],'FontSize',10);
    text(0.1,0.1,['迭代次数：' num2str(output_args.itrNum)],...
    'Color','k','FontSize',10);
    return;
end
text(0.1,0.3,['迭代' num2str(maxNum) '次后未收敛'],...
    'Color','k','FontSize',10);
fprintf('迭代%d次后暂未收敛\n',maxNum);
end
  • 测试代码
1
2
3
4
xishu_mat = [20 8 -5 -5;7 9 -2 -10;-5 -1 9 -1;-7 1 -9 10];
b_mat = [-1 -5 2 -1];
init_X = zeros(size(b_mat));
out_Jacobi = cqj_Jacobi(init_X,xishu_mat,b_mat,1000);

如有任何疑问或建议欢迎下方留言 :-)