Есть ли ошибка logi c в моей реализации odeint? - PullRequest
Новая
       23

Есть ли ошибка logi c в моей реализации odeint?

0 голосов
/ 30 апреля 2020

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

Я пытаюсь смоделировать простую химическую реакцию первого порядка, решая систему ODE. С точки зрения логики c мои функции верны, и я могу решить эту проблему в MATLAB без особых проблем. Я хотел бы также быть в состоянии сделать эту работу в python. Я думаю, что odeint - это инструмент для работы, но я могу ошибаться. Мое решение не должно сходиться при независимой переменной = 10 каждый раз, но оно всегда происходит независимо от входных данных. 

from matplotlib.pyplot import (plot,grid,xlabel,ylabel,show,legend)
import numpy as np 
from scipy.integrate import odeint

wght= np.linspace(0,20)
# reaction is A -> B
def PBR(fun,W):

    X,y = fun

    P_0=20;#%bar
    v_0=5; #%m^3/min
    y_A0=1; #unitless
    k=.005; #m^3/kg/min
    alpha=0.1; #1/kg
    epi=.13; #unitless
    R=8.314; #J/mol/K

    F_A0= .5 ;#mol/min

    ra = -k *y*(1-X)/(1+epi*X)


    dX = (-ra)/F_A0
    dy = -alpha*(1+epi*X)/(2*y)

    return [dX,dy]

X0 = 0.0
y0 = 1.0


sol = odeint(PBR, [X0, y0],wght)


plot(wght, sol[:, 0], 'b', label='X')
plot(wght, sol[:, 1], 'g', label='y')
legend(loc='best')
xlabel('W')
grid()
show()
print(sol)

График вывода Output graph

1 Ответ

0 голосов
/ 30 апреля 2020

Ваш график c не воспроизводится, в python 3.7, scipy 1.4.1, я получаю расхождение практически до бесконечности около 12,5, и после очистки рабочего пространства оба графика перемещаются в ноль после времени 10.

Проблема заключается в вашем делении на 2*y во втором уравнении. Фактически это означает, что y является квадратом root какой-либо другой функции, и этот тип функции ведет себя плохо при небольших значениях, когда тангенс становится вертикальным, а вскоре после этого изменяется на ноль.

Однако можно десингуляризовать деление простым способом 

    dy = -alpha*(1+epi*X)*y/(1e-10+2*y*y)

или более сложным способом, который оставляет решение далеко от нуля действительно неизменным 

    dy = -alpha*(1+epi*X)*y/(y*y+max(1e-12,y*y))

в результате получается график

enter image description here

Если это все еще не ожидаемое поведение, то в вашей системе ODE что-то другое отличается от версии Matlab.

То, что особая точка находится на отметке w=10, почти полностью зависит от alpha=0.1. Поскольку epi мало, а X изначально мало, второе уравнение близко к 

d(y^2)/dW = -alpha  ==>  y^2 = 1 - alpha*W

, и оно достигает нуля примерно на W=10, где решение должно заканчиваться как квадрат y не может принимать отрицательные значения.

...