Imposing a constraint on a system of ODEs to be solved in Mathematica

I have a system of three ODEs which I am trying to solve, let’s call them x(t), y(t) and z(t). Two constrains must be placed on the solutions: x(t)+y(t)+z(t)=1 & {x(t),y(t),z(t)} ∈ [0,1]. This has proved to be a difficult ta
See here the original piece of code which was used to attempt to solve the aforementioned math problem:

1st code block. StackOverflow wouldn’t accept it as a code block, so I had to add an image. Apologies.

I was working with this for the longest time, trying to figure out how to implement the “x(t)+y(t)+z(t)=1” constraint without ending up with an under-determined system (4 equations, yet only 3 unknowns), until it occurred to me that I could just substitute one variable for the others. This produced a new piece of code:

2nd code block

The good news is that all 3 variables do indeed add up to 1. The bad news is that some variables become negative in order to compensate for other variables becoming much greater than 1. Hence the need for the {x(t),y(t),z(t)} ∈ [0,1], which is ultimately what I am hung up on right now. This is the best I could come up with on my own:

3rd code block

This still produces values outside of the [0,1] range, which I do not desire. Any assistance would be appreciated.

  • Text for the first code block eqn1=Θh'[t]==-2 km1 Θh[t]^2+km1(Θs[t]+Θ2h[t])+(k3+k4)pch4 Θs[t]-(km3+km4)pch3 Θh[t]; bc1=Θh[0]==1; eqn2=Θ2h'[t]==2 k1 Θh[t]^2-km1(Θs[t]+Θ2h[t])-k2 Θ2h[t]+km2 ph2 Θs[t]; bc2=Θ2h[0]==0; eqn3=Θs'[t]==2 k1 Θh[t]^2-km1(Θs[t]+Θ2h[t])+k2 Θ2h[t]-km2 ph2 Θs[t]-(k3+k4)pch4 Θs[t]+(km3+km4)pch3 Θh[t]; bc3=Θs[0]==0; eqns={eqn1,eqn2,eqn3,bc1,bc2,bc3}; sol=NDSolve[eqns,{Θh[t],Θ2h[t],Θs[t]},{t,0,10^-9}] Without any of the values for your constants it is difficult to experiment with this to see if ideas work or not. Clip can enforce ranges but NDSolve probably won’t like it


Leave a Comment