Python方程组获取雅克比矩阵和海塞矩阵

前面讲过当我们处理一个常微分方程组（一般对应于一个物理系统的求解）的时候，可以直接采用matlab中的ode函数比如ode15s和ode45，python的scipy中也有对应的函数像odeint和odebvp。但是，如果你想硬刚，自己去实现这些东西的话你就需要造很多的轮子，这里给出一个可以根据你输入的代数方程组获取雅克比矩阵和海塞矩阵的函数，并且我还特地将他做成了一个类哦，你可以直接进行实例化然后调用，非常方便，而且我还给出了测试函数。特别适合想硬刚底层算法的铁子。

主要的核心就是给出sympy中方程组的定义，以及sympy中将符号转化为可以计算的表达式的方式

# -*- coding: utf-8 -*-"""Created on 2022/1/25 10:45@Author: 123@Email: 985040320@qq.com@Describe: **加入文件描述, 要实现的功能, 注意事项等**"""class j_h_m:

    def __init__(self, funcs, vars):
        self.funcs = funcs
        self.vars = vars
        self.jac_m = sy.zeros(len(vars), len(vars))
        self.His_m = sy.zeros(len(vars), len(vars))

    # 获取雅克比矩阵
    def Ja(self):
        self.jac_m = self.funcs.jacobian(self.vars)
        return self.jac_m

    # 获取海塞矩阵
    def His(self):
        for i, fi in enumerate(self.funcs):
            for j, r in enumerate(self.vars):
                for k, s in enumerate(self.vars):
                    self.His_m[k, j] = sy.diff(sy.diff(fi, r), s)
        return self.His_mif __name__ == "__main__":
    import sympy as sy

    # 符号化
    vdot1, vdot2, y1, y2 = sy.symbols("vdot1, vdot2, y1, y2")
    # 构建方程组
    vdot1 = y2
    vdot2 = 1000 * (1 - y1 ** 2) * y2 - y1
    # 封装成sympy的符号矩阵
    funcs = sy.Matrix([vdot1, vdot2])
    args = sy.Matrix([y1, y2])
    # 实例化类并计算雅克比矩阵
    h = j_h_m(funcs, args)
    jac = h.Ja()
    his = h.His()
    # 将sympy的符号转化为函数表达式
    f = sy.lambdify([y1, y2], his, 'numpy')  # 通过这段话转化为可以计算的函数表达式
    print(jac)
    print(his)