import numpy as np
import matplotlib.pyplot as plt
from sympy import *
# 定义符号
x1, x2 = symbols('x1, x2')
# 定义所求函数
f = 0.2*x1**2 + x2**2
# 求解梯度值
def get_grad(f, X):
# 计算一阶导数
f1 = diff(f, x1)
f2 = diff(f, x2)
X = X.tolist()
X1 = X[0][0]
X2 = X[1][0]
# 代入具体数值计算
grad = np.array([ [f1.subs({x1: X1, x2: X2}).evalf()],
[f2.subs({x1: X1, x2: X2}).evalf()] ])
return grad
# 求解Hession矩阵
def get_hess(f, X):
# 计算二次偏导
f1 = diff(f, x1)
f2 = diff(f, x2)
f11 = diff(f,x1,2)
f22 = diff(f,x2,2)
f12 = diff(f1,x2)
f21 = diff(f2,x1)
# 计算具体数值计算
hess = np.array([[f11.subs([(x1,X[0]), (x2,X[1])]),
f12.subs([(x1,X[0]), (x2,X[1])])],
[f21.subs([(x1,X[0]), (x2,X[1])]),
f22.subs([(x1,X[0]), (x2,X[1])])]])
# 转换数值类型为了后续求逆矩阵
hess = np.array(hess, dtype = 'float')
return hess
# 牛顿迭代
def newton_iter(X, epsilon, max_iter):
print('初始值','x1=',X[0][0],'x2=',X[1][0])
count = 0
while count< max_iter:
grad = get_grad(f, X)
grad_1_value = grad[0][0]
grad_1_value = grad[1][0]
if abs(grad_1_value) + abs(grad_1_value) >= epsilon:
hess = get_hess(f, X0)
# 得到Hession矩阵的逆
hess_inv = np.linalg.inv(hess)
# 牛顿迭代公式!!!!
X = X - np.dot(hess_inv, grad)
count += 1
print('第',count,'次迭代:','x1=',X[0][0],'x2=',X[1][0])
else:
break
print('迭代次数为:',count)
# 设置初始点
X0 = np.array([[1],[1]])
epsilon = 0.00001
max_iter = 50
newton_iter(X0, epsilon, max_iter)
def f(x, y):
return 0.2*(x**2)+y**2
X = np.array([1, 0])
Y = np.array([1, 0])
x_values = np.linspace(-1, 1, 500)
y_values = np.linspace(-1, 1, 500)
xx, yy = np.meshgrid(x_values, y_values)
fig = plt.figure(figsize=(8, 8))
contour_line = plt.contour(xx, yy, f(xx, yy), levels=20, cmap=plt.cm.gray)
plt.clabel(contour_line, inline=1, fontsize=10)
plt.plot(X, Y, color='r', linestyle='-', marker='*', linewidth =2.0)
plt.show()
07-06
9698
![](https://csdnimg.cn/release/blogv2/dist/pc/img/readCountWhite.png)
03-05