# How can I make assumptions in NDSolve?

Suppose I try to solve the following differential equation numerically. (My example is much more complicated, but for simplicity, let's consider this example.)

NDSolve[{z'[t]^2 == z[t], z[0] == 1}, {z}, {t, 0, 10}]


Then, Mathematica spits out the following two answers:

I actually want the second answer. I think I need the condition that z'[0]>0. How can I enforce this condition?

• In this case, you could solve z'[t]==Sqrt[z[t]]. Can you rearrange your example in a similar way? Dec 23, 2022 at 18:51
• That's the reason why I am saying that my example is much more complicated. It is impossible to rearrange like this, because z[t], z'[t], y[t], y'[t] are related by two complicated equations that are expressed by exponential functions and polynomials. It's impossible to express y'[t] and z'[t] analytically in terms of y[t] and z[t]. Dec 23, 2022 at 18:58
• Select[sol, And @@ (Positive[z[#]] & /@ {0, 5, 10}) /. # &] Dec 23, 2022 at 19:20

I do not think you can do this using assumptions. But you could always filter out the solutions that you do not want as follows.

I tried to use FunctionMonotonicity to check but this does not seem to work on interpolation functions.

Manually, you could always obtain the grid values from the interpolation function and use Differences then check that all are >=0 or not. This will tell you which function is the one you want. Differences is the numerical equivalent to D used on analytical functions.

ClearAll[z, t]
ode = z'[t]^2 == z[t];
ic = z[0] == 1;
sol = z /. NDSolve[{ode, ic}, z, {t, 0, 1}];


y1 = sol[[1]]["ValuesOnGrid"] // Chop


dy1 = Differences[y1];
y2 = sol[[2]]["ValuesOnGrid"] // Chop


dy2 = Differences[y2];
AllTrue[dy1, # >= 0 &]
AllTrue[dy2, # >= 0 &]


So it is the second one you want since that is what gives True

Ofcourse the above code can be automated more and made into a function that you pass it the output of NDSolve and it returns back only the one which is Monotonically increasing. This is left as an exercise. But the above is the idea.

Update

Here is a small function to automate this. This can be improved more and made shorter if needed.

filterSolutionBy[sol_List, f_Function] := Module[{res, n, y},
res = Reap@Do[
y = Chop[sol[[n]]["ValuesOnGrid"]];
If[AllTrue[Differences[y], f], Sow[sol[[n]]]]
, {n, 1, Length@sol}
];
Last@res
]


Now you can use it as follows

ClearAll[z, t]
ode = z'[t]^2 == z[t];
ic = z[0] == 1;
sol = z /. NDSolve[{ode, ic}, z, {t, 0, 1}]
filterSolutionBy[sol, # >= 0 &]


filterSolutionBy[sol, # <= 0 &]


One way, if you don't mind the IDA method of numerical integration:

NDSolve[{z'[t]^2 == z[t], z[0] == 1, z'[0] == 1}, {z}, {t, 0, 10},
Method -> {"EquationSimplification" -> "Residual"}]


Another way (see this tutorial):

states =
NDSolveProcessEquations[{z'[t]^2 == z[t], z[0] == 1}, {z}, {t, 0,
10}];


There are two states, one for each solution:

NDSolveProcessSolutions[states[[1]], "Forward"]

{z[0.] -> 1., Derivative[1][z][0.] -> -1.}

NDSolveProcessSolutions[states[[2]], "Forward"]

{z[0.] -> 1., Derivative[1][z][0.] -> 1.}


Pick the desired one:

{state} =
Pick[states,
Positive[
z'[0.] /. NDSolveProcessSolutions[#, "Forward"] & /@ states]];


Integrate and retrieve the solution:

NDSolveIterate[state,
NDSolveSolutionDataComponent[state@"VariableRanges", "Time"]]
sol = NDSolveProcessSolutions[state]

{z -> InterpolatingFunction[{{0., 10.}}, ...]


Third way, using an internal hack:

InternalInheritedBlock[{Solve},
Unprotect@Solve;
Solve[eq_, v_, opts___] /; ! TrueQ[$$in] := Block[{$$in = True},
Solve[eq, v,
Assumptions -> v ∈ PositiveReals && z ∈ Reals,
opts]];
Protect@Solve;
NDSolve[{z'[t]^2 == z[t], z[0] == 1}, {z}, {t, 0, 10}]
]


The trouble with using Assuming[...] etc. without the hack is that z'[t] is rewritten as a Module variable with an unpredictable module number appended to it (for example, z\$64135).