Как я могу решить систему дифференциальных уравнений с константами в 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) являются промежуточными продуктами,

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