Решение системы дифференциальных уравнений первого и второго порядка в Python - PullRequest
Новая
       107

Решение системы дифференциальных уравнений первого и второго порядка в Python

0 голосов
/ 30 марта 2020

Мне нужно решить следующую систему дифференциальных уравнений: 

$\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$

Начальные условия при t = 0, x2 = 1. И в момент времени t = 1, соединение y вводится в y2 компартмент, y2 = 10. Значение KR равно 1e-3.

Я решил гораздо более простую систему с использованием возведения в степень матрицы, и мне было интересно, возможно ли решить вышеуказанную систему с использованием аналогичного подхода .

У меня есть секционная модель системы X, упрощенная версия которой выглядит следующим образом:

this

Система дифференциала уравнения тогда:

the following

Я могу решить эту систему уравнений, используя следующий матричный подход.

Сначала я пишу скорость матрица [R]. Из [R] можно получить новую матрицу [A], сначала заменив каждый диагональный элемент [R] на отрицательную сумму суммы каждого из элементов строки, а затем транспонировав ее:

enter image description here

Я могу рассчитать сумму в каждом отсеке, выполнив следующие действия:

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

Этот метод хорошо работает для меня.

Модель выше (первоначальный вопрос) немного сложнее, чем просто система X. В этой модели реагенты в отделения 1 и 2 систем X и Y объединяются для получения продукта в системе Z.

X + Y -> Z с константой реакции KR.

enter image description here

, и соответствующая система дифференциальных уравнений будет иметь вид:

enter image description here

Я борюсь с методом, чтобы решить эту проблему Система дифференциальных уравнений (1-го и 2-го порядка) для расчета суммы в каждом отсеке в определенный момент времени t, учитывая исходное условие ns, KR и скорости передачи k1, k2, k3 и т. д. c ...

Могу ли я решить эту проблему, используя матричный метод, подобный приведенному выше, для системы дифференциальных уравнений первого порядка? Какие еще опции в Python у меня есть?

Заранее спасибо!

1 Ответ

1 голос
/ 01 апреля 2020

Ну, как указано в комментариях, ваш (более сложный) ODE нелинейный. Следовательно, матричный экспоненциальный подход больше не будет работать.

В целом, существует два основных подхода к решению ODE. Во-первых, вы можете попытаться найти решение Symboli c. В большинстве случаев вы придерживаетесь некоторого подхода, основанного на обоснованном предположении. Существует несколько типов ODE, для которых известны решения Symboli c.

Однако это не относится к подавляющему большинству ODE. Следовательно, мы обычно сталкиваемся с численным решением c, по существу, численно интегрируя ODE на основе правой части.

Результат - не явная функция, а приближение к значениям функции при определенных точка. В python вы можете использовать scipy для решения ODE таким способом. Исходя из вашей правой части (исключая любые ошибки с моей стороны), это будет выглядеть примерно так: 

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

Это дает следующие значения x1: 

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

на основе точки сетки 

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

Если вы хотите увидеть, как ведет себя функция, может быть достаточно интегратора. В противном случае, я бы порекомендовал прочитать учебники по символическим c подходам к ОДУ ...

...