数值积分法求解常微分方程组

import matplotlib
import numpy as np
import sympy
from scipy import integrate
from matplotlib import pyplot as plt



def f(xy, t, a, b, c, d):
    """x(t) 是猎物的数量， y(t) 是捕食者的数量。
    a是猎物的出生速度，d是猎物的死亡速度，b是捕食者消耗猎物的速度，c是捕食者种群的增长速度"""
    x, y = xy
    return [a*x - b*x*y, c*x*y -d*y]



if __name__ == '__main__':


    xy0 = 600, 400
    t = np.linspace(0, 50, 250)

    a,b = 0.2, 0.002
    c,d = 0.001, 0.7
    xy_t = integrate.odeint(f, xy0, t, args=(a, b, c, d))

    matplotlib.rcParams["font.sans-serif"] = ["SimHei"]
    figure, axes = plt.subplots(1, 2, figsize=(20, 10))

    axes[0].plot(t, xy_t[:, 0], "green", label="猎物")
    axes[0].plot(t, xy_t[:, 1], "red", label="捕食者")
    axes[0].set_xlabel("时间")
    axes[0].set_ylabel("动物数量")
    axes[0].legend()

    axes[1].plot(xy_t[:, 0], xy_t[:, 1],"blue")
    axes[1].set_xlabel("猎物数量")
    axes[1].set_ylabel("捕食者数量")
    axes[1].set_title("猎物捕食者模型 相空间")

    t = sympy.symbols('t')
    x = sympy.Function('x')
    y = sympy.Function('y')
    axes[0].set_title(f"${sympy.latex(sympy.Eq(x(t).diff(t), a* x(t) - b* x(t)*y(t)))}$ \n ${sympy.latex(sympy.Eq(y(t).diff(t), c* x(t)*y(t) - d* y(t)))}$" )
    plt.show()