算法---使用牛顿迭代法求极值并绘制曲线

import numpy as np
import matplotlib.pyplot as plt

# 定义函数 jacobian，用于计算目标函数的雅可比矩阵
def jacobian(x):
    return np.array([-400 * x[0] * (x[1] - x[0] ** 2) - 2 * (1 - x[0]), 200 * (x[1] - x[0] ** 2)])


# 定义函数 hessian，用于计算目标函数的海森矩阵
def hessian(x):
    return np.array([[-400 * (x[1] - 3 * x[0] ** 2) + 2, -400 * x[0]],
                     [-400 * x[0], 200]])


# 生成坐标网格，用于后续绘制函数等高线图
# 在 x1 方向上，从 -1.5 到 1.5（包含），步长为 0.05 生成一组数值
x1 = np.arange(-1.5, 1.5 + 0.05, 0.05)
# 在 x2 方向上，从 -3.5 到 2（包含），步长为 0.05 生成一组数值
x2 = np.arange(-3.5, 2 + 0.05, 0.05)
# 使用 meshgrid 函数将 x1 和 x2 这两组一维数组转换为二维坐标网格，方便后续计算和绘图
[x1, x2] = np.meshgrid(x1, x2)

# 定义目标函数
# 这里的目标函数是一个二元函数，形式为 100 * (x2 - x1 ** 2) ** 2 + (1 - x1) ** 2
f = 100 * (x2 - x1 ** 2) ** 2 + (1 - x1) ** 2
# 绘制目标函数的等高线图，等高线数量为 20 条，通过等高线图可以直观地看到函数的形状和极值情况等
plt.contour(x1, x2, f, 20)


# 定义函数 newton，使用牛顿法对给定的初始点进行迭代优化
def newton(x0):
    print('初始点为:')
    print(x0, '\n')
    # 用于存储迭代过程中的点，初始化一个较大的二维数组，形状为 (2, 10**3)，
    # 这里假设最多迭代 1000 次，实际迭代结束后会截取用到的部分
    W = np.zeros((2, 10 ** 3))
    i = 1  # 迭代次数计数器，初始化为 1，因为第一次迭代前已经给定了初始点
    imax = 1000  # 最大迭代次数限制，防止迭代不收敛而陷入无限循环
    W[:, 0] = x0  # 将初始点赋值给存储迭代点的数组 W 的第一列
    x = x0  # 将初始点赋值给当前迭代点变量 x
    delta = 1  # 用于衡量相邻两次迭代点之间的距离变化，初始化为 1
    alpha = 1  # 步长，这里暂时固定为 1，实际应用中可以根据情况进行调整
    while i < imax and delta > 10 ** (-5):
        # 计算牛顿方向，通过求解海森矩阵的逆与雅可比矩阵的乘积得到下降方向
        p = -np.dot(np.linalg.inv(hessian(x)), jacobian(x))
        x0 = x  # 记录当前迭代点，用于后续计算距离变化
        x = x + alpha * p  # 根据牛顿方向和步长更新当前迭代点
        W[:, i] = x  # 将更新后的迭代点存入存储迭代点的数组 W 中
        delta = np.sum((x - x0) ** 2)  # 计算当前迭代点与上一次迭代点之间距离的平方和，衡量迭代收敛情况
        print('第', i, '次迭代结果:')
        print(x, '\n')
        i = i + 1  # 迭代次数加 1
    # 截取实际用到的迭代点，去除初始化时多余的部分，只保留从初始点到最后一次有效迭代的点
    W = W[:, 0: i]
    return W


# 设置初始点，这里初始点为一个包含两个元素的一维数组，对应函数自变量的初始取值
x0 = np.array([-1.2, 1])
W = newton(x0)
# 绘制牛顿法迭代路径，将迭代过程中的点以绿色星号（'g*'）的形式绘制在之前绘制的等高线图上，
# 这样可以直观地看到牛顿法的迭代走向以及最终是否收敛到函数的极值点附近
plt.plot(W[0, :], W[1, :], 'g*')
# 显示图形，将绘制好的等高线图和迭代路径展示出来
plt.show()
返回结果：