Cholesky分解法 | matlab

function [X]=Mchol(A,b) 
%Cholesky分解法
[N, N] = size(A); 
X = zeros(N, 1); 
Y = zeros(N, 1); 
for i = 1:N 
    A(i, i) = sqrt(A(i, i)-A(i, 1:i-1) * A(i, 1:i-1)'); 
    if A(i, i) == 0 
        fprintf('矩阵是奇异的')
        break 
    end 
    for j = i+1:N
        A(j, i) = (A(j, i) - A(j, 1:i-1) * A(i, 1:i-1)') / A(i, i); 
    end 
end 
A 
b 
%前代法
for j = 1:N 
    Y(j) = (b(j) - A(j, 1:j-1) * Y(1:j-1)) / A(j, j); 
end 
Y 
% 
A=A' 
for k = N:-1:1 
    X(k) = (Y(k) - A(k, k+1:N) * X(k+1:N)) / A(k, k); 
end 

A = [4, 2, 4; 2, 10, -1; 4, -1, 6]
b = [4; 17; 0]
m_chol(A, b)

%=> x = [-1, 2, 1]