数值积分法求解常微分方程组
求解 方程组
2023-06-13 09:15:27 时间
Scipy 的 integrate 模块的 odeint 函数也可以用来以数值积分法求解常微分方程组。下面的代码以 猎物-捕食者模型为例讲解其用法。
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()
图中可以看出,狼的数量快速增长的时候,羊的数量急速下降。
相关文章
- matlab求解不定方程组_matlab解参数方程组
- matlab求解时滞微分方程_matlab延迟环节传递函数
- 小站工具|求解课题万能公式“A基因通过B通路调控C疾病的D功能”与模块使用教程~
- 四阶龙格库塔(Runge-Kutta)求解微分方程-多种编程语言
- 论文拾萃 | PISTS算法求解obnoxious p-median problem (附Python代码)
- 【运筹学】线性规划数学模型 ( 线性规划求解 | 根据非基变量的解得到基变量解 | 基解 | 基可行解 | 可行基 )
- 【组合数学】递推方程 ( 有重根递推方程求解问题 | 问题提出 )
- 【组合数学】生成函数 ( 使用生成函数求解不定方程解个数 )
- 【组合数学】指数生成函数 ( 指数生成函数求解多重集排列示例 2 )
- Lingo V18中文版电脑安装,Lingo优化求解软件下载安装教程
- Lingo优化求解器软件V18激活版电脑下载安装,Lingo软件下载激活
- 求解Oracle时间差:求出两个日期之间的天数(oracle时间差天数)
- 高效求解MySQL问题,瞬间秒杀上机答案(mysql上机答案)
- Oracle 中两数相除快速求解你要的结果(oracle 中两数相除)
- 差异求解Oracle中两日期月份之差(oracle两日期月)
- oracle一列的总数求解(oracle 一列总数)
- 用贪心法求解背包问题的解决方法