通过前代法后代法列主元Cholesky求解Ax=b

首先设出要求解的矩阵

function A = matrix_Builder(n)
%对于一个n，我们按课本120页要求生成(n-1)^2*(n-1)^2的矩阵
A = diag(repmat([4], 1, (n-1)^2))+diag(repmat([-1], 1, n*(n-2)), 1)+diag(repmat([-1], 1, n*(n-2)), -1)-diag(repmat([1],1,(n-1)*(n-2)),1-n)-diag(repmat([1],1,(n-1)*(n-2)),n-1);
for i = 1:n-2
    A((n-1)*i,(n-1)*i+1)=0;
    A((n-1)*i+1,(n-1)*i)=0;
end

简述前代后代法

前代法后代法主要是用来求解A为上三角或者下三角的情况，我们首先通过列主元消去或者是Cholesky分解将一个一般的矩阵转化为两个上下三角矩阵的乘积，然后通过上下三角线性方程组的求解算法，即可解出解

前代法

function b = Predecessor(L,b)
n = length(b);
for j=1:n
    b(j) = b(j)/L(j,j);
    b(j+1:n) = b(j+1:n) - b(j)*L(j+1:n,j);
end
b(n) = b(n)/L(n,n);

后代法

function y = Descendant(U,y)
n = length(y);
for j =n:-1:2
    y(j) = y(j)/U(j,j);
    y(1:j-1) = y(1:j-1)-y(j)*U(1:j-1,j);
end
y(1) = y(1)/U(1,1);

矩阵分解

我们给出的是通过列主元消去和Cholesky分解方法对矩阵进行分解：

列主元消去

function  Column_principal_elimination(n)
%首先我们定义一个脚本matrix_Builder生成（n-1）^2维的方阵
%我们首先采用列主元消去进行求解方程Ax=b,计算量为2/3*n^3
A = matrix_Builder(n);
b = randn([(n-1)*(n-1),1]);tic
P = eye((n-1)^2);
for i = 1:(n-1)^2
    %首先确定最大值的位置,row_max为此时的列最大值，Line_indicator为最大值的行指标
    A_remaining = A(i:(n-1)^2,i:(n-1)^2);
    A_row_i = A(:,i);
    row_max = max(abs(A_row_i(i:(n-1)^2)));
    Line_indicator = find(A_row_i==row_max);
    %然后我们设置容器变量temp，并利用Line_indicator进行交换行
    temp = A(Line_indicator,:);
    A(Line_indicator,:)=A(i,:);
    A(i,:)=temp;
    temp = P(Line_indicator,:);
    P(Line_indicator,:)=P(i,:);
    P(i,:)=temp;
    if abs(A(i,i))
        A(i:(n-1),i) = A(i+1:n,i)/A(i,i);
        A(i+1:n,i+1:n) = A(i+1:n,i+1:n)-A(i+1:n,i)*A(i,i+1:n);
    else
        break
    end
end
%这一算法将LU分别存储在A的上三角和下三角部分，最后P即为置换矩阵.
%计算Ly = Pb，Ux = y；
L = tril(A,-1)+eye((n-1)^2);
U = triu(A);
y = Descendant(L,P*b);
x = Predecessor(U,y);
toc

Cholesky分解

function cholesky(n)
%首先我们定义一个脚本matrix_Builder生成（n-1）^2维的方阵
%我们首先采用cholesky进行求解方程Ax=b,计算量为1/3*n^3
A = matrix_Builder(n);
b = randn([(n-1)*(n-1),1]);
temp=A\b;
tic
v = zeros([1,(n-1)^2]) ;
n=(n-1)^2;
for j = 1:n
    for i =1:j-1
        v(i) = A(j,i)*A(i,i);
    end
    A(j,j) = A(j,j)-A(j,1:j-1)*v(1:j-1)';
    A(j+1:n,j) = (A(j+1:n,j)-A(j+1:n,1:j-1)*v(1:j-1)')/A(j,j);
end
%最后我们解Ly=b,DL'x=y
L = tril(A,-1)+eye(n);
D = diag(diag(A));
y = Predecessor(L,b);
x = Descendant(D*L',y);norm(x-temp);
toc

理论分析结果为：
列主元消去计算量为2/3n^3
Cholesky分解计算量为1/3n^3
实际计算通过计算时间得到同样结果！
于是就结束了！