I am using JiTCDDE to solve DDE's on a modified Oregantor model. The issue I am having is near the bifurcation point it will return negative values. While I know these are mathematically valid solutions, the Oregantor represents a chemical system. Hence, negative answers are not realistic to a real system. Is there a way to set up the code to return a minimum value for a variable when it is <=0. Below is the main part of the code I have so far.
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')
There is no direct way to restrict the signs of solutions in JiTCDDE. (Source: I am the author; I would know.)
What you need to do is to find out why your solutions become negative. So far I can think of three possible reasons (and I cannot tell you more without knowing your exact setup):
Numerical noise kicks the solution from a mathematically correct positive one to a mathematically correct negative one (but the transition between the solutions is not mathematically correct). In that case, you should be able to avoid this by decreasing the absolute tolerance (atol) or possibly the maximum step size (max_step). An alternative is to work in the logarithmic domain (also see this answer to a similar question), i.e., you make zero unreachable by construction (also, you weigh errors close to zero more strongly).
Your model allows for this to happen when it shouldn’t. In that case, and if you really need to avoid negative values, avoid a negative derivative of dynamical variables that are already pretty close to zero by wrapping a function around it that does exactly this – by using sigmoids to implement this logic (also see this answer of mine to another question). Resist the temptation to use step functions for this because integrators do not like this (further reading).
Your initial conditions are positive but not plausible. In that case, use different initial conditions. If there is no straightforward way to tell what initial conditions are plausible, it may be appropriate to pick random ones until you find a good one.
Finally, I recommend that you make your mind up what you want to happen if the solution would become negative and whether this would really be better.