7
$\begingroup$

I just discovered the following "works":

NDSolveValue[{ (* Update note: 
     the ODE code below evaluates to 
         y'[t] == {-y[t], -2y[t], -3y[t]}
     before NDSolve[] is called -- this is made clear
     below but maybe worth pointing out at the start *)
   y'[t] == -{1, 2, 3}*y[t], 
   y[0] == {1, 2, 3}}
  , y, {t, 0, 1}] // ListLinePlot

The result returned by NDSolve[] is the vector-valued solution for ${\bf y}=(y_1,y_2,y_3)$ to the system $$\eqalign{ y_1'(t) &= -y_1(t)\,,\ y_1(0)=1 \cr y_2'(t) &= -2y_2(t)\,,\ y_2(0)=2 \cr y_3'(t) &= -3y_3(t)\,,\ y_3(0)=3 \cr }$$

Astute readers may object that the solution cannot satisfy the input ODE, since its actual form is

y'[t] == {-y[t], -2y[t], -3y[t]}

and y[t] and y'[t] are vectors. What it actually satisfies (approximately) is

If[VectorQ[y[t0]] && VectorQ[y'[t0]},
 y'[t0] == {-1, -2, -3} * y[t0]
 ]

That is, you hold up the threading of Times[] until y[t] and y'[t] are explicit vectors. We may test this in code as follows:

wp = MachinePrecision;
sol = NDSolveValue[{
    y'[t] == -{1, 2, 3}*y[t],
    y[0] == {1, 2, 3}}
   , y, {t, 0, 1}, WorkingPrecision -> wp];
And @@ MapThread[
  Block[{y, t, Derivative,
     Internal`$EqualTolerance = (wp + 1)/2.}, (* pg = wp/2 plus tol. *)
    y[t] = #1; y'[t] = #2;
    y'[t] == -{1, 2, 3}*y[t]
    ] &,
  {sol["ValuesOnGrid"], sol'["ValuesOnGrid"]}
  ]
(*  True  *)

My inference is that we can use NDSolve[] in this way to solve uncoupled systems of the form $${y}_k'(t) = (f_1(t,{y}_k),\dots,f_n(t,{y}_k)) \,,\ k=1,\dots,n\,;\quad {\bf y}(t_0) = (a_1,\dots,a_n)$$ with code for the IVP having the form

{y'[t] == {f1[t, y[t]],…,fn[t, y[t]]}, y[0] == {a1,…,an}}

Questions:

  1. Is this usage documented somewhere? (I couldn't find it.) If not, we could document it here as a service to the community, I suppose.
  2. Are there restrictions (conditions that have to be true in order for it to work properly)?
  3. If undocumented, does anyone know the status of this usage? Might it go away in the future?
Improve this question
$\endgroup$
8
  • 1
    $\begingroup$ From documentation of NDSolveValue under "Details and Options": "The c0, dc0, etc. can be lists, specifying that u[x] is a function with vector or general list values." Then under "Scope - Ordinary Differential Equations" third example is "Solve for a vector-valued function:", where we can see even more complicated example of such usage. $\endgroup$ Commented 2 days ago
  • $\begingroup$ This usage is not documented there. I am well aware that NDSolve[] and DSolve[] solve matrix-vector forms of ODEs and have given many answers using such. And I already knew that the above form is equivalent to the corresponding diagonal matrix form with vector solution (duh). But the above is different from the documented form. $\endgroup$ Commented 2 days ago
  • $\begingroup$ @MichaelE2 But has this "undocumented form" any sense? I don't think so. What is NDSolveValue really solving is the form shown in azerbajdzan's answer. It is more a bug of NDSolveValue then a feature to me. On the other hand if nonsensical input is provided should we be expecting any reasonable output? $\endgroup$ Commented 2 days ago
  • $\begingroup$ Ah finally here comes a post about this behavior. I also noticed this here, here and here but never looked into it. @three777 I believe this feature is, in short, a remedy, because NonThreadable isn't introduced until v14.1. $\endgroup$ Commented yesterday
  • 1
    $\begingroup$ Changing any of the coefficients above to zero yields a buggy example: NDSolve[{y'[t] == -{1, 0, 3}*y[t], y[0] == {1, 2, 3}}, y, {t, 0, 1}] emits an error and returns unevaluated. $\endgroup$ Commented yesterday

2 Answers 2

Reset to default
7
$\begingroup$

I've known about the ability of NDSolve[] to handle vector forms of ODEs for years. But the OP concerns a different form.

{
  (* undocumented form *)
  sol1 = NDSolveValue[{
     ode1 = y'[t] == -{1, 2, 3}*y[t],
     ic = y[0] == {1, 2, 3}},
    y, {t, 0, 1}] // ListLinePlot,
  (* documented form *)
  sol2 = NDSolveValue[{
     ode2 = y'[t] == -DiagonalMatrix[{1, 2, 3}] . y[t],
     ic},
    y, {t, 0, 1}] // ListLinePlot} // GraphicsRow

enter image description here

They're equivalent -- no doubt. Indeed, @azerbajdzan shows they're the same -- not surprising. This all seems obvious from the description in the OP of NDSolve[] setting up a diagonal system from an IVP like {ode1, ic}.

Here's what's different: The solution sol1 is not a solution to the equation ode1, but both sol1 and sol2 are solutions of ode2:

ode1 /. t -> 0 /. y -> sol1
ode2 /. t -> 0 /. y -> sol1
ode2 /. t -> 0 /. y -> sol2
(*
{-1., -4., -9.} == {{-1., -2., -3.}, {-2., -4., -6.}, {-3., -6., -9.}}
True
True
*)

Why is that? Well look at ode1, which is the actual input to NDSolveValue[]:

ode1
(*  y'[t] == {-y[t], -2 y[t], -3 y[t]}  *)

If you plug a vector such as y[0.] == {1.,2.,3.} in for y[t], you get a matrix on the right-hand side. As can be seen previous code output, if y is a vector-valued function, so is y'; and further, the left-hand side of ode1 will be a vector and the right-hand side will be a matrix. There is no way to satisfy ode1.

This problem of Plus[] and Times[] threading vectors and matrices in setting up vector/matrix forms of ODEs has a long history on this site, to which I've contributed many answers. (The threading is due to Plus[] and Times[] being Listable.) Long-time community members interested in will know that often the solution has required the construction of ?NumericQ-protect functions to prevent threading. The issue in ode1 is the premature threading in {1,2,3} y[t]. The Times-threading does not occur after y[t] has a vector value as @azerbajdzan seems to think, but before it. In times past, this has yielded an error, with which site visitors have sought help. Indeed, DSolve[{ode1, ic}, y, t] still gives an error, which seems correct since ode1 cannot have a solution.

So if the input with -{1, 2, 3}*y[t] doesn't make sense, unlike the DiagonalMatrix[] form, what is NDSolve[] doing? It is easy to see, if you know the inner workings of NDSolve[] as explained in Advanced Numerical Differential Equation Solving in the Wolfram Language and other NDSolve[] manuals. NDSolve[] is un-threading the right-hand side of the ODE. It is somewhat remarkable, since threading is not strictly invertible. It is also worth remarking, though less remarkable, that it has bugs.

The unthreading converts ode1 to a diagonal system as explained in the OP, although it need not be a linear one as is the simple example in the OP. Consequently, it is equivalent to the explicitly diagonal system in ode2.

For instance, the first example gives the sort of error I'd expect, but the second, which changes a coefficient 3 to a 2, works fine:

NDSolve[{y'[t] == {2 t y[t] - y[t]^3/5, 3 t y[t] + y[t]^2/4}, 
  y[0] == {2, 1}}, y, {t, 0, 1}]

NDSolve::ndfdmc: Computed derivatives do not have dimensionality consistent with the initial conditions.

NDSolve[{y'[t] == {2 t y[t] - y[t]^3/5, 2 t y[t] + y[t]^2/4}, 
  y[0] == {2, 1}}, y, {t, 0, 1}]

(*  {{y -> InterpolatingFunction[{{0., 1.}}, "<>"]}}  *)

The problem seems to go away because of identical terms 2 t y[t] in both components in the second example. Bugs sometimes don't make a lot of sense.

Here is a way to inspect the unthreading:

ones // ClearAll;
ones[a_] := ConstantArray[1, Dimensions@a];
checkRHSFN // ClearAll;
checkRHSFN[rhss_, ics_ : 1] := Module[{state, ode},
   {state} = NDSolve`ProcessEquations[{
      ode = y'[t] == (rhss /. {y[t] -> y[t], y -> y[t]})
      , y[0] == Replace[ics, n_?NumericQ :> n*ones[rhss]]}
     , y, {t, 0, 1}];
   Grid[
    Join[
     {{"ODE", ode}},
     {#, state["NumericalFunction"]@#} & /@ {"FunctionExpression", 
       "ArgumentDimensions", "ResultDimensions"}
     ], Alignment -> Left, Spacings -> {1, Automatic}]
   ];

The right-hand side rhss in checkRHSFN[] can be given without the [t] arguments to y.

Two of the examples we've seen show that NDSolve[] does some work to get the job done. The NumericalFunction's "FunctionExpression" represents the right-hand side of the ODE. It is, in effect, evaluated on numeric quantities of the dimensions shown. This means that the threading happens after the y argument is given a vector value.

checkRHSFN[-{y, 2 y, 3 y}, {1, 2, 3}]
checkRHSFN[{2 t y - y^3/5, 2 t y + y^2/4}, {2, 1}]

output

Here's an example where NDSolve[] fails. The unthreaded expression in the function should be something like t y + {-1, -2, -3} y instead of just t y, and the solutions should decrease (right) instead of increase (left).

checkRHSFN[{t y - y, t y - 2 y, t y - 3 y}, {2, 1, 1/2}]

output

graphs of the incorrect and correct solutions

I'm pretty sure that this example and several others like it are bugs in the unthreading process, but I was wondering whether the proper behavior was documented. Maybe there is some reason why they should be dismissed.

Improve this answer
$\endgroup$
2
$\begingroup$

It is documented in both NDSolve and DSolve.

What you wrote as -{1, 2, 3}*y[t] can be more understandable in the form -{{1, 0, 0}, {0, 2, 0}, {0, 0, 3}} . y[t] as you can see in the code below. nds1 is code from OP, nds2 and ds are its more understandable equivalents.

nds1 = NDSolveValue[{y'[t] == -{1, 2, 3}*y[t], 
   y[0] == {1, 2, 3}}, y[t], {t, 0, 1}]

nds2 = NDSolveValue[{y'[t] == -{{1, 0, 0}, {0, 2, 0}, {0, 0, 3}} . y[t], 
   y[0] == {1, 2, 3}}, y[t], {t, 0, 1}]

%% === %

ds = DSolveValue[{y'[t] == -{{1, 0, 0}, {0, 2, 0}, {0, 0, 3}} . y[t],
   y[0] == {1, 2, 3}}, Element[y[t], Vectors[3]], t]

{nds1, nds2, ds} /. t -> 0.1

enter image description here

Updated Conclusion:

nds1 is not solution to the mathematically nonsensical system of equations y'[t] == {-y[t], -2y[t], -3y[t] but to the system of equation y'[t] == -{{1, 0, 0}, {0, 2, 0}, {0, 0, 3}} . y[t]. It could be a coincidence that NDSolve internally transforms nonsensical input into something meaningful or it is simply a bug of NDSolve.

Original incorrect Addendum:

(true only if y[t] were substituted by vector before the multiplication took place in {1, 2, 3}*y[t])

{1, 2, 3}*y[t] and {{1, 0, 0}, {0, 2, 0}, {0, 0, 3}} . y[t] are numerically identical if y[t] is vector of length 3. NDSolve, being a numerical method, does not care whether it is the former or the latter form (unlike DSolve).

Improve this answer
$\endgroup$
6
  • $\begingroup$ As for DSolve[], try DSolveValue[{y'[t] == -{1, 2, 3}*y[t], y[0] == {1, 2, 3}}, y[t], {t, 0, 1}]. $\endgroup$ Commented 2 days ago
  • $\begingroup$ @MichaelE2 Have you read the addendum at the bottom? {1, 2, 3}*y[t] is not same as {1, 2, 3}.y[t]. $\endgroup$ Commented 2 days ago
  • $\begingroup$ Isn't the addendum wrong? They are not numerically identical. {1, 2, 3}*y[t] is a 3x3 matrix and {{1, 0, 0}, {0, 2, 0}, {0, 0, 3}} . y[t] is a 3-vector. Right? $\endgroup$ Commented 2 days ago
  • 1
    $\begingroup$ y[t] is vector {y1, y2, y3} in both cases. How is {1, 2, 3}*y[t] 3x3 matrix when {1, 2, 3}*{y1, y2, y3} produces vector {y1, 2 y2, 3 y3}??? $\endgroup$ Commented 2 days ago
  • 4
    $\begingroup$ "How is {1, 2, 3}*y[t] 3x3 matrix..." The key point is the evaluation order. If it's y[t]={y1, y2, y3}; {1, 2, 3}*y[t] then it won't be a matrix of course, but what's discussed here is amount to Clear[y]; {1,2,3} y[t]/.y[t]->{y1,y2,y3}. (NDSolve doesn't have HoldAll, HoldFirst, etc. attribute, so something like {1, 2, 3}*y[t] won't hold, it immediately evaluates to {y[t],2y[t],3y[t]}, but NDSolve internally reverses the procedure, this is quite unusual, and that's what Michael E2 wants to talk about in this post. ) $\endgroup$ Commented yesterday

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.