# Adding a constant vector to a vector differential equation seems to break NDSolve. Why?

I'm trying to solve a differential equation that's phrased in terms of matrices and vectors.

My minimum working example is this:

k = 1.5*IdentityMatrix;
x0 = {1, 3};
v0 = {0, 0};

soln = NDSolve[{x''[t] == k.x[t], x == x0, x' == v0}, x, {t, 0, 2}]
Plot[Evaluate[x[t] /. soln], {t, 0, 2}]


It works fine, and plots two lines against time. AFAIU, this is the solution the system of differential equations:

$\left[\begin{array}{c} \ddot{x}_{1}\\ \ddot{x}_{2} \end{array}\right]=\left[\begin{array}{cc} 1.5 & 0\\ 0 & 1.5 \end{array}\right]\left[\begin{array}{c} x_{1}\\ x_{2} \end{array}\right].$

But it seems to break if I add a constant vector g to the RHS:

g = {1, 2}
k = 1.5*IdentityMatrix;
x0 = {1, 3};
v0 = {0, 0};
soln = NDSolve[{x''[t] == k.x[t] + g, x == x0, x' == v0}, x, {t, 0, 2}]
Plot[Evaluate[x[t] /. soln], {t, 0, 2}]


Now this should be solving the system:

$\left[\begin{array}{c} \ddot{x}_{1}\\ \ddot{x}_{2} \end{array}\right]=\left[\begin{array}{cc} 1.5 & 0\\ 0 & 1.5 \end{array}\right]\left[\begin{array}{c} x_{1}\\ x_{2} \end{array}\right]+\left[\begin{array}{c} 1\\ 2 \end{array}\right]$

But instead,I now get the error "NDSolve::ndfdmc: Computed derivatives do not have dimensionality consistent with the initial conditions. >>". I don't understand why this happens. It seems like the dimensionality of the differential equation matches that of the initial conditions, matches that of the added vector.

Why am I getting this error?

EDIT: Thanks for the workarounds! I think this is approaching bug territory. I've filed a complaint here ( http://www.wolfram.com/support/contact/email/ ) and encourage others to do the same.

• Thanks, Nasser. I saw those. I posted this because it seems to be the simplest case, and it seemed weird that it would accept x[t] as a vector in one case, but not another. – BenB Jul 26 '13 at 6:03
• I think this bug has been solved since. At least your second example does not produce the error message in v.11. Be alerted though that vectorial NDSolve is still not entirely reliable. – István Zachar Aug 28 '18 at 9:56

This is a short-coming of how the arguments are evaluated. The symbols k.x[t] is treated as a single term, while g is treated as a list; Plus automatically threads over the list creating a little mess:

x''[t] == k.x[t] + g
(* x''[t] == {1 + {{1.5, 0.}, {0., 1.5}}.x[t], 2 + {{1.5, 0.}, {0., 1.5}}.x[t]} *)


If a 2-vector value is substituted for x[t], this expression becomes a 2x2 matrix; hence the error message. To get around this, make g into a constant function g[t] that is evaluated only when an actual number is plugged into t:

Clear[g];
g[t_?NumericQ] := {1, 2};
k = 1.5*IdentityMatrix;
d = 2*Sqrt[k];
x0 = {1, 3};
v0 = {0, 0};
soln = NDSolve[{x''[t] == k.x[t] + g[t], x == x0, x' == v0}, x, {t, 0, 2}]
(* {{x -> InterpolatingFunction[{{0., 2.}}, <>]}} *)


The ability to recognize vectorial unknowns such as x[t] in NDSolve is a relatively new feature that doesn't seem to work reliably for my applications, either. So I usually find it much safer to do things in a slightly more "old-fashioned" way, by declaring all unknown functions individually using Array. That can be done relatively efficiently and doesn't add a whole lot of typing:

g = {1, 2}
k = 1.5*IdentityMatrix;
d = 2*Sqrt[k];
x0 = {1, 3};
v0 = {0, 0};

With[
{
x = Through[
Array["x", 2][t]]
},
soln = First@
NDSolve[Join[Thread[D[x, t, t] == k.x + g],
Thread[(x /. t -> 0) == x0], Thread[(D[x, t] /. t -> 0) == v0]],
x, {t, 0, 2}];
Plot[Evaluate[x /. soln], {t, 0, 2}]
] In the With block, I define x as vector function with two components using Array. The individual entries are labeled by String "x" to avoid conflicts, but you can also use regular unused symbols like xEntry. The Through command means that each array entry gets the argument [t].

The major additional modification to the equations then is to wrap each == expression in Thread which has the effect of applying the == to each pair of list members on either side of the vector equation. Since I also defined x to have the time-dependence built into each entry, I have to specify the initial conditions with the replacement rule /. t -> 0.

• for large systems that's going to be slow. Think of all the processing that needs to be done. Could you show examples where the vector version does not work reliably? Possibly what you are referring to is the cause of evaluation semantics, or are there other issues? – user21 Jul 26 '13 at 7:06
• @ruebenko another failure example is here. I just don't like to have to hack around such issues on a case-to-case basis. My approach can be used on a more routine basis. For large systems, I'd just pre-process the system some more before invoking NDSolve. – Jens Jul 26 '13 at 16:41
• I am still not quite sure why you call this a hack - this is a consequence of the fact that Plus is Listable. – user21 Aug 2 '13 at 15:57
• @ruebenko Actually, if I had to give the fundamental reason for these issues, I'd probably say it's the fact that NDSolve doesn't hold its first argument unevaluated. If it did, it could potentially analyze the structure of the equations beforehand and avoid treating Dot and similarly Cross as scalars in additions (which is caused by what you said). I refer to "hacking" because we have to hide the intended meaning of g in the original problem to "trick" the evaluator into holding the sum unevaluated. – Jens Aug 2 '13 at 16:55
• So in this sense, my approach is more "honest" semantically. I do of course prefer the shortened notation without Array whenever it works, but if I ask a student to use it, I'd be very worried about encountering these evaluation order issues. Of course they all have a perfectly logical explanation, but it's a distraction from what I really want to achieve. – Jens Aug 2 '13 at 17:02