5
$\begingroup$

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:

enter image description here

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

$\endgroup$
3
  • $\begingroup$ In this case, you could solve z'[t]==Sqrt[z[t]]. Can you rearrange your example in a similar way? $\endgroup$
    – mikado
    Dec 23, 2022 at 18:51
  • $\begingroup$ 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]. $\endgroup$ Dec 23, 2022 at 18:58
  • $\begingroup$ Select[sol, And @@ (Positive[z[#]] & /@ {0, 5, 10}) /. # &] $\endgroup$
    – Bob Hanlon
    Dec 23, 2022 at 19:20

2 Answers 2

5
$\begingroup$

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}];

Mathematica graphics

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

Mathematica graphics

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

Mathematica graphics

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

Mathematica graphics

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

Mathematica graphics

filterSolutionBy[sol, # <= 0 &]

Mathematica graphics

$\endgroup$
5
$\begingroup$

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).

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.