Solving coupled differential equations with sympy

Question

I am trying to solve the following system of first order coupled differential equations:

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

The knowed initial conditions are: X(0)=0, Y(0)=0
The knowed constants values are: k1=6500, k2=0.9

This equations defines a kinetick model and I need to solve them to get the Y(t) function to fit my data and find k3 and k4 values. In order to that, I have tried to solve the system simbologically with sympy. There is my code:

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

There is the answer:

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

I have tried with dsolve too, but I get the same.
Is there any other solver I can use or some way of doing this that will allow me to get the function for the fitting? I'm using python 3.8 in Spider with Anaconda in windows64.

Thank you!

# Update
Following

You are saying "experiment". So you have data and want to fit the model to them, find appropriate values for k3 and k4 at least, and perhaps for all coefficients and the initial conditions (the first measured data point might not be the initial condition for the best fit)? See stackoverflow.com/questions/71722061/… for a recent attempt on such a task. – Lutz Lehmann 23 hours

There is my new code:

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

But the output are not correct, as are the same as initial condition p0:

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

Graph

I understand that the function 'integration' gives packed values of y for each function at each instant of time, but I don't know how to unpack them to make the curve_fitt separately. Maybe I don't quite understand how it works.

Thank you!

Answer 1

As you observed, sympy is not able to solve this system. This might mean that

• the procedure to classify ODE in sympy is not complete enough, or
• some trick/method is needed above the standard set of methods that is implemented in sympy, or
• that there just is no symbolic solution.

The last case is the generic one, take a symbolically solvable ODE, add some random term, and almost certainly the resulting ODE is no longer symbolically solvable.

---

As I understand with the comments, you have an model via ODE system with state space (cX,cY,cR,cL,cI) with equations with 4 parameters k1,k2,k3,k4 and, by the structure of a reaction system R+I <-> X, R+L <-> Y, the sums cR+cX+cY, cL+cY, cI+cX are all constant.

For some other process that is approximately represented by the model, you have time series data t[k],y[k] for the Y component. Also you have partial information on the initial state and the parameter set. If there are sufficiently many data points one could also forget about these, fit for all parameters, and compare how far away the given parameters are to the computed ones.

There are several modules and packages that solve this fitting task in a more or less abstract fashion. I think pyomo and gekko can both be used. More directly one can use the facilities of scipy.odr or scipy.optimize.

Define the forward function that transforms time and parameters

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

Prepare some data plus noise

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)

Prepare the function that splits the combined parameter vector in initial state and coefficients, call the curve fitting function

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)

It turns out that curve_fit goes into strange regions with large negative values. A likely reason is that the Y component is, in the end, not coupled strongly enough to all the other components, meaning that large changes in some of the parameters have minimal influence on Y, so that minimal noise in Y can lead to large deviations in these parameters. Here this apparently happens (first) to k3.