Не удается импортировать проблему XЖесткий ODE решатель для модели Oregonator - PullRequest
Новая
       107

Не удается импортировать проблему XЖесткий ODE решатель для модели Oregonator

0 голосов
/ 18 октября 2019

Ошибка возникает при попытке импортировать метод Радау из scipy.integrate (необходим, потому что модель Oregonator является жесткой системой).

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

Простите меня за неопытность, я новичок в Python.

Ошибка: ImportError: не могу импортировать имя 'radau' из 'scipy.integrate'

По моему незнанию, япостроил метод Рунге-Кутта 4-го порядка с нуля. Найдя цены на акции вместо химических колебаний, я перешел к использованию odeint. Это все еще не удалось. Только после этого я открыл идею жесткой системы, поэтому я работал над методом Радау. 

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate. import radau

# Dimensionless parameters
e = 0.04
q = 0.0008
f = 1.0

# Oregonator model
def Oregonator(Y, t):
    return [((Y[0] * (1 - Y[0])  - ((Y[0] - q) * f * Y[1]) // (q + Y[0]))) 
    // e, Y[0] - Y[1]]

# Time span and inital conditions
ts = np.linspace(0, 10, 100)
Y0 = [1, 3]

# Numerical algorithm/method
NumSol = radau(Oregonator, 0, Y0, t_bound=30)
x = NumSol[:,0]
z = NumSol[:,1]

Ожидаемые результаты должны быть такими же, как в (стр. 12): https://pdfs.semanticscholar.org/0959/9106a563e9d88ce6442e3bb5b242d5ccbdad.pdf только для x и z. Отсутствие y связано с использованным мною стационарным приближением.

1 Ответ

0 голосов
/ 19 октября 2019

Используйте solve_ivp в качестве однострочного интерфейса для классов решателей, таких как RK45 или Radau. Используйте правильную заглавную букву класса. Используйте правильный порядок аргументов в функции ODE (вы можете использовать tfirst=True в odeint, чтобы использовать одну и ту же функцию во всем). Избегайте целочисленного деления там, где вы намереваетесь использовать деление с плавающей запятой. 

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

# Dimensionless parameters
eps = 4e-2
q = 8e-4
f = 2.0/3

# Oregonator model
def Oregonator(t,Y):
    x,z = Y;
    return [(x * (1 - x) + (f*(q-x)*z) / (q + x)) / eps, x - z]

# Time span and inital conditions
ts = np.linspace(0, 10, 100)
Y0 = [1, 0.5]

# Numerical algorithm/method
NumSol = solve_ivp(Oregonator, [0, 30], Y0, method="Radau")
x, z = NumSol.y
y = (f*z)/(q+x)
t = NumSol.t
plt.subplot(221);
plt.plot(t,x,'b'); plt.xlabel("t"); plt.ylabel("x");
plt.subplot(222);
plt.plot(t,y,'r'); plt.xlabel("t"); plt.ylabel("y");
plt.subplot(223);
plt.plot(t,z,'g'); plt.xlabel("t"); plt.ylabel("z");
plt.subplot(224);
plt.plot(x,z,'k'); plt.xlabel("x"); plt.ylabel("z");
plt.tight_layout(); plt.show()

Затем получается график

enter image description here

с периодическими колебаниями.

Дальнейшие действия могут заключаться в использовании опции tspan или «плотного вывода» для получения выборок решения в определенных пользователем точках выборки. Для получения достоверных результатов установите допуски ошибок вручную.

f=0.51262 находится близко к точке перехода от сходящегося к колебательному поведению. enter image description here

...