Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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

My minimum working example is this:

g = {1, 2}
k = 1.5*IdentityMatrix[2];
d = 2*Sqrt[k];
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 to the RHS:

g = {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, 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.

share|improve this question
    
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

2 Answers 2

up vote 5 down vote accepted

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.}}, <>]}} *)
share|improve this answer

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}]
 ]

soln

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.

share|improve this answer
    
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
1  
@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

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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