上机实验四：共轭梯度法程序设计（MATLAB）

共轭梯度法是求解线性方程组的一种有效方法，MATLAB中提供了相关的函数可以使用。这里我们将介绍如何手动编写共轭梯度法程序。

首先，我们需要明确共轭梯度法的基本思路。共轭梯度法是一种迭代算法，其核心思想是通过利用前一次迭代得到的信息来加速收敛。具体来说，对于一个对称正定的系数矩阵A和一个右端向量b，共轭梯度法的基本步骤如下：

1. 初始化迭代：取一个初始向量x0和残差向量r0=b-Ax0，设置迭代次数k=0。

2. 计算搜索方向：对于第k次迭代，计算搜索方向p_k=r_k，如果k=0，则p_0=r_0；否则，p_k=r_k+beta_k*p_{k-1}，其中beta_k是一个系数，定义为beta_k=(r_k'*r_k)/(r_{k-1}'*r_{k-1})。

3. 计算步长：沿着搜索方向p_k移动一个步长alpha_k，使得x_{k+1}=x_k+alpha_k*p_k。

4. 更新残差：计算新的残差r_{k+1}=b-Ax_{k+1}。

5. 检查停止条件：如果满足停止条件（如残差的范数小于某个阈值），则停止迭代；否则，令k=k+1，返回步骤2。

基于以上思路，我们可以编写共轭梯度法的MATLAB程序。下面是一个简单的例子：

function [x, k] = conjgrad(A, b, x0, tol, maxiter)
% 输入参数：
% A: 系数矩阵
% b: 右端向量
% x0: 初始向量
% tol: 残差的阈值
% maxiter: 最大迭代次数

% 输出参数：
% x: 迭代结束时的解向量
% k: 实际迭代次数

r = b - A * x0;      % 初始残差
p = r;               % 初始搜索方向
x = x0;              % 初始解向量
k = 0;               % 初始迭代次数

while norm(r) > tol && k < maxiter
    alpha = (r'*r) / (p' * A * p);     % 计算步长
    x = x + alpha * p;                 % 更新解向量
    r_new = r - alpha * A * p;         % 计算新的残差
    beta = (r_new' * r_new) / (r' * r);% 计算beta
    p = r_new + beta * p;              % 计算新的搜索方向
    r = r_new;                         % 更新残差
    k = k + 1;                         % 更新迭代次数
end
end

这个程序定义了一个名为conjgrad的函数，输入参数包括系数矩阵A、右端向量b、初始向量x0、残差阈值tol和最大迭代次数maxiter。输出参数包括迭代结束时的解向量x和实际迭代次数k。

在函数内部，我们首先计算了初始残差r、初始搜索方向p和初始解向量x。然后进入迭代循环，直到满足停止条件为止。在每次迭代中，我们先计算了步长alpha，然后更新了解向量x。接着计算了新的残差r_new和新的搜索方向p，最后更新了残差和迭代次数。在迭代结束后，函数返回最终的解向量和实际迭代次数。

使用这个函数非常简单，只需要提供系数矩阵A、右端向量b和初始向量x0，以及其他参数即可。例如：

A = [4 -1 0; -1 4 -1; 0 -1 4];
b = [10; 10; 10];
x0 = [0; 0; 0];
tol = 1e-6;
maxiter = 100;

[x, k] = conjgrad(A, b, x0, tol, maxiter);

disp(['迭代次数：', num2str(k)]);
disp(['解向量：', num2str(x')]);

这个例子是求解线性方程组Ax=b，其中系数矩阵A是一个对称正定矩阵。我们设置了初始向量x0=[0; 0; 0]，残差阈值tol=1e-6和最大迭代次数maxiter=100，然后调用conjgrad函数求解。最后输出了迭代次数和解向量。