用sympy求解耦合的微分方程

Question

我正在尝试解决以下一阶耦合微分方程的系统：

 - dR/dt=k2*Y(t)-k1*R(t)*L(t)+k4*X(t)-k3*R(t)*I(t)
 - dL/dt=k2*Y(t)-k1*R(t)*L(t)
 - dI/dt=k4*X(t)-k3*R(t)*I(t)
 - dX/dt=k1*R(t)*L(t)-k2*Y(t)
 - dY/dt=k3*R(t)*I(t)-k4*X(t)   

已知的初始条件是：x（0）= 0，y（0）= 0
已知的常数值为：k1 = 6500，k2 = 0.9

此方程式定义了kinetick模型，我需要求解它们以获取y（t）功能以适合我的数据并找到K3和K4值。为此，我试图用SYSPY在Symbologologologologologologologology上解决。有我的代码：

import matplotlib.pyplot as plt
import numpy as np
import sympy 
from sympy.solvers.ode.systems import dsolve_system
from scipy.integrate import solve_ivp
from scipy.integrate import odeint
    
k1 = sympy.Symbol('k1', real=True, positive=True)
k2 = sympy.Symbol('k2', real=True, positive=True)
k3 = sympy.Symbol('k3', real=True, positive=True)
k4 = sympy.Symbol('k4', real=True, positive=True)
t = sympy.Symbol('t',real=True, positive=True)
L = sympy.Function('L')
R = sympy.Function('R')
I = sympy.Function('I')
Y = sympy.Function('Y')
X = sympy.Function('X')

f1=k2*Y(t)-k1*R(t)*L(t)+k4*X(t)-k3*R(t)*I(t)
f2=k2*Y(t)-k1*R(t)*L(t)
f3=k4*X(t)-k3*R(t)*I(t)
f4=-f2
f5=-f3

eq1=sympy.Eq(sympy.Derivative(R(t),t),f1)
eq2=sympy.Eq(sympy.Derivative(L(t),t),f2)
eq3=sympy.Eq(sympy.Derivative(I(t),t),f3)
eq4=sympy.Eq(sympy.Derivative(Y(t),t),f4)
eq5=sympy.Eq(sympy.Derivative(X(t),t),f5)

Sys=(eq1,eq2,eq3,eq4,eq5])
solsys=dsolve_system(eqs=Sys,funcs=[X(t),Y(t),R(t),L(t),I(t)], t=t, ics={Y(0):0, X(0):0})

有答案：

NotImplementedError: 
The system of ODEs passed cannot be solved by dsolve_system.

我也尝试过DSOLVE，但我得到了相同的。
我是否可以使用其他求解器或某种方法可以使我获得配件的功能？我在Windows64中使用Anaconda在蜘蛛中使用Python 3.8。

谢谢你！

＃update
下列的

您说的是“实验”。因此，您有数据并希望将模型适合它们，至少找到K3和K4的适当值，也许对于所有系数和初始条件（第一个测量的数据点可能不是最佳拟合的初始条件）？有关此类任务的最新尝试，请参见stackoverflow.com/questions/71722061/4。 - 卢兹·莱曼 23小时

有我的新代码：

t=[0,0.25,0.5,0.75,1.5,2.27,3.05,3.82,4.6,5.37,6.15,6.92,7.7,8.47,13.42,18.42,23.42,28.42,33.42,38.42,43.42,48.42,53.42,58.42,63.42,68.42,83.4,98.4,113.4,128.4,143.4,158.4,173.4,188.4,203.4,218.4,233.4,248.4]
yexp=[832.49,1028.01,1098.12,1190.08,1188.97,1377.09,1407.47,1529.35,1431.72,1556.66,1634.59,1679.09,1692.05,1681.89,1621.88,1716.77,1717.91,1686.7,1753.5,1722.98,1630.14,1724.16,1670.45,1677.16,1614.98,1671.16,1654.03,1661.84,1675.31,1626.76,1638.29,1614.41,1594.31,1584.73,1599.22,1587.85,1567.74,1602.69]
def derivative(S, t, k3, k4):
    k1=1798931
    k2=0.2629
    x, y,r,l,i = S
    doty = k1*r*l+k2*y
    dotx = k3*r*i-k4*x
    dotr = k2*y-k1*r*l+k4*x-k3*r*i
    dotl = k2*y-k1*r*l
    doti = k4*x-k3*r*i
    return np.array([doty, dotx, dotr, dotl, doti])
