# Simple question in NDSolveValue (Want to set maximum in the condition)

Thank you for reading my question.

I'm trying to solve kinetic process using NDSolveValue.

The equations are simply

Pa'[t]=ra[t]*Pb[t]-rb[t]*Pb[t],
Pb[t]=rb[t]*Pa[t]-ra[t]*Pb[t],
Pa=0.5, Pb=0.5


However, I want Pa[t] and Pb[t] to stay within 0 to 1. But I don't know how to introduce this boundary condition to the code. Can you help me with this?

Below is the real code I'm using. ra, rb, are actually having another parameter in here but they are constant that I put for my convenience.

B[t_]:=91*Cos[0.000075*pi*t];
Ea[B_] := (24.5 + B^2/(4*24.5) - B)*(1 - HeavisideTheta[B-2*24.5])-(24.5 + B^2/(4*24.5) + 24.5)*HeavisideTheta[-B - 2*24.5]
γA[t] := 10*EXP[-Ea[B[t]]];
γB[t] := 10*EXP[-Ea[-B[t]]];

sol :=
sol =
NDSolveValue[{
Pa'[t] == (γB[t] * Pb[t] - γA[t] * Pa[t]),
Pb'[t] == (-γB[t] * Pb[t] + γA[t] * Pa[t]),
Pa == 0.5, Pb == 0.5},
{Pa, Pb}, {t, 3/0.000075}]

• Can you give explicit values for your parameters so your code can be executed? – MarcoB Jun 25 '18 at 21:55
• Thank you for your interest. They are actually having multi-variables but I changed those parts into numbers. Is it possible to make Pa[t] and Pb[t] not exceed over 1 or below 0? Thank you. – JUNSANG MOON Jun 26 '18 at 1:28
• Your code still does not execute correctly, even with the explicit constants. Please fix it so it gives a result, even though it is not yet what you want. Others may be able to build upon it then. If the code does not work, though, you ask us to do double duty, i.e. to fix the code AND to solve your problem. That's rather uninviting. – MarcoB Jun 26 '18 at 1:38
• If you add your two equations, you'll find that Pa'[t]+Pb'[t]==0, so you should automatically have Pa[t]+Pb[t]==1 for all time based on your initial conditions. Incidentally, I think there's a typo in the first equation (in the text, not in your code). – Chris K Jun 26 '18 at 2:29

B[t_] := 91*Cos[0.000075*Pi*t];
Ea[B_] := (24.5 + B^2/(4*24.5) - B)*(1 -
HeavisideTheta[B - 2*24.5]) - (24.5 + B^2/(4*24.5) + 24.5)*
HeavisideTheta[-B - 2*24.5]
\[Gamma]A[t_] := 10*Exp[-Ea[B[t]]];
\[Gamma]B[t_] := 10*Exp[-Ea[-B[t]]];
Pa1[t_] :=
If[0 < Pa[t] < 1, (\[Gamma]B[t]*Pb[t] - \[Gamma]A[t]*Pa[t]), 0]
Pb1[t_] :=
If[0 < Pb[t] < 1, (-\[Gamma]B[t]*Pb[t] + \[Gamma]A[t]*Pa[t]), 0]

sol = NDSolveValue[{Pa'[t] == Pa1[t], Pb'[t] == Pb1[t], Pa == 0.5,
Pb == 0.5}, {Pa[t], Pb[t]}, {t, 0, 3/0.000075}];
{Plot[sol, {t, 0, 1/0.000075}, PlotLegends -> {Pa, Pb},
PlotRange -> All],
Plot[sol, {t, 0, 2}, PlotLegends -> {Pa, Pb}, PlotRange -> All]} 