数值计算方法——Jacobi迭代法和G-S迭代法

没完成的代码，仅做储存

function [L,D,U]=LU(A)
n=length(A);
L=zeros(n,n);%矩阵初始化
U=zeros(n,n);
for i=2:n
    for j=1:i-1
        L(i,j)=A(i,j);
    end
end
for i=1:n-1
    for j=i+1:n
        U(i,j)=A(i,j);
    end
end
D=diag(diag(A));

A=[1 -2 2
    -1 1 -2
    -2 -2 1];
[L,D,U]=LU(A);

disp(-L);  %输出L
disp(D);  %输出D
disp(-U); %输出U
[x,y]=Jacobi(A);

function [x,y]=Jacobi(A,b)
n=length(A);
[L,D,U]=LU(A); %调用函数
G=D\(L+U);
[x,y]=eig(G);
% if max(abs(eig(G)))>1
%     disp('不收敛,请重新调用函数')
% end

矩阵常用操作

A = [1,-1,2,3;5,1,-4,2;3,0,5,2;11,15,0,9];
B = [0.43,43,2;-8.9,4,21];
%(1)方阵
diag(A) %主对角元素
triu(A) %上三角
tril(A) %下三角
inv(A)  %矩阵的逆
det(A)  %行列式的值
rank(A) %矩阵的秩
norm(A) %矩阵的范数
cond(A) %矩阵的条件数
trace(A)%矩阵的迹
%(2)非方阵
diag(B) %主对角元素
triu(B) %上三角
tril(B) %下三角
pinv(B)  %矩阵的逆
%B不是方阵，不存在行列式的值
rank(B) %矩阵的秩
norm(B) %矩阵的范数
cond(B) %矩阵的条件数
sum(diag(B))%矩阵的迹

Jacobi迭代实现

文件名：Jacobi_solve

function [x,error,iter]=Jacobi_solve(A,b,epsilon,max_iter)
n=length(b);
x=zeros(n,1);
error=1;
iter=0;    %输出的值初始化
B=zeros(1,n);%用于对角元素置0
L=-tril(A);%取矩阵的下三角
U=-triu(A);%取矩阵的上三角
D=diag(diag(A));%取矩阵的对角元素
L(logical(eye(size(L))))=B;%对角元素置0
U(logical(eye(size(U))))=B;
disp('k x^k err');
disp(D);
disp(L);
disp(U);
while error>epsilon && iter<max_iter
    for i=1:max_iter
        y=x;%提前将旧值赋值赋给另一个量
        x=D\(L+U)*x+D\b;%x应该是一个向量，进行迭代运算
        error=norm(x-y);%误差为两个量的范数
        if error>epsilon
            iter=iter+1;%输出达到精度的最小迭代次数
            break
        end
    end
end

文件名：Jacobi_test

A=[8 -3 2
    4 11 -1
    6 3 12];
b=[20;33;36];
epsilon=10e-6;
max_iter=15;
[x,error,iter]=Jacobi_solve(A,b,epsilon,max_iter);
disp(x);
disp('最大迭代次数下的误差：');
disp(error);
disp('最小迭代次数：');
disp(iter+1);

输出结果：

G-S迭代实现

function [x,error,iter]=Jacobi_solve(A,b,epsilon,max_iter)
n=length(b);
x=zeros(n,1);
error=1;
iter=0;    %输出的值初始化
B=zeros(1,n);%用于对角元素置0
L=-tril(A);%取矩阵的下三角
U=-triu(A);%取矩阵的上三角
D=diag(diag(A));%取矩阵的对角元素
L(logical(eye(size(L))))=B;%对角元素置0
U(logical(eye(size(U))))=B;
disp('k x^k err');
disp(D);
disp(L);
disp(U);
while error>epsilon && iter<max_iter
    for i=1:max_iter
        y=x;%提前将旧值赋值赋给另一个量
        x=(D-L)\U*x+(D-L)\b;%G-S迭代就是Jacobi直接修改此行就好
        error=norm(x-y);%误差为两个量的范数
        if error>epsilon
            iter=iter+1;%输出达到精度的最小迭代次数
            break
        end
    end
end

A=[8 -3 2
    4 11 -1
    6 3 12];
b=[20;33;36];
epsilon=10e-6;
max_iter=15;
[x,error,iter]=Jacobi_solve(A,b,epsilon,max_iter);
disp(x);
disp('最大迭代次数下的误差：');
disp(error);
disp('最小迭代次数：');
disp(iter+1);

输出结果 ：