发篇博客纪念一下,折磨了我一个小时。。。
利用逆Broyden方法求方程组的数值解,迭代到某误差范围内结束,精确解x=(0.5,1)转秩
import numpy as np
import pprint
f_l = []
h_l = []
x_l = []
for i in range(12):
if i == 0:
x1 = 0.75
x2 = 0.75
# print(x1, x2)
x_l.append([x1, x2])
A = np.array([[8 * x1 + 2 * x2, 2 * x2 + 2 * x1 - 1], [4 * x1 + 3 * x2, 3 * x1 + 2 * x2]]) # F(X)的导数 A0
H = np.linalg.inv(A) # H0
h_l.append(H)
X = np.array([[x1], [x2]]) # Xk
F = np.array([[x1 * x1 * 4 + x2 * x2 + 2 * x1 * x2 - x2 - 2], [2 * x1 * x1 + 3 * x1 * x2 + x2 * x2 - 3]]) # F_Xk
f_l.append(F)
if i >= 1:
Y = f_l[-1] - f_l[-2]
H = h_l[-1] + (
(R - np.dot(h_l[-1], Y)) * np.dot(R.T, h_l[-1]) / (np.dot(np.dot(R.T, h_l[-1]), Y)))
h_l.append(H)
X_new = X - np.dot(H, F) # x1
R = X_new - X # r0
x1 = float(X_new[0])
x2 = float(X_new[1])
x_l.append([x1, x2])
# print(x1, x2)
pprint.pprint(x_l)
for i in range(10):
a = max(abs(x_l[i + 1][0] - x_l[i][0]), abs(x_l[i + 1][1] - x_l[i][1]))
print(a, bool(a <= 0.5 / 100000), i + 1)