NDSolveValue: lost solutions?

Question

Here a very easy examplary ode-problem with two solutions:

sol = NDSolve[{x'[t]^2 == x[t]^2, x[0] == 1}, x, {t, 0, 1}];
Plot[x[t] /. sol, {t, 0, 1}, PlotRange -> {0, Automatic}]

as expected.

If I try to solve the same problem with NDSolveValue

X = NDSolveValue[{x'[t]^2 == x[t]^2, x[0] == 1}, x, {t, 0, 1}] 
Plot[X[t] , {t, 0, 1}, PlotRange -> {0, Automatic}]

Mathematica only evaluates one solution!!!

I'm wondering..., any idea what could be the reason? Thanks

Answer 1

NDSolveValue is a convenience wrapper around NDSolve, that makes it easy to get the value of the ODE solution (an InterpolatingFunction object). A simple example:

sol = NDSolveValue[{s'[x] - s[x] == 0, s[0] == 1}, s, {x, 0, 1}]
sol[.5]

1.64872

If NDSolveValue were made to return a list of multiple solutions when they exist, then the for consistency, the output should also return a list with a single solution when only one solution exists. That would mean that one would always have to do a part extraction of the result to get one of the solutions, making it far less convenient to use.

One might ask why not have just a single InterpolatingFunction object (not a list) when there is only one solution, but a list of InterpolatingFunction objects when there is more than one solution. That would make NDSolveValue even more inconvenient to use, as one would have to check whether the output were a list to decide if one solution or multiple solutions were returned. Even more problematic, the second argument of NDSolveValue could be a list of solutions. In that case, one would need to know the input of the NDSolveValue in order to determine whether multiple solutions or a single solution with multiple outputs were being returned.

Basically, in order to be a convenience function, NDSolveValue should only return a single solution even when multiple solutions exist. If you need to worry about multiple solutions, use NDSolve.

Answer 2

This sort of post hoc reasoning behind a decision made by a group of developers is a little dangerous. I think a better phrasing of the question might be, What are the advantages of NDSolveValue over NDSolve, given that NDSolveValue returns only one of the solutions computed by NDSolve? Encouraged by @Szabolcs and by the upvotes on my oversimplified comment, here goes....

I think the primary use of NDSolveValue (singular, not plural NDSolveValues) is to return a single function ifn (or function expression ifn[x]) that can be evaluated on values of the independent variable(s). @Carl Woll has already discussed some aspects of the convenience of this, but I would add that this use probably occurs many more times in practice than all other uses put together, because so many useful ODEs are of the form $y^{(n)} = F\big(x,y,y',\dots,y^{(n-1)}\big)$, usually just $y'=F(x,y)$ or $y''=F(x,y,y')$. NDSolve will return a single solution to a valid IVP for such an ODE, and NDSolveValue therefore will also return the only solution.

My comment suggested that Values@NDSolve[..] would be equivalent to a "plural" form of NDSolveValue, but that's not exactly right. NDSolveValue is much more general in what it can return, except that its return value will be based on the first solution returned by NDSolve, other solutions being ignored. AFAICT, the following two lines are equivalent:

NDSolveValue[<IVP>, expr, <INTERVAL>]
expr /. First@NDSolve[<IVP>, <VARIABLES>, <INTERVAL>]

It seems the expression expr can be anything. If expr = 2, then 2 is returned; if expr = Plot[Sin[t], {t, 0, 10}], then the plot is returned, warnings aside. As far as Values@NDSolve[..] goes, the equivalence limited:

vnds = Values@ NDSolve[{x'[t]^2 == x[t]^2, x[0] == 1}, {x}, {t, 0, 1}];
ndsv = NDSolveValue[{x'[t]^2 == x[t]^2, x[0] == 1}, {x}, {t, 0, 1}];
First[vnds] === ndsv
(*  True  *)

vnds = Values@ NDSolve[{x'[t]^2 == x[t]^2, x[0] == 1}, x, {t, 0, 1}];
ndsv = NDSolveValue[{x'[t]^2 == x[t]^2, x[0] == 1}, x, {t, 0, 1}];
First[vnds] === ndsv
First@First[vnds] === ndsv
(*
  False
  True
*)

They are equivalent if expr = {x} -- that is, if NDSolveValue returns the value of the list {x}, but not if NDSolveValue returns just the function x.

Instead of focusing on the issue of the lost solutions, one should rather think of the advantages of having some syntactic sugar equivalent to expr /. First@NDSolve[..]. A couple favorites:

ListLinePlot@NDSolveValue[ivp, x, {t,..}]  (* equiv. to Plot[x[t], {t,..}] *)
NSolveValue[{y'[x] == f[x], y[a] == 0}, (* equiv. to NIntegrate[f[x], {x, a, b}] *)
 y[b], {x, a, b}]

Those are aside from the highly useful InterpolatingFunction returned by NDSolveValue[ivp, x, {t, 0, 1}]. I think almost everyone who uses ODEs a lot hates always having to remember to type First in x[t] /. First[sol].

OTOH, if you are analyzing a solution, I find sol = NDSolve[ivp, x,..] convenient for plugging into the ODE, invariants and so forth.

Workarounds: (1) If you want an NDSolveValues[] or NDSolveValueList[], try this:

ClearAll[NDSolveValueList];
NDSolveValueList[sys_, expr_, indeps : {_, _, _} .., rest___] := 
  Module[{sol, deps},
   deps = Internal`ProcessEquations`FindDependentVariables[
     Flatten@{sys}, {indeps}[[All, 1]]];
   sol = NDSolve[sys, deps, indeps, rest];
   (expr /. sol) /; FreeQ[sol, NDSolve]
   ];

I don't know much about Internal`ProcessEquations`FindDependentVariables, so it might not work in all cases. Alternatively, for non-PDE systems, one might use the following, which uses documented functions:

ClearAll[NDSolveValueList];
NDSolveValueList[sys_, expr_, {t_, a_, b_}, opts___?OptionQ] :=
  Module[{state, states},
   states = NDSolve`ProcessEquations[sys, expr, {t, a, b}, opts];
   Table[
       state = states[[i]];
       NDSolve`Iterate[state, {a, b}];
       First@Values@NDSolve`ProcessSolutions@state,
       {i, Length@states}] /; FreeQ[states, $Failed]
   ];

If there's a way to get the time-integration interval {a, b} from state, it does not seem to be documented. If someone finds or knows it, then this could probably be extended to PDEs with ease. The data is inside state.

(2) If you want simply to get a different solution from NDSolveValue, there are two approaches. Both determine the solution by the highest derivative value. One uses the method "EquationSimplification" -> "Residual", and the other differentiates the ODE so that the derivative value becomes an initial condition.

NDSolveValue[
 {x'[t]^2 == x[t]^2,
  {x[0], x'[0]} == Values@
     Solve[{x'[t]^2 == x[t]^2, x[0] == 1} /. t -> 0, {x'[0], x[0]}][[2]]},
 x, {t, 0, 1}, 
 Method -> {Automatic, "EquationSimplification" -> "Residual"}]

NDSolveValue[
 {D[x'[t]^2 == x[t]^2, t],
  {x[0], x'[0]} == Values@
    Solve[{x'[t]^2 == x[t]^2, x[0] == 1} /. t -> 0, {x'[0], x[0]}][[2]]},
 x, {t, 0, 1}]

Answer 3

Look at https://reference.wolfram.com/language/ref/NDSolveValue.html in Possible issues/Multiple Solutions. They are indeed saying that NDSolveValue can't make it.

There are two possible differential équations in ur one : $x'(t)=x(t)$ and $x'(t)=-x(t)$.

X = NDSolveValue[{x'[t] == x[t], x[0] == 1}, x, {t, 0, 1}] 
Plot[X[t] , {t, 0, 1}, PlotRange -> {0, Automatic}]

X = NDSolveValue[{x'[t]== -x[t], x[0] == 1}, x, {t, 0, 1}] 
Plot[X[t] , {t, 0, 1}, PlotRange -> {0, Automatic}]

will give you your two solutions