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

Question

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[2];
x0 = {1, 3};
v0 = {0, 0};

soln = NDSolve[{x''[t] == k.x[t], x[0] == x0, x'[0] == 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[2];
x0 = {1, 3};
v0 = {0, 0};    
soln = NDSolve[{x''[t] == k.x[t] + g, x[0] == x0, x'[0] == 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.

Answer 1

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[2];
d = 2*Sqrt[k];
x0 = {1, 3};
v0 = {0, 0};
soln = NDSolve[{x''[t] == k.x[t] + g[t], x[0] == x0, x'[0] == v0}, x, {t, 0, 2}]
(* {{x -> InterpolatingFunction[{{0., 2.}}, <>]}} *)

Answer 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[2];
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.