前面講過當(dāng)我們處理一個(gè)常微分方程組（一般對(duì)應(yīng)于一個(gè)物理系統(tǒng)的求解）的時(shí)候，可以直接采用matlab中的ode函數(shù)比如ode15s和ode45，python的scipy中也有對(duì)應(yīng)的函數(shù)像odeint和odebvp。但是，如果你想硬剛，自己去實(shí)現(xiàn)這些東西的話你就需要造很多的輪子，這里給出一個(gè)可以根據(jù)你輸入的代數(shù)方程組獲取雅克比矩陣和海塞矩陣的函數(shù)，并且我還特地將他做成了一個(gè)類哦，你可以直接進(jìn)行實(shí)例化然后調(diào)用，非常方便，而且我還給出了測(cè)試函數(shù)。特別適合想硬剛底層算法的鐵子。

主要的核心就是給出sympy中方程組的定義，以及sympy中將符號(hào)轉(zhuǎn)化為可以計(jì)算的表達(dá)式的方式

# -*- coding: utf-8 -*-"""Created on 2022/1/25 10:45@Author: 123@Email: 985040320@qq.com@Describe: **加入文件描述, 要實(shí)現(xiàn)的功能, 注意事項(xiàng)等**"""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

    # 符號(hào)化
    vdot1, vdot2, y1, y2 = sy.symbols("vdot1, vdot2, y1, y2")
    # 構(gòu)建方程組
    vdot1 = y2
    vdot2 = 1000 * (1 - y1 ** 2) * y2 - y1
    # 封裝成sympy的符號(hào)矩陣
    funcs = sy.Matrix([vdot1, vdot2])
    args = sy.Matrix([y1, y2])
    # 實(shí)例化類并計(jì)算雅克比矩陣
    h = j_h_m(funcs, args)
    jac = h.Ja()
    his = h.His()
    # 將sympy的符號(hào)轉(zhuǎn)化為函數(shù)表達(dá)式
    f = sy.lambdify([y1, y2], his, 'numpy')  # 通過這段話轉(zhuǎn)化為可以計(jì)算的函數(shù)表達(dá)式
    print(jac)
    print(his)

喜歡的朋友可以給個(gè)關(guān)注或者聯(lián)系我