Проблемы с функцией и odeint в Python
В течение нескольких месяцев я начал работать с Python, учитывая его большие преимущества. Но недавно я использовал odeint от scipy для решения системы дифференциальных уравнений. Но в процессе интеграции реализованная функция работает не так, как ожидалось.
В этом случае я хочу решить систему дифференциальных уравнений, в которой одно из начальных условий (x[0]) изменяется (между 4-5) в зависимости от значения, которое переменная достигает в процессе интеграции (она программируется внутри функция с помощью структуры if).
#Control of oxygen
SO2_lower=4
SO2_upper=5
if x[0]<=SO2_lower:
x[0]=SO2_upper
Когда функция используется odeint, некоторые строки кода внутри функции исключаются, даже когда функции изменяют значение x[0]. Вот весь мой код:
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
plt.ion()
# Stoichiometric parameters
YSB_OHO_Ox=0.67 #Yield for XOHO growth per SB (Aerobic)
YSB_Stor_Ox=0.85 #Yield for XOHO,Stor formation per SB (Aerobic)
YStor_OHO_Ox=0.63 #Yield for XOHO growth per XOHO,Stor (Aerobic)
fXU_Bio_lys=0.2 #Fraction of XU generated in biomass decay
iN_XU=0.02 #N content of XU
iN_XBio=0.07 #N content of XBio
iN_SB=0.03 #N content of SB
fSTO=0.67 #Stored fraction of SB
#Kinetic parameters
qSB_Stor=5 #Rate constant for XOHO,Stor storage of SB
uOHO_Max=2 #Maximum growth rate of XOHO
KSB_OHO=2 #Half-saturation coefficient for SB
KStor_OHO=1 #Half-saturation coefficient for XOHO,Stor/XOHO
mOHO_Ox=0.2 #Endogenous respiration rate of XOHO (Aerobic)
mStor_Ox=0.2 #Endogenous respiration rate of XOHO,Stor (Aerobic)
KO2_OHO=0.2 #Half-saturation coefficient for SO2
KNHx_OHO=0.01 #Half-saturation coefficient for SNHx
#Other parameters
DT=1/86400.0
def f(x,t):
#Control of oxygen
SO2_lower=4
SO2_upper=5
if x[0]<=SO2_lower:
x[0]=SO2_upper
M=np.matrix([[-(1.0-YSB_Stor_Ox),-1,iN_SB,0,0,YSB_Stor_Ox],
[-(1.0-YSB_OHO_Ox)/YSB_OHO_Ox,-1/YSB_OHO_Ox,iN_SB/YSB_OHO_Ox-iN_XBio,0,1,0],
[-(1.0-YStor_OHO_Ox)/YStor_OHO_Ox,0,-iN_XBio,0,1,-1/YStor_OHO_Ox],
[-(1.0-fXU_Bio_lys),0,iN_XBio-fXU_Bio_lys*iN_XU,fXU_Bio_lys,-1,0],
[-1,0,0,0,0,-1]])
R=np.matrix([[DT*fSTO*qSB_Stor*(x[0]/(KO2_OHO+x[0]))*(x[1]/(KSB_OHO+x[1]))*x[4]],
[DT*(1-fSTO)*uOHO_Max*(x[0]/(KO2_OHO+x[0]))*(x[1]/(KSB_OHO+x[1]))* (x[2]/(KNHx_OHO+x[2]))*x[4]],
[DT*uOHO_Max*(x[0]/(KO2_OHO+x[0]))*(x[2]/(KNHx_OHO+x[2]))*((x[5]/x[4])/(KStor_OHO+(x[5]/x[4])))*(KSB_OHO/(KSB_OHO+x[1]))*x[4]],
[DT*mOHO_Ox*(x[0]/(KO2_OHO+x[0]))*x[4]],
[DT*mStor_Ox*(x[0]/(KO2_OHO+x[0]))*x[5]]])
Mt=M.transpose()
MxRm=Mt*R
MxR=MxRm.tolist()
return ([MxR[0][0],
MxR[1][0],
MxR[2][0],
MxR[3][0],
MxR[4][0],
MxR[5][0]])
#ODE solution
t=np.linspace(0.0,3600,3600)
#Initial conditions
y0=np.array([5,176,5,30,100,5])
Var=odeint(f,y0,t,args=(),h0=1,hmin=1,hmax=1,atol=1e-5,rtol=1e-5)
Sol=Var.tolist()
plt.plot(t,Var[:,0])
Большое спасибо заранее!!!!!
1 ответ
Короткий ответ:
Вы не должны изменять вектор входного состояния внутри вашей функции ODE. Вместо этого попробуйте следующее и проверьте свои результаты:
x0 = x[0]
if x0<=SO2_lower:
x0=SO2_upper
# use x0 instead of x[0] in the rest of this function body
Я полагаю, что это ваша проблема, но я не уверен, так как вы не объяснили, что именно не так с результатами. Более того, вы не меняете "начальные условия". Начальное условие
y0=np.array([5,176,5,30,100,5])
Вы просто меняете вектор входного состояния.
Подробный ответ:
Ваш интегратор odeint, вероятно, использует один из адаптивных методов Рунге-Кутты высшего порядка. Этот алгоритм требует нескольких оценок функции ODE для расчета одного шага интегрирования, поэтому изменение вектора входного состояния может привести к неопределенным результатам. В C++ boost::odeint это даже невозможно сделать, потому что входной переменной является "const". Python, однако, не так строг, как C++, и я полагаю, что такую ошибку можно сделать непреднамеренно (хотя я и не пробовал).
РЕДАКТИРОВАТЬ:
Хорошо, я понимаю, чего ты хочешь достичь.
Ваша переменная x[0] ограничена модульной алгеброй, и ее невозможно выразить в виде
x' = f(x,t)
которое является одним из возможных определений обыкновенного дифференциального уравнения, которое должна решить библиотека ondeint. Однако здесь можно использовать несколько возможных "хаков", чтобы обойти это ограничение.
Одна из возможностей - использовать фиксированный шаг и младший порядок (потому что для решателей более высокого порядка вам нужно знать, в какой части алгоритма вы находитесь, см., Например, RK4) решатель и изменить уравнение dx[0] (в вашем коде это элемент MxR[0][0]) для:
# at the beginning of your system
if (x[0] > S02_lower): # everything is normal here
x0 = x[0]
dx0 = # normal equation for dx0
else: # x[0] is too low, we must somehow force it to become S02_upper again
dx0 = (x[0] - S02_upper)/step_size // assuming that x0_{n+1} = x0_{n} + dx0*step_size
x0 = S02_upper
# remember to use x0 in the rest of your code and also remember to return dx0
Однако я не рекомендую этот метод, потому что он сильно зависит от алгоритма, и вы должны знать точный размер шага (хотя я могу порекомендовать его для ограничений насыщения). Другая возможность - выполнять один шаг интеграции за раз и корректировать x0 каждый раз, когда это необходимо:
// 1 do_step(sys, in, dxdtin, t, out, dt);
// 2 do something with output
// 3 in = out
// 4 return to 1 or finish
Извините за синтаксис C++, здесь приведена исчерпывающая документация ( степперы C++ odeint), а также ее эквивалент в языке Python ( класс одонов Python). Интерфейс C++ odeint лучше для вашей задачи, однако вы можете добиться того же в Python. Просто посмотрите на:
integrate(t[, step, relax])
set_initial_value(y[, t])
в документах.