本文参考链接：牛顿法实现
偏导数求解

import numpy as np
from sympy import *


#获取变量
def get_vars(*args):
    arg_str = ''
    i=0
    for arg in args:
        i += 1
        if i < len(args):
            arg_str += arg +  ','
        else:
            arg_str += arg
    return arg_str

#获取hessian矩阵
def get_hessian(a,b):
    hessian = zeros(2, 2)  
    for i,fi in enumerate(f):
        for j,r in enumerate(vars):
            for k, s in enumerate(vars):
                hessian[j,k] = diff(diff(fi, r),s).subs({vars[0]:a, vars[1]:b})
    return hessian

#牛顿法迭代
def newton(max_step, x_init):
    i = 1
    while i < max_step:
        if i == 1: 
            #第一次迭代
            grandient = np.array([diff(f1,vars[0]).subs({vars[0]:x_init[0], vars[1]:x_init[1]}), 
                                 diff(f1,vars[1]).subs({vars[0]:x_init[0], vars[1]:x_init[1]})]) 
            hessian = get_hessian(x_init[0], x_init[1])
            #此处需要对hessian进行求逆，需要先将其转为类型为float的矩阵再进行求解，否则可能报错
            new_ab = x_init - np.matmul(np.linalg.inv(np.mat(hessian,dtype='float')), grandient)
        else: 
            grandient = np.array([diff(f1,vars[0]).subs({vars[0]:new_ab[0,0], vars[1]:new_ab[0,1]}), 
                             diff(f1,vars[1]).subs({vars[0]:new_ab[0,0], vars[1]:new_ab[0,1]})])
            hessian = get_hessian(new_ab[0,0], new_ab[0,1])
            new_ab = np.array(new_ab - np.matmul(np.linalg.inv(np.mat(hessian,dtype='float')), grandient))
        print('迭代第%d次：%.5f %.5f' %(i, new_ab[0,0], new_ab[0,1]))
        i = i + 1
    return new_ab

var = get_vars('a', 'b')
vars = symbols(var)
f1 = vars[0]**4 + vars[1]**4 + vars[0]*vars[1]
f = sympify([str(f1)])
x_init = np.array([30,30])
max_step = 30
newton(max_step, x_init)

输出结果

函数图像

从输出结果可以看出最终收敛至（0,0）,即f在（0,0）处取得极小值。