用条件语句求解一组 ODE

Question

我正在尝试用Python求解一组微分方程。基本上，我想绘制下图：

该图显示了营养物质（Am）和浮游植物（Phy）的浓度如何随时间变化。 我必须将两个变量绘制为时间函数的 Python 代码：

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

# parameters
dil = 0.2            
har_f = 20           
har_pc = 0.95       
ext_Am = 100
umax_Phy = 0.693
kAm_Phy = 14
t0 = 0           
tf = 100             
ns = 200
dt = (tf-t0)/(ns-1)

params = [dil, har_f, har_pc, ext_Am, umax_Phy, kAm_Phy, t0, tf, ns, dt]

def odes(t,x):
    dx = [0,0]
    if t == 0:
        har_dil = 0
    elif t > 0 and t % har_f == 0:
        har_dil = har_pc/dt
    else:
        har_dil = 0
        
    har_dil = 0
    time_dil = params[0] + har_dil
    in_Am = params[3] * time_dil
    u_Phy = params[4]*x[0]/(x[0]+params[5])
    gro_Phy = u_Phy * x[1]
    out_Am = x[0] * time_dil
    out_Phy = x[1] * time_dil
    dx[0] = in_Am - out_Am - gro_Phy
    dx[1] = gro_Phy - out_Phy

    return dx

# initial conditions:
x0 = [99,1]

# declare a time vector(time window)
tspan = (0,100)

# solving the set of ODEs
sol = solve_ivp(odes,tspan,x0)

fig1,axs_lst = plt.subplots(nrows=1,ncols=1,squeeze=False)

axs_lst[0,0].plot(sol.t,sol.y[0, :],'-',color='red',label='Am')
axs_lst[0,0].plot(sol.t,sol.y[1, :],color='green',label='Phy')
axs_lst[0,0].legend(loc='best',frameon=True,edgecolor='black',fancybox=False)
axs_lst[0,0].set_xlabel('Time [day]',size='12')
axs_lst[0,0].set_ylabel('Concentration [$g\:biomass\:m^{-3}$]',size='12')
axs_lst[0,0].tick_params(axis='both', which='both', direction='in')
axs_lst[0,0].grid(True)

但是生成的图表如下所示：

正如你所看到的，营养物质的浓度并没有按其应有的那样减少，因为变量har_dilaccount的值只是条件语句的最终值。如果有人可以帮助我制定解决此问题的策略，我将不胜感激。

Answer 1

你写道：

    if t == 0:
        har_dil = 0
    elif t > 0 and t % har_f == 0:
        har_dil = har_pc/dt
    else:
        har_dil = 0
        
    har_dil = 0

好吧，那里发生了一些事情。

(1.) 你在计算

har_dil

时遇到了一些麻烦， 然后你无条件地将其归零。 这不可能是你想要的。 建议您将

if

作为辅助函数， 并无条件分配该函数的结果。

(2.)

t % har_f == 0

表达式显然是 试图注入一些周期性，如你的第一张图所示。 好的。 但是，公平警告，询问是否 FP 数量恰好等于零是相当偶然的。 建议您将

int(t)

放入该表达式中， 或者询问模结果是否“接近”零。 也就是说，选择一个小的 epsilon 并询问是否

abs(t % har_f) < epsilon

。