如何在python中求解带常数的微分方程组？

Question

我想在不定义常量的情况下，使用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的情况下求解这个系统

注释：我正在描述（3）、（4）、（5）为中间产物的DNA链置换反应的动力学方程

（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和常数之间的函数。你知道吗

Answer 1

定义一个向量函数来包含动态

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()