Matlab中使用共轭梯度 Jacobi G-S迭代法解决Ax=b

在上一篇我们通过对矩阵A的分解，将A分解为上三角矩阵与下三角矩阵的乘积，并通过简单地前代法和后代法得到x的值。
在本篇文章中，我们主要使用迭代法，通过不同的xk的递推关系，达到逼近x的效果，主要依靠的课本是北京大学徐树方老师等人编著的数值线性代数，下面开始正文部分。

设矩阵A b

我们是自己设出了一个比较复杂的A矩阵：

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

b为对应维数的单位列向量

共轭梯度法

function Conjugate_gradient_method(n)
%首先我们定义一个脚本matrix_Builder生成（n-1）^2维的方阵
%我们首先采用共轭梯度进行求解方程Ax=b
x0=ones((n-1)*(n-1),1);
tic
A = matrix_Builder(n);
b = ones((n-1)^2,1);
x=x0;
k=0;
r=b-A*x;
maxit=10000;
res=1e-6;
R=chol(A);
M=R'*R;
while sqrt(r'*r)>res*norm(b,2) && k<maxit
    z=inv(M)*r;
    k=k+1;
    if k==1
        p=z;
        rho=r'*z;
    else
        rhoo=rho;
        rho=r'*z;
        beta=rho/rhoo;
        p=z+beta*p;
    end
    w=A*p;
    alpha=rho/(p'*w);
    x=x+alpha*p;
    r=r-alpha*w;
end
delta=norm(x-inv(A)*b);
toc
end

Jacobi迭代

function  Jacobi(n)
%首先我们定义一个脚本matrix_Builder生成（n-1）^2维的方阵
%我们采用jacobi进行求解方程Ax=b
x0 = ones([(n-1)^2,1]);
A = matrix_Builder(n);b = randn([(n-1)*(n-1),1]);
D = diag(diag(A));U = -(triu(A)-D);L = -(tril(A)-D);
tic
B = inv(D)*(L+U);g = inv(D)*b;
xk = B*x0+g;k=0;
while norm(xk-x0)>0.000001
    x0=xk;
    xk = B*x0+g;
    k=k+1;
    %fprintf('在%d轮',k);
    %fprintf('x(k+1)-x(k)为%f\n',norm(xk-x0));
end
toc

Gauss-Seidel迭代法

function  Gauss_Seidel(n)
%首先我们定义一个脚本matrix_Builder生成（n-1）^2维的方阵
%我们Gauss_Seidel采用进行求解方程Ax=b
x0 = ones([(n-1)^2,1]);
A = matrix_Builder(n);b = randn([(n-1)*(n-1),1]);
D = diag(diag(A));U = -(triu(A)-D);L = -(tril(A)-D);
tic
B = inv(D-L)*U;g = inv(D-L)*b;
xk = B*x0+g;k=0;
while norm(xk-x0)>0.000001
    x0=xk;
    xk = B*x0+g;
    k=k+1;
    %fprintf('在%d轮，x(k+1)-x(k)为%f\n',k,norm(xk-x0));
end
toc

在实际编程中，笔者认为首先需要将算法步骤写出来，即伪代码，这样我们在实际写这种有点复杂的过程时，便于对着算法的流程修改bug，非常方便~