How can I make assumptions in NDSolve?

Question

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?

Answer 1

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 &]

Answer 2

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 = 
   NDSolve`ProcessEquations[{z'[t]^2 == z[t], z[0] == 1}, {z}, {t, 0, 
     10}];

There are two states, one for each solution:

NDSolve`ProcessSolutions[states[[1]], "Forward"]

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

NDSolve`ProcessSolutions[states[[2]], "Forward"]

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

Pick the desired one:

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

Integrate and retrieve the solution:

NDSolve`Iterate[state, 
 NDSolve`SolutionDataComponent[state@"VariableRanges", "Time"]]
sol = NDSolve`ProcessSolutions[state]

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

Third way, using an internal hack:

Internal`InheritedBlock[{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).