NDSolve::ndcf: Repeated convergence test failure. How to solve?

Question

I am trying to simulate a system of $n$ pendulums with some friction in Mathematica 9. This is the code I am using:

nPendulos = 3;
tiempoFinal = 20;
fps = 30;
g = 9.81;
r = 0.5;
x[0][t_] := 0
y[0][t_] := 0
CondicionesIniciales = 
           Join[Table[x[n][0] == n, {n, 1, nPendulos}], Table[y[n]'[0] == 0, {n, 1, nPendulos}]];
Restricciones = Table[(x[n][t] - x[n-1][t])^2 + (y[n][t] - y[n - 1][t])^2 == 1, {n, 1, nPendulos}];
EqNewton = Join[Table[x[n]''[t] == λ[n][t]     (x[n][t] - x[n - 1][t]) - 
                      λ[n + 1][t] (x[n + 1][t] - x[n][t]) - r x[n]'[t], {n, 1, nPendulos - 1}], 
               Table[y[n]''[t] == λ[n][ t] (y[n][t] - y[n - 1][t]) - λ[n + 1][t] 
                                  (y[n + 1][t] - y[n][t]) - g - r y[n]'[t], {n, 1, nPendulos - 1}], 
{x[nPendulos]''[t] == λ[nPendulos][ t] (x[nPendulos][t] - x[nPendulos - 1][t]) - r x[nPendulos]'[t]}, 
{y[nPendulos]''[t] == λ[nPendulos][t] (y[nPendulos][t] - y[nPendulos - 1][t]) - g 
                     - r y[nPendulos]'[t]}];
Vars = Flatten@Table[{λ[n], x[n], y[n]}, {n, 1, nPendulos}];

Sol = NDSolve[{EqNewton, Restricciones, CondicionesIniciales}, Vars, 
              {t, 0, tiempoFinal}, AccuracyGoal -> 2, PrecisionGoal -> 2, MaxStepSize -> 0.01, 
              Method -> {"IndexReduction" -> {True, "ConstraintMethod" -> {"Projection", 
                                                "Invariants" -> Restricciones}}}];

It works fine when the friction is high or nPendulos (number of pendulus) is low. But for example for nPendulos = 4 and r = 0.5 or nPendulos = 3 and r = 0.15, I get things like

NDSolve::ndcf: Repeated convergence test failure at t == 0.93296875`; unable to continue.

NDSolve::ndsz: At t == 2.5562500630525498`, step size is effectively zero; singularity or stiff system suspected.

I am almost sure that the physics behind the system is right, because the results for example when nPendulos = 2 or 3 are nice (see https://dl.dropboxusercontent.com/u/35192406/3_con_rozamiento.gif with friction or https://dl.dropboxusercontent.com/u/35192406/test2.gif with no friction)

Why is NDSolve failing? How can I make it work?

Using Subscript[] makes the code very unpleasant to read here. — Dr. belisarius, Nov 25 '13 at 04:53
Read point number 5 here http://mathematica.stackexchange.com/a/18395/193 — Dr. belisarius, Nov 25 '13 at 05:29
I think it's an issue with stiffness. Using the "StiffnessSwitching" option for the problem with 3 pendulums and 0.15 friction got me a reasonably good looking solution, without any messages. I don't have anything near an answer for you, though. It was mostly poking at options at random, many of which seem to have changed in Mathematica 9.0. — Pillsy, Dec 06 '13 at 01:37
See the n-pendulum example: http://reference.wolfram.com/mathematica/tutorial/NDSolveDAE.html#1781598843 — Michael E2, Dec 06 '13 at 02:04
@MichaelE2 The "IndexGoal" -> 0 option seems to be the crucial one for making the solver succeed in this instance. Everything else seems dispensable. It worked fine without "StiffnessSwitching" or "EquationSimplification" -> "Residual". — Pillsy, Dec 06 '13 at 14:48
@Pillsy Yes, "IndexGoal" -> 0 seems key. Setting gravity = 1 as in the example helps a lot, too. But I would not always get a complete solution depending on how I changed the parameters. Also solutions found with different settings did not always agree very well. I didn't think I found a definitive answer, and the best I could do was give the OP something to play with. — Michael E2, Dec 06 '13 at 15:40

score 2 · Accepted Answer · answered Nov 03 '14 at 23:24

Ok dude, if you want to get convergence on this type of stiffness you need a "more well posed problem", you can't have too much pendulus and little friction at the same time,for avoid confusion this is just a numerical issue. In other words if you want more Pendulos take there biggers! The problem converge on n=4 if r=10.. I personally think this is a numerical limit with 64bit and NDSolve. Engineer advice: model it bigger and scale down smaller.

nPendulos = 4;
tiempoFinal = 20;
fps = 30;
g = 9.81;
r = 10;
x[0][t_] := 0
y[0][t_] := 0
CondicionesIniciales = 
  Join[Table[x[n][0] == n, {n, 1, nPendulos}], 
   Table[y[n]'[0] == 0, {n, 1, nPendulos}]];
Restricciones = 
  Table[(x[n][t] - x[n - 1][t])^2 + (y[n][t] - y[n - 1][t])^2 == 
    1, {n, 1, nPendulos}];
EqNewton = 
  Join[Table[
    x[n]''[t] == \[Lambda][n][
        t] (x[n][t] - x[n - 1][t]) - \[Lambda][n + 1][
        t] (x[n + 1][t] - x[n][t]) - r x[n]'[t], {n, 1, 
     nPendulos - 1}], 
   Table[y[n]''[
      t] == \[Lambda][n][
        t] (y[n][t] - y[n - 1][t]) - \[Lambda][n + 1][
        t] (y[n + 1][t] - y[n][t]) - g - r y[n]'[t], {n, 1, 
     nPendulos - 1}], {x[nPendulos]''[
      t] == \[Lambda][nPendulos][
        t] (x[nPendulos][t] - x[nPendulos - 1][t]) - 
      r x[nPendulos]'[t]}, {y[nPendulos]''[
      t] == \[Lambda][nPendulos][
        t] (y[nPendulos][t] - y[nPendulos - 1][t]) - g - 
      r y[nPendulos]'[t]}];
Vars = Flatten@Table[{\[Lambda][n], x[n], y[n]}, {n, 1, nPendulos}];
<< DifferentialEquations`NDSolveUtilities`
Sol = NDSolve[{EqNewton, Restricciones, CondicionesIniciales}, 
  Vars, {t, 0, tiempoFinal}, AccuracyGoal -> 8, PrecisionGoal -> 10, 
  MaxStepSize -> 0.001, 
  Method -> {"IndexReduction" -> {Automatic, 
      "ConstraintMethod" -> {"Projection", 
        "Invariants" -> Restricciones}}}, SolveDelayed -> True]
StepDataPlot[Sol]

NDSolve::ndcf: Repeated convergence test failure. How to solve?

1 Answers1