Search code examples
pythonmatrixdifferential-equationsodeintexponentiation

Solving a system of first and second order differential equations in Python

I need to solve the following system of differential equations:

$\frac{dx_1}{dt} = -k_1x_1+k_2x_2-(K_R)x_1y_1$
$\frac{dx_2}{dt} =  k_1x_1-k_2x_2-k_3x_2-(K_R)x_2y_2$
$\frac{dx_3}{dt} =  k_3x_3$
$\frac{dy_1}{dt} = -k_1y_1+k_2y_2-(K_R)x_1y_1$
$\frac{dy_2}{dt} =  k_1y_1-k_2y_2-k_3y_2-(K_R)x_2y_2$
$\frac{dy_3}{dt} =  k_3y_3$
$\frac{dz_1}{dt} = -k_1z_1+k_2z_2+(K_R)x_1y_1$
$\frac{dz_2}{dt} =  k_1z_1-k_2z_2-k_3z_2+(K_R)x_2y_2$
$\frac{dz_3}{dt} =  k_3z_3$

The initial conditions at t = 0, is x2 = 1. And at time t = 1, a compound y is introduced in y2 compartment, y2 = 10. The value of KR is 1e-3.

I have solved a much simpler system using exponentiation of matrix, and was wondering whether it's possible to solve the above system using similar approach.

I have a compartmental model system X, a simplified version of which, looks like this:

this

The system of differential equations is then:

the following

I can solve this system of equations using the following matrix approach.

First, I write the rate matrix [R]. From [R] one can obtain a new matrix [A] by first replacing each diagonal element of [R] by the negative of the sum of each of row elements, and then transposing it:

enter image description here

I can calculate the amount in each compartment by doing the following:

enter image description here

In python:

RMatrix = model_matrix.as_matrix()
row, col = np.diag_indices_from(RMatrix)
RMatrix[row, col] = -(RMatrix.sum(axis=1)-RMatrix[row,col])
AMatrix = RMatrix.T

def content(t):
    cont = np.dot(linalg.expm(t*AMatrix), x0))

This method is working well for me.

The model above (the original question) is a little more complicated than just System X. In this model, reactants in compartments 1 and 2 of Systems X and Y combine to get product in System Z.

X + Y --> Z, with a reaction constant of KR.

enter image description here

, and the corresponding system of differential equations would be:

enter image description here

I am struggling with a method to solve this system of differential equations (1st and 2nd order) to calculate the amount in each compartment at a certain time t, given the initial conditions, KR, and the transfer rates k1, k2, k3, etc...

Can I solve it using the matrix method like the one above for a system of first order differential equations? What other options in Python do I have?

Thanks in advance!

Solution

  • Well, as pointed out in the comments, your (more complicated) ODE is nonlinear. Therefore, a matrix exponential approach will not work anymore.

    In general, there are two general approaches to solving ODEs. First, you can try to find a symbolic solution. In most cases, you follow some approach based on an educated guess. There are several types of ODEs for which symbolic solutions are known.

    However, this is not the case for the vast majority of ODEs. Hence, we generally contend ourselves with a numeric solution, essentially numerically integrating the ODE based on the right hand side.

    The result is not an explicit function, but instead an approximation on the function values at certain point. In python, you can use scipy to solve ODEs this way. Based on your right hand side (barring any errors on my part), this would look something like this:

    import numpy as np

import scipy.integrate 

k_1 = 1
k_2 = 1
k_3 = 1
K_R = 1

def eval_f(v, t):
    [x, y, z] = np.split(v, [3, 6])

    return np.array([-k_1*x[0] +k_2*x[1] - (K_R)*x[0]*y[0],
                     k_1*x[0] - k_2*x[1] - k_3*x[1] - (K_R)*x[1]*y[1],
                     k_3*x[2],
                     - k_1*y[0] + k_2*y[1] - (K_R)*x[0]*y[0],
                     k_1*y[0] - k_2*y[1] - k_3*y[1] - (K_R)*x[1]*y[1],
                     k_3*y[2],
                     - k_1*z[0] + k_2*z[1] + (K_R)*x[0]*y[0],
                     k_1*z[0] - k_2*z[1] - k_3*z[1] + (K_R)*x[1]*y[1],
                     k_3*z[2]])

initial = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9])

t = np.linspace(0, 1, 10)

values = scipy.integrate.odeint(eval_f, initial, t)

# variable x[0]
print(values[:,0])

    This yields the following values of x1:

    [1.         0.70643591 0.49587121 0.35045691 0.25034256 0.1809533
 0.13237994 0.09800056 0.07338967 0.05557138]

    based on the grid points

    [0.         0.11111111 0.22222222 0.33333333 0.44444444 0.55555556
 0.66666667 0.77777778 0.88888889 1.        ]

    If you want to see how the function behaves, an integrator may be sufficient. Otherwise, I would recommend reading up on symbolic approaches to ODEs in a text book...