python - 如何在 python 中求解具有常数的微分方程组?
问题描述
我想在不定义常量的情况下使用 Python 解决以下系统。
dx1(t)/dt = - kf1*x1(t)*x2(t) + kr1*x3(t)
dx2(t)/dt = - kf1*x1(t)*x2(t) + kr1*x3(t)
dx3(t)/dt = kf1*x1(t)*x2(t) - kr1*x3(t) - k2*(x3(t) - x4(t))
dx4(t)/dt = k2*(x3(t) - x4(t)) + kf3*x5(t)*x6(t) - kr3*x4(t)
dx5(t)/dt = -kf3*x5(t)*x6(t) + kr3*x4(t)
dx6(t)/dt = -kf3*x5(t)*x6(t) + kr3*x4(t)
x1(0)=x1_0, x2(0)=x2_0 and x(3)=x(4)=x(5)=x(6)=0
我想解决系统而不用kf1,kr1,k2,kf3,kr3,x1_0 and x2_0
实数替换
评论:我正在描述 DNA 链置换反应的动力学方程,其中 (3)、(4) 和 (5) 是中间产物,
(1) + (2) <--> (5) + (6)
我尝试使用 sympy 并将我的常量定义为符号但没有成功
from sympy import *
x1,x2,x3,x4,x5,x6 =symbols('x1 x2 x3 x4 x5 x6', cls=Function)
kf1,kr1,k2,kf3,kr3 = symbols("kf1 kr1 k2 kf3 kr3")
diffeqq1=Eq(x1(t).diff(t), - kf1*x1*x2 + kr1*x3)
diffeqq2=Eq(x2(t).diff(t), - kf1*x1*x2 + kr1*x3)
diffeqq3=Eq(x3(t).diff(t), kf1*x1*x2 - kr1*x3 - k2*(x3 - x4))
diffeqq4=Eq(x4(t).diff(t), k2*(x3 - x4) + kf3*x5*x6 - kr3*x4)
diffeqq5=Eq(x5(t).diff(t), -kf3*x5*x6 + kr3*x4)
diffeqq6=Eq(x6(t).diff(t), -kf3*x5*x6 + kr3*x4)
dsolve(system,[x1,x2,x3,x4,x5,x6])
我想要的结果是 x1,x2,x5,x6 和常量之间的函数。
解决方案
您定义一个向量函数来包含动态
def reactions(t,u):
x1,x2,x3,x4,x5,x6 = u
F1 = kf1*x1*x2
R1 = kr2*x3
FR2 = k2*(x3-x4)
F3 = kf3*x5*x6
R3 = kr3*x4
return [ -F1+R1, -F1+R1, F1-R1-FR2, FR2+F3-R3, -F3+R3, -F3+R3]
然后scipy.integrate.solve_ivp
在足够长的时间跨度内调用求解系统,
res = solve_ivp(reactions, [t0, tf], x_init, dense_output=True, atol=1e-8, rtol=1e-6)
t = np.linspace(t0, tf, 702)
u = res.sol(t)
x1,x2,x3,x4,x5,x6 = u
plt.plot(t,x1,lw=3,label='$x_1$')
plt.plot(t,x2,lw=3,label='$x_2$')
plt.plot(t,x6,lw=3,label='$x_6$')
plt.grid(); plt.legend(); plt.show()
推荐阅读
- javascript - 如何做到这一点,当您按下 png 图标时,它会将您带到 NavLink 中“/”的主页,而无需重新加载页面?
- javascript - 悬停时更改底部边框的颜色
- angular5 - 纱线链接,对链接包的更改未反映在主机应用程序中
- javascript - 将主要 Typescript 代码及其依赖项组合在一个文件中
- activiti - 如何在不取消用户任务的情况下通过信号边界事件启动子进程
- python - Python:带有自定义标头的请求失败
- java - 扩展类的构造函数被调用两次
- javascript - 角度路由不显示模板
- dictionary - Thrift IDL 语法错误 - 接受并返回映射
- c# - 如何配置 swashbuckle 以显示 api 版本而不是 v{version} 变量?