def solver(XY,t,para): 
    return odeint(derivative, XY, t, args = para, atol=1e-8, rtol=1e-11)
def integration(XY_arr,*para):
    XY0 = para[:5]
    para = para[5:]
    T = np.arange(len(XY_arr))
    res0 = solver(XY0,T, para)
    res1 = [ solver(XY0,[t,t+1],para)[-1] 
             for t,XY in enumerate(XY_arr[:-1]) ]
    
    return np.concatenate([res0,res1]).flatten()
XData =yexp
YData = np.concatenate([ yexp,yexp,yexp,yexp,yexp,yexp[1:],yexp[1:],yexp[1:],yexp[1:],yexp[1:]]).flatten()
p0 =[0,0,100,10,10,1e8,0.01]

params, info = curve_fit(integration,XData,YData,p0=p0, maxfev=5000)
XY0, para = params[:5], params[5:]
print(XY0,tuple(para))

t_plot = np.linspace(0,len(t),500)
x_plot = solver(XY0, t_plot, tuple(para))  

但是输出不正确，与初始条件P0相同：

[  0.   0. 100.  10.  10.] (100000000.0, 0.01)

图

我知道在每个时间瞬间，功能“集成”给出了每个函数的Y值，但是我不知道如何解开它们以将curve_fitt分别制成。也许我不太了解它是如何工作的。

谢谢你！

Answer 1

正如您观察到的那样，Sympy无法解决该系统。这可能意味着将

• Sympy中的ODE分类的过程不够完整，或者
• 需要某些技巧/方法，而不是Sympy实现的标准方法集，或者
• 只是没有符号解决方案。

最后一个情况是通用的情况，采用象征性解决的颂歌，添加一些随机术语，几乎可以肯定，所产生的颂歌不再象征性地解决。

---

据我了解，您可以通过状态空间（CX，CY，CR，CL，CI）带有4个参数K1，K1，K2，K2，K3，K4的ode系统具有模型。，并且通过反应系统的结构r+i＆lt; - ＆gt; x，r+l＆lt; - ＆gt; y，总和cr+cx+cy，cl+cy，ci+cx都是恒定的。

对于大约由模型表示的其他过程，您有时间序列数据t [k]，y [k] y component。另外，您还拥有有关初始状态和参数集的部分信息。如果有足够多的数据点，人们也可以忘记这些数据，适合所有参数，并将给定参数与计算的参数相比。

有几个模块和软件包以或多或少的抽象方式解决了这项合适的任务。我认为可以使用Pyomo和Gekko。更直接的人可以使用scipy.odr或scipy.ptimize的设施。

正向函数

def model(t,u,k1,k2,k3,k4):
    X,Y,R,L,I = u
    dL = k2*Y - k1*R*L
    dI = k4*X - k3*R*I
    dR = dL+dI
    dX = -dI
    dY = -dL
    return dX,dY,dR,dL,dI

def solver(t,u0,k):
    res = solve_ivp(model, [0, t[-1]], u0, args=tuple(k), t_eval=t, 
                           method="DOP853", atol=1e-7, rtol=1e-10)
    return res.y

定义转换时间和参数准备一些数据的

k1o = 6.500; k2o=0.9
T = np.linspace(0,0.05,21)
U = solver(T, [0,0,50,40,25], [k1o, k2o, 5.400, 0.7])
Y = U[1] # equilibrium slightly above 30
Y += np.random.uniform(high=0.05, size=Y.shape)

，并准备以初始状态和系数分配组合参数向量的功能，调用曲线拟合功能，

from scipy.optimize import curve_fit

def partial(t,R,L,I,k3,k4): 
    print(R,L,I,k3,k4)
    U = solver(t,[0,0,R,L,I],[k1o,k2o,k3,k4])
    return U[1]

params, info = curve_fit(partial,T,Y, p0=[30,20,10, 0.3,3.000])
R,L,I, k3,k4 = params
print(R,L,I, k3,k4)

结果证明curve_fit进入奇怪的区域具有较大的负值。可能的原因是，y组件最终与所有其他组件相结合不足，这意味着某些参数的很大变化对y ，因此y中的最小噪声可以导致这些参数中的较大偏差。这显然发生（首先）到k3。