python实现共轭梯度法

以下是Python实现共轭梯度法的示例代码：

import numpy as np

def conjugate_gradient(A, b, x0, tol=1e-6, max_iter=100):
    # A: 系数矩阵
    # b: 常数向量
    # x0: 初始解向量
    # tol: 迭代停止的残差阈值
    # max_iter: 最大迭代次数
    
    x = x0
    r = b - A @ x # 残差
    p = r # 初始搜索方向
    for k in range(max_iter):
        alpha = np.dot(r, r) / np.dot(p, A @ p) # 步长
        x = x + alpha * p # 更新解向量
        r_new = r - alpha * A @ p # 更新残差
        if np.linalg.norm(r_new) < tol: # 判断是否满足停止条件
            break
        beta = np.dot(r_new, r_new) / np.dot(r, r) # 计算beta
        p = r_new + beta * p # 更新搜索方向
        r = r_new # 更新残差
    
    return x

使用示例：

A = np.array([[4, 1], [1, 3]])
b = np.array([1, 2])
x0 = np.array([0, 0])
x = conjugate_gradient(A, b, x0)
print(x)

输出：

[ 0.09090909  0.63636364]

其中，A和b分别表示线性方程组Ax=b的系数矩阵和常数向量，x0表示初始解向量，tol表示迭代停止的残差阈值，max_iter表示最大迭代次数。