4
$\begingroup$

Consider the following code

NDSolveValue[{
  x'[t] == x[t]^a,
  x[0] == 0
  }, x, {t, 0, 10}]

which solves the equation $x'(t)=x(t)^a$, for some exponent $a$, with $x(0)=0$. If I want to solve two identical equations, I can simply add a new initial condition

NDSolveValue[{
  x'[t] == x[t]^a,
  x[0] == {0, 0}
  }, x, {t, 0, 10}]

I want to solve this system when $a=1/2$ (square root), however this doesn't work for more than one equation and I get, among others, the error

"NDSolveValue::nbnum1: The function value {0,0}==0 is not True or False when the arguments are {0.,{0.,0.},{0.,0.}}."

Why is that? If I pick any $a\geq 1$ it works as intended. Why is it suddenly interpreting the right hand side as a scalar for $a<1$? In the one equation case, square root works. Any ideas?

Edit and motivation: I want to mimic the following code

NDSolveValue[{x1'[t] == Sqrt[x1[t]], x2'[t] == Sqrt[x2[t]], 
  x1[0] == x2[0] == 0}, {x1, x2}, {t, 0, 10}]

and the reason I'm using this is because I would like to generate a big list of equations. For $a\geq 1$, at least, the main difference between this approach and mine is that here I get a list of two interpolating functions with scalar outputs, and in my case I get one interpolating function with 2 output dimensions.

In other words, what is the best way to solve $$ x_1'(t)=\sqrt{x_1(t)}, x_2'(t)=\sqrt{x_2(t)}, \cdots,x_n'(t)=\sqrt{x_n(t)}\\ $$ for an arbitrary $n$? The main reason I'm interested in this approach is because it's what I use to solve a system like this $$ \begin{pmatrix} x_1'(t)\\ x_2'(t)\\ \vdots\\ x_n'(t) \end{pmatrix}=A f\begin{pmatrix} x_1(t)\\ x_2(t)\\ \vdots\\ x_n(t) \end{pmatrix} $$ where $A$ is some $n\times n$ matrix and $f:\mathbb{R}^n\to\mathbb{R}^n$. In this specific case, $A$ is the identity matrix and $f$ maps the square root.

$\endgroup$
10
  • 1
    $\begingroup$ "If I want to solve two identical equations, I can simply add a new initial condition" - I don't understand that at all, as I don't understand what you mean by x[0] = {0, 0} in your initial conditions. I suspect that NDSolveValue is having the same problem interpreting what you mean. Also, exactly what code generates the error you saw? $\endgroup$
    – MarcoB
    Apr 19 at 14:48
  • $\begingroup$ @MarcoB Here's my motivation: Imagine I want to solve a system of $n$ ODEs that are dependent on each other (say we're measuring the signal of $n$ cells interacting in a tissue, for example). I noticed that a 'compact' way of doing that in Mathematica, when the system of ODEs is really big, is to define an adjacency matrix and treating all functions as a list. For example, see what I did in this question: mathematica.stackexchange.com/questions/224470/… Also, see the edit section, to see what I mean $\endgroup$
    – sam wolfe
    Apr 19 at 14:55
  • $\begingroup$ It sounds like what you want to do is to solve the ODE for many different initial conditions "all at once"? If so, you can do something like: a = 1/2; NDSolveValue[{x'[t] == x[t]^a, x[0] == #}, x, {t, 0, 10}] & /@ Range[0, 1, 0.1] where the Range is just used here to specify all the different initial conditions. $\endgroup$
    – bill s
    Apr 19 at 14:55
  • 1
    $\begingroup$ I don't understand. Even your very first system: NDSolveValue[{x'[t] == x[t]^a, x[0] == 0}, x, {t, 0, 10}] does not work without assigning a value to $a$. $\endgroup$
    – bill s
    Apr 19 at 14:58
  • 3
    $\begingroup$ I think the problems is that Sqrt is a multivalued function. If you use Surd (this gives a real value) instead of Sqrt, then the following works: NDSolveValue[{x'[t] == Surd[x[t], 2], x[0] == {0, 1}}, x, {t, 0, 10}] $\endgroup$ Apr 19 at 15:22
1
$\begingroup$

For a = 1/2, the problem is that the discontinuity processing automatically handles the singularity in the x[t]^a power (Sqrt[x[t]]) by setting up an event when x[t] == 0. Unfortunately, it doesn't notice that x[t] is a vector quantity, so that the event testing causes a slew of errors. The most telling one is NDSolveValue::ecboo, which mentions an "event" and clearly points to some trouble with the discontinuity processing.

A fix is to turn "DiscontinuityProcessing" off:

ClearAll[x, t, a];
Block[{a = 1/2},
 foo = NDSolveValue[{x'[t] == x[t]^a, x[0] == {0, 0}}, x, {t, 0, 10}, 
   Method -> {Automatic, "DiscontinuityProcessing" -> None}]
 ]
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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