# NDSolve returns unevaluated when system gets bigger

I have a bridge, which I want to discretize on $n$ stiff rods. To find a solution I wrote the following code in mathematica:

    n = 4;

variables = Table[ToExpression["φ" <> ToString[i]], {i, 1, n}];

xComponents = Table[Cos[variables[[i]][t]], {i, 1, n}];

yComponents = Table[-Sin[variables[[i]][t]], {i, 1, n}];

endsOfRods = Table[{xComponents[[i]], yComponents[[i]]}, {i, 1, n - 1}];

xTotal = Table[
l (Sum[endsOfRods[[i, 1]], {i, 1, j - 1, 1}] + 1/2 endsOfRods[[j, 1]]), {j,
2, n - 1, 1}];

yTotal = Table[
l (Sum[If[i <= n/2, endsOfRods[[i, 2]], (-1)*endsOfRods[[i, 2]]], {i, 1,
j - 1, 1}] +
If[j > n/2, (-1)*1/2 endsOfRods[[j, 2]], 1/2 endsOfRods[[j, 2]]]), {j,
2, n - 1, 1}];

centreOfMass = Table[{xTotal[[i]], yTotal[[i]]}, {i, 1, Length[xTotal], 1}];

velocityCentreOfMass =
Table[D[centreOfMass[[i]], t].D[centreOfMass[[i]], t], {i, 1,
Length[xTotal], 1}];

kineticEnergy =
Sum[If[i == 1 || i == n, 1/2*(1/3 m l^2)*D[variables[[i]][t], t]^2,
1/2 m velocityCentreOfMass[[i - 1]] +
1/2*(1/12 m l^2)*D[variables[[i]][t], t]^2], {i, 1, n, 1}];

potentialEnergy =
Sum[If[i == 1 || i == n, -m g l 1/2 Sin[variables[[i]][t]] ,
m g yTotal[[i - 1]]], {i, 1, n, 1}](*gravity*)+
Sum[If[i < n/2, 1/2 k (variables[[i]][t] - variables[[i + 1]][t])^2,
If[i == n/2, 1/2 k (variables[[i]][t] + variables[[i + 1]][t])^2,
1/2 k (variables[[i + 1]][t] - variables[[i]][t])^2]], {i, 1, n - 1,
1}](*rotational springs*);

lagrange = kineticEnergy - potentialEnergy;

Nonconservative part

F = {0, d D[(-1)*(yTotal[[n/2 - 1]] - l/2 Sin[variables[[n/2]][t]]), t]};

r = {(xTotal[[n/2 - 1]] + l/2 Cos[variables[[n/2]][t]]), (yTotal[[n/2 - 1]] -
l/2 Sin[variables[[n/2]][t]])};

generalizedForces = Table[F.D[r, variables[[i]][t]], {i, 1, n/2, 1}];

Numerical solution

equations =
Table[If[i <= n/2,
D[D[lagrange, D[variables[[i]][t], t]], t] -
D[lagrange, variables[[i]][t]] - generalizedForces[[i]] == 0,
D[D[lagrange, D[variables[[i]][t], t]], t] -
D[lagrange, variables[[i]][t]] == 0], {i, 1, n - 1, 1}];

AppendTo[equations, (-1)*(yTotal[[n/2]] + l/2 Sin[variables[[n/2 + 1]][t]] +
l Sin[variables[[n]][t]]) == 0];

m = 250(*kg*);
l = 0.75(*m*);
d = 4000(*Ns/m*);
EE = 2.1*10^(11)(*Pa*);
II = 1.5*10^(-6)(*m^4*);
g = 9.81(*m/s^2*);
k = EE II/l;

AppendTo[equations, Table[variables[[i]][t] == 0 /. t -> 0, {i, 1, n, 1}]];
AppendTo[equations,
Table[D[variables[[i]][t], t] == 0 /. t -> 0, {i, 1, n, 1}]];

solution =
NDSolve[equations, Table[variables[[i]][t], {i, 1, n, 1}], {t, 0, 5}];


Now the problem with this code is, that it works perfectly in case of disretization to $n=4$ parts. But changing $n$ to something bigger ($n=10$ for example) ends rather quickly where NDSolve[] only returns the same equations that were given as input.

Are there any other methods or how should one solve a system like this?

PS: This is an upgrade to THIS PROBLEM.

• Seems like you're trying a method of discretization of continuous problem into n arbitrary sections and then solving with continuous functions over them? Something about this approach doesn't seem right to me. If you're doing that you might as well use FEM, FDM, DQM or so many other methods based on matrix solution of algebraic equations. I mean there might be a reason people are not using this method and that could be the same reason why Mathematica can't solve the system when it becomes larger. – MathX Mar 20 '16 at 21:00
• @BehzadNazari You can look at this bridge I mentioned in my OP as if it was a chain with $n$ elongated links. Now imagine, the chain is in a straight horizontal line at time $t=0$. The problem wants me to calculate the movement of the chain links (center of mass) due to the gravity field. There are other forces in the system but this should give you an idea what I am doing. Those variables from the code in the OP are in fact angles between the chain link and horizontal line. Applying equations of motion to each rod should give me a result. I can't see any FEM here, or am I missing something? – skrat Mar 20 '16 at 21:18
• Unfortunately I'm sorry to say that I don't know this subject good enough to give you accurate answers (I switched to more fluids related research 10 years ago) nor do I have the time to read about it or try to learn it because I have my own research to worry about. I'm just suggesting that you look at your approach more objectively. Depending on the step size Mathematica uses, the precision, and the time interval for which you're running the simulation, the equations might become unmanageable. Using larger n values could easily make messy non-banded augmented matrices. Formulation is the key. – MathX Mar 20 '16 at 22:08

This seems to be the same problem (bug?) as mentioned in this comment. Add SolveDelayed->True to NDSolve will resolve the problem. (This option is red, but don't worry. )
BTW, when n >= 10, NDSolve spits out ndsz warning and stops before t reaches 5, and it seems to be the nature of the equation set. If the singularity isn't expected, probably there's something wrong with the definition of equations.