# Solving a constrained vector system with NDSolve

I'm trying to replicate the pendulum DAE example in the Documentation Center (example/ModelConstrainedSystemsAsDAEs) using Mathematica 9's vector support with NDSolve.

So far I have:

deqns = {m X''[t] == \[Lambda][t]-{0, m 9.81}};
constraints = {Norm[X[t]] == 1};
initialConditions = {X[0] == {1, 0}, X'[0] == {0, 1}};
params = {g -> 9.81, m -> 1, l -> 1};

pendulumSol =
First[NDSolve[{deqns, constraints, initialConditions} /.
params, {X, \[Lambda]}, {t, 0, 15},
Method -> {"IndexReduction" -> Automatic}]];


but am getting:

NDSolve::ndincd: "Initial conditions for derivatives of the
function X[t] do not have consistent dimensions."


Is this an issue with NDSolve properly identifying the vector nature of the equations, or have I overlooked something? Any help would be appreciated.

-
Your deqns is not quite right to begin with: where are the X[t] and -m g? But this is not the cause of the error message, however. I'm wondering how can MMA determine the dimentionality of Lamda? –  asterix314 Feb 7 at 12:29
I suppose I thought that Mathematica would assume Lambda would have the same dimensionality as X. (Obviously not the case?) –  bridges2 Feb 7 at 14:25
Lamda could be scalar or a vector from the way you write the equations. How can MMA tell? –  asterix314 Feb 8 at 7:55

I think you have to specify the dimension of X when you solve this kind of equations, mathematica can't just deduce it from the fact that {0,9.81} has dimension 2, however the behaviour of NDSolve is still quite curious, try with:

X[t_] := {x1[t], x2[t]};
vec = {0, 9.81};
conditions1 = {Norm[X[t]]^2 == 1, X[0] == {1, 0}, X'[0] == {0, 1}};
conditions2 = {x1[t]^2 + x2[t]^2 == 1, x1[0] == 1, x2[0] == 0,
x1'[0] == 0, x2'[0] == 1};
deqns = X''[t] == \[Lambda][t] X[t] - vec;
sol = NDSolve[{deqns,conditions2}, {X[t], \[Lambda][t]}, {t, 0, 15},
Method -> {"IndexReduction" -> Automatic}];


If you use conditions2 the output of NDSolve is what should be: plot it with ParametricPlot[Evaluate[X[t] /. sol], {t, 0, 15}] you'll see that it is what you expect from a pendulum, however with conditions1, which is substantially the same, the kernel crashes and mathematica gives no output.

Note that deqns is written using vectors (like conditions1) in both cases but doesn't make the system crash, I can't explain why this happens only with initial conditions or constraints.

EDIT:

conditions3={x1[t]^2 + x2[t]^2 == 1, X[0] == {1, 0}, X'[0] == {0, 1}};
conditions4={X[t].X[t] == 1, X[0] == {1, 0}, X'[0] == {0, 1}};


work as fine as conditions2, so the problem must be with the function Norm.

EDIT 2:

As asterix314 noted my first answer was essentialy a trick to avoid the use of NDSolve with symbolic vectors like X[t] since I think that this is in general a more robust approach: NDSolve presents many bugs when dealing with vectors, in the solution of this system we encounter two of them:

1) NDSolve can't handle equations with both symbolic vectors (like X[t]) and lists (like {0, 9.81}) this is because the expression is evaluated BEFORE NDSolve starts operating on it, for example: X''[t]==X[t]+{a,b} is immediately evaluated to X''[t]=={X[t]+a,X[t]+b}, only after the evaluation NDSolve looks at the initial conditions and realizes that X[t] is a vector.

Obviously now on the lhs we have a 2-dimensional vector while on the rhs we have a nested list, because X[t] is a vector itself. This triggers an error, to avoid it I usually define the constant vector in this way:

c[t_?NumericQ]:={a,b}
NDSolve[{X''[t]==X[t]+c[t],X[0]=={x01,x02},X'[0]=={x'01,x'02}},X,{t,0,10}]


Now this works, why? because mathematica doesn't know c[t] is a list until he tries to numerically evaluate it thanks to the conditional pattern _?NumericQ which checks its argument (maybe this is not the most efficient way to go, but it's the only thing that came to my mind).

2) If you try to solve the OP's equation with this trick you'll notice that NDSolve still issues the same error: I think it happens because of the 1-dimensional condition Norm[X[t]]==1, however I don't think there's a real solution for this: when using NDSolve with symbolic vectors all the initial conditions regarding them must be vectors, at least to my experience. It is a limit of NDSolve and this is why I suggest to define vectors as lists and then solving the equation like I proposed at the beginning of this answer.

-
But your answer does not use the vector capabilities of NDSolve -- X will be reduced to {x1, x2} before NDSolve sees them. –  asterix314 Feb 8 at 7:52
@asterix314 {x1,x2} is a vector, but I got your point, I think it has something to do with the vector {0, 9.81}: X[t]-{0,9.81} is evaluated to {X[t], X[t]-9.81} before the equation is solved, that is, before mma sees that X[t] is a vector. –  John Feb 8 at 10:20
I mean the solution will be for x1 and x2 as scalars, not for X as a vector. The question is how to solve for X in and of itself as a vector. –  asterix314 Feb 8 at 15:18
Yes, this is a known problem that happens when you add known constant constant vectors (which NDSolve treats as lists). I will edit my post with a workaround, however I think the solution I presented is the simplest one, after all in this case vectors are used just to simplify the notation. –  John Feb 8 at 15:41