python实现最速下降法求二元函数极值并绘制3D图像

目标

原函数： f(x) = 2x1x1 + x2x2 -2x1*x2 - 4x1 + 4
起始点 为（0，0）
精度 ：误差要求小于 10**（-3）

过程：

直接上代码看注释吧（）

import numpy as np
import math
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D


def f(x):
    return 2*x[0]**2 + x[1]**2 - 2*x[0]*x[1]-4*x[0]+4   

def step(x):      # 计算每次迭代的步长，H为该函数的hessian矩阵
    H = np.matrix([[4, -2], [-2, 2]])
    a = np.matrix(grad(x))
    b = np.matrix(grad(x).reshape(2, 1))    # array.reshape() 是把原列表按顺序变成指定的形状，本来是行向量形状是（2，1），经过处理后就变成了列向量形状是（2，1）

    return a*b/(a*np.matrix(H)*b)     


def grad(x):      #在该点的梯度，返回一个向量
    return np.array([4*x[0]-2*x[1]-4, 2*x[1]-2*x[0]])


def accuracy(x):   #计算在该点的精度
    return math.sqrt(x[0]*x[0]+x[1]*x[1])


if __name__=="__main__":
    x = np.array([0, 0])  # 定义初始点
    trace = []  # 用来记录迭代的轨迹
    deta = accuracy(grad(x)) # 精度
    while deta > 10**-3: # 一直循环，直到达到精度要求
        trace.append(np.append(x, f(x)))  
        t = float(step(x)) # step(x)为只有一个元素的矩阵（二维列表），所以先把他转化为浮点型
        x = x - t*grad(x)  # 更新点
        deta = accuracy(grad(x)) # 计算新的误差
    print(len(trace))  #看一下他迭代了多少次得到解
    
	# 一下是画图环节
	
    fig = plt.figure()
    ax = Axes3D(fig)

    X = np.arange(0, 3, 0.25)  # 先输出一下trace列表看看大致的范围再确定坐标，范围是（0，3）步长为 0.25

    Y = np.arange(0, 3, 0.25)

    X, Y = np.meshgrid(X, Y) # meshgrid是 网格 的意思，用来形成网格
    Z = 2*X**2 + Y**2 - 2*X*Y-4*X+4   # 求Z坐标的值   X，Y，Z都是列表
    plt.xlabel('x') #设置坐标轴的名称
    plt.ylabel('y')

    for p in trace:
        ax.scatter(*p)  # 显示轨迹点
    ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap='summer') #绘制函数在设置的坐标范围的图像 XYZ是坐标列表，后面两个是画图坐标的间隔，默认值就是1，越大跨度越大点越稀疏，cmap 是图的风格（颜色渐变等等），有很多自己去网上找
    plt.show()  # 显示

效果：
迭代次数为23

只看到两个点是因为被图覆盖了，如果只显示点，是这样的：

极值点为（2，2）