使用 JiTCDDE 将解决方案限制为正值

Question

我正在使用 JiTCDDE 对修改后的 Oregantor 模型求解 DDE。我遇到的问题是在分叉点附近它将 return 负值。虽然我知道这些是数学上有效的解决方案，但 Oregantor 代表了一个化学系统。因此，否定答案对于真实系统来说是不现实的。有没有一种方法可以将代码设置为 return 变量的最小值 <=0。以下是我到目前为止的代码的主要部分。

def P1(k):
        return(
            ((H*y(k))/(k01+H*y(k)+kl*H*H*A))*phi
        )

def C(i,j):

    return(
        M1 * ( y(j,t-tau1)-y(i) ) 
      + M2 * ( y(j,t-tau2)-y(i) ) 
    )




MO4 = [
 k1*A*y(1)-k2*y(0)*y(1)+ k3*A*y(0)-2.0*k4*y(0)*y(0)-(y(0)-xsur)*kf,           #HBrO2
 -k1*A*y(1)-k2*y(0)*y(1)+f1*k5*y(2)-(y(1)-ysur)*kf+P1(3)+C(2,6),              #Bromide
 2*k3*A*y(0)-k5*y(2)+P1(3)+C(2,6),                                            #Cataylst
 k1*A*y(1)+2*k2*y(0)*y(1)+k4*y(0)*y(0)-k6*y(3)-(y(3)-vsur)*kf-P1(3)-C(2,6),   #BrMa
 k1*A*y(5)-k2*y(4)*y(5)+ k3*A*y(4)-2.0*k4*y(4)*y(4)-(y(4)-xsur)*kf,           #HBrO2
 -k1*A*y(5)-k2*y(4)*y(5)+f2*k5*y(6)-(y(5)-ysur)*kf+P1(7)+C(6,2),              #Bromide
 2*k3*A*y(4)-k5*y(6)+P1(7)+C(6,2),                                            #Cataylst
 k1*A*y(5)+2*k2*y(4)*y(5)+k4*y(4)*y(4)-k6*y(7)-(y(7)-vsur)*kf-P1(7)-C(6,2),   #BrMa
]

I = jitcdde(MO4)
I.set_integration_parameters(rtol=1e-7,atol=1e-7)
I.constant_past ([0,1.0e-6,0,0,1.0e-6,1.0e-6,1.0e-6,1.0e-6], time=0.0)
I.step_on_discontinuities(max_step=.00001)


data=[]
for time in times:
    data.append( I.integrate(time))
np.savetxt('peaks_%d.dat'%(i), data,)

data1=np.loadtxt('peaks_%d.dat'%(i),dtype = float,delimiter=' ',skiprows=200,usecols=(2,6)).T #,skiprows=80
plt.plot(data1[0],'r')
plt.plot(data1[1],'-.b')
plt.title( 'Catalyst ' )
plt.xlabel('time(sec)')
plt.ylabel('Amplitude')
plt.show()
print('DONE')

Answer 1

JiTCDDE 中没有直接的方法来限制解的符号。（来源：我是作者；我会知道的。）

你需要做的是找出为什么你的解会变成负数。到目前为止，我可以想到三个可能的原因（如果不知道您的确切设置，我无法告诉您更多）：

数值噪声将解从数学上正确的正解踢到数学上正确的负解（但解之间的转换在数学上不正确）。在这种情况下，您应该能够通过减小绝对容差 (atol) 或可能的最大步长 (max_step) 来避免这种情况。另一种方法是在对数域中工作（另请参阅），即通过构造使零无法达到（同时，您更强烈地权衡接近零的误差）。
您的模型允许这种情况在不应该发生的时候发生。在那种情况下，如果你真的需要避免负值，通过在它周围包装一个函数来避免已经非常接近于零的动态变量的负导数——通过使用 sigmoids 来实现这个逻辑（另见). Resist the temptation to use step functions for this because integrators do not like this (further reading).
您的初始条件是肯定的，但不合理。在这种情况下，使用不同的初始条件。如果没有直截了当的方法来判断哪些初始条件是合理的，那么在找到合适的初始条件之前，随机选择可能是合适的。

最后，我建议您下定决心，如果解决方案变为负面，您希望发生什么，以及这是否真的会更好。

使用 JiTCDDE 将解决方案限制为正值

Limit solutions to a positive value using JiTCDDE

python

differential-equations

jitcode-jitcdde-jitcsde