三对角矩阵的类雅可比迭代程序


三对角矩阵的类雅可比迭代算法解线性方程组

  • 函数实现
function [ output_args ] = cqj_QuasiJacobi(init_X,coefficient_mat,b_mat,maxNum,epsilon)
% 类雅可比迭代法 只不过对角矩阵D由三对角矩阵代替
% 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

%% 三对角线
a = diag(coefficient_mat,-1)';
b = diag(coefficient_mat)';
c = diag(coefficient_mat,1)';
D = diag(a,-1)+diag(b)+diag(c,1);

% 采用冉瑞生等人《三对角矩阵求逆的算法》中的算法
[~,D_inv] = cqj_InverseMatrixOfTridiagonalMatrices(a,b,c);
if isempty(D_inv)
    disp('系数矩阵不可逆');
    return;
end
%% 开始迭代
% 下面代码段也行 但是对内存的消耗比较大 
% 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
  • 测试程序
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_QuasiJacobi = cqj_QuasiJacobi(init_X,xishu_mat,b_mat,1000);

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


评论