我正在使用 JiTCDDE 在修改后的 Oregantor 模型上解决 DDE。我遇到的问题在分岔点附近,它将返回负值。虽然我知道这些是数学上有效的解决方案,但 Oregantor 代表了一个化学系统。因此,否定答案对于真实系统来说是不现实的。有没有办法设置代码以在 <= 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')