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 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?

share|improve this question
3  
Using Subscript[] makes the code very unpleasant to read here. –  belisarius Nov 25 '13 at 4:53
    
@belisarius What can I use instead? –  Trollkemada Nov 25 '13 at 5:12
    
Perhaps x[n] instead of Subscript[x, n] –  belisarius Nov 25 '13 at 5:23
2  
Done! I have modified it –  Trollkemada Nov 25 '13 at 5:37
1  
See the n-pendulum example: reference.wolfram.com/mathematica/tutorial/… –  Michael E2 Dec 6 '13 at 2:04

1 Answer 1

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

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.