python实现共轭梯度算法（含误差与运算次数的折线图）

西安交通大学《数值分析》第三章课后题3.2

import numpy as np
import math
from math import log
import matplotlib.pyplot as plt
def generate_matrix(n):
    # 使用对角矩阵相加得到三对角矩阵A
    array_a = np.diag([-2] * n)
    array = np.diag([1] * (n-1))
    a = np.zeros((n-1))
    b = np.zeros(n)
    array_b = np.insert(array, 0, values=a, axis=0)# 添加行
    array_b = np.insert(array_b, (n-1), values=b, axis=1)# 添加列
    array_c = np.insert(array, (n-1), values=a, axis=0)
    array_c = np.insert(array_c, 0, values=b, axis=1)
    matrix_A = array_a + array_b + array_c

    return matrix_A

def Norm(r):
    sum = 0
    for i in r:
        sum += i*i
    norm_r = '%.15f' % math.sqrt(sum)
    print('norm',norm_r)
    return norm_r

def Conjugate_Gradient(ri,di,A,xi,b):
    print('---------------------------------------------------------')
    global beta
    beta = (-1) * (np.dot(ri.T,np.dot(A,di)) / np.dot(di.T,np.dot(A,di)))
    print('beta',beta)
    print('d',di)
    print('r',ri)
    global d
    d = ri + beta*di
    print('d',d)
    global a
    a = np.dot(ri.T,d) / np.dot(d.T,np.dot(A,d))
    print('a',a)
    global x
    x = xi + a*d
    print('x',x)
    global r
    r = b - np.dot(A,x)
    print('r',r)
    # 计算r的2-范数
    return Norm(r)


if __name__ == '__main__':
    # 设置系数矩阵A的阶数
    n = 50
    # 生成系数矩阵A
    A = generate_matrix(n)
    print('A',A)
    # 生成矩阵b
    b = np.zeros(n)
    b[0] = -1
    b[(n-1)] = -1
    print('b',b)
    # 设置初始值
    global x
    x = np.zeros(n)
    # 第一次迭代
    global r
    r = b - np.dot(A,x)
    print('r',r)
    global d
    d = r
    global a
    a = np.dot(r.T,d) / np.dot(d.T,np.dot(A,d))
    print('a',a)
    x = x + a*d
    print('x',x)
    r = b - np.dot(A,x)
    print('r',r)
    global norm_r
    norm_r = Norm(r)
    log_norm = log(float(norm_r), 10)  # 画图 以10为底取对数
    list_norm = [] # 画图
    list_norm.append(log_norm) # 画图 加入第一次计算的对数范数
    # 设置误差上限alpha
    alpha = 1e-8
    count = 1 # 迭代次数
    list_count = [] # 画图
    list_count.append(count) # 画图 加入第一次迭代计数
    while float(norm_r) >= alpha:
        # 第2~N次迭代
        count += 1
        norm_r = Conjugate_Gradient(r,d,A,x,b)
        log_norm = log(float(norm_r),10) # 画图 以10为底取对数
        list_norm.append(log_norm) # 画图 加入后续的对数范数
        list_count.append(count) # 画图
        # print(norm_r)
        # print('r',r)
        # print('d',d)
        # print('x',x)
    print('x*',x)
    print('count',count)
    plt.plot(list_count,list_norm)
    plt.show()

结果输出：答案为全是1

A [[-2  1  0 ...  0  0  0]
 [ 1 -2  1 ...  0  0  0]
 [ 0  1 -2 ...  0  0  0]
 ...
 [ 0  0  0 ... -2  1  0]
 [ 0  0  0 ...  1 -2  1]
 [ 0  0  0 ...  0  1 -2]]
b [-1.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.
  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.
  1.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.
  2.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.
  3.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.
  4.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.
  5.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.
  6.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.
  7.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.
  8.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.
  9.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.
  10. -1.]
  
无数次循环后
norm 0.000000000000195
x* [1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
 1. 1.]
count 25

另一篇笔记 快速下降法需要计算快8k次，我人没了。