最优化-共轭梯度法
import numpy as np
def CG(A,b):
n = b.shape[0]
xs = []
rs = []
ps = []
alphas = []
x0 = np.array([2,1])
xs.append(x0)
r0 = b - np.dot(A,x0)
rs.append(r0)
p0 = r0
ps.append(p0)
alpha0 = p0.dot(p0)/p0.dot(A).dot(p0)
alphas.append(alpha0)
print(rs)
print(alphas)
print(ps)
for i in range(n):
r = rs[i] - alphas[i] * A.dot(ps[i])
rs.append(r)
beta = np.dot(r,r)/(rs[i].dot(rs[i]))
alpha = ps[i].dot(rs[i])/(ps[i]).dot(A).dot(ps[i])
alphas.append(alpha)
x = xs[i] + alpha * ps[i]
xs.append(x)
p = r + beta * ps[i]
ps.append(p)
return xs
A = np.array([[4,1],[1,3]])
b = np.array([1,2])
c1 = np.linalg.inv(A).dot(b)
print("the math sol is",c1)
c2 = CG(A,b)
print("the numerical sol is",c2)