Notice removed Draw attention by Community
    Bounty Ended with no winning answer by Community
    Notice added Draw attention by José D.
    Bounty Started worth 150 reputation by José D.
    Notice removed Draw attention by Community
    Bounty Ended with no winning answer by Community
    Notice added Draw attention by José D.
    Bounty Started worth 50 reputation by José D.
5 Code make-up added.
source | link

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

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

I am almost sure that the physics behind the system is right, because the results for example when nPendulos = 2nPendulos = 2 or 33 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 NDSolveNDSolve failing? How can I make it work?

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

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:

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?

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

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

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?

4 deleted 97 characters in body
source | link
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[
Join[Table[x[n]''[t] == λ[n][t]  Table[x[n]''[
   (x[n][t] - x[n - t]1][t]) ==- \[Lambda][n][
        t] (x[n][t] - x[n - 1][t]) - \[Lambda][n + 1][
      λ[n + t]1][t] (x[n + 1][t] - x[n][t]) - r x[n]'[t], {n, 1, 
     nPendulos - 1}], 
   Table[y[n]''[
      t] == \[Lambda][n][
     Table[y[n]''[t] == λ[n][ t] (y[n][t] - y[n - 1][t]) - \[Lambda][nλ[n + 1][1][t] 
        t]                          (y[n + 1][t] - y[n][t]) - g - r y[n]'[t], {n, 1, 
     nPendulos - 1}], {x[nPendulos]''[
      t]{x[nPendulos]''[t] == \[Lambda][nPendulos][
       λ[nPendulos][ t] (x[nPendulos][t] - x[nPendulos - 1][t]) - 
      r x[nPendulos]'[t]}, {y[nPendulos]''[
{y[nPendulos]''[t] == λ[nPendulos][t] (y[nPendulos][t] - y[nPendulos t]- ==1][t]) \[Lambda][nPendulos][
- g  
      t] (y[nPendulos][t] - y[nPendulos - 1][t]) - g - 
       - r y[nPendulos]'[t]}];
Vars = Flatten@Table[{\[Lambda][n]λ[n], x[n], y[n]}, {n, 1, nPendulos}]; 

Sol = NDSolve[{EqNewton, Restricciones, CondicionesIniciales}, Vars, 
   Vars,           {t, 0, tiempoFinal}, AccuracyGoal -> 2, PrecisionGoal -> 2, 
   MaxStepSize -> 0.01, 
              Method -> {"IndexReduction" -> {True, "ConstraintMethod" -> {"Projection", 
       "ConstraintMethod" -> {"Projection", 
                                       "Invariants" -> Restricciones}}}];
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] == \[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}];
Sol = NDSolve[{EqNewton, Restricciones, CondicionesIniciales}, 
   Vars, {t, 0, tiempoFinal}, AccuracyGoal -> 2, PrecisionGoal -> 2, 
   MaxStepSize -> 0.01, 
   Method -> {"IndexReduction" -> {True, 
       "ConstraintMethod" -> {"Projection", 
          "Invariants" -> Restricciones}}}];
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}}}];
    Tweeted twitter.com/#!/StackMma/status/404877794179842048
3 added 16 characters in body
source | link
nPendulos = 3;

tiempoFinal = 20; fps = 30; g = 9.81; r = 0.5;

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] == \[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}];
Sol = NDSolve[{EqNewton, Restricciones, CondicionesIniciales}, 
   Vars, {t, 0, tiempoFinal}, AccuracyGoal -> 2, PrecisionGoal -> 2, 
   MaxStepSize -> 0.01, 
   Method -> {"IndexReduction" -> {True, 
       "ConstraintMethod" -> {"Projection", 
         "Invariants" -> Restricciones}}}];
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] == \[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}];
Sol = NDSolve[{EqNewton, Restricciones, CondicionesIniciales}, 
   Vars, {t, 0, tiempoFinal}, AccuracyGoal -> 2, PrecisionGoal -> 2, 
   MaxStepSize -> 0.01, 
   Method -> {"IndexReduction" -> {True, 
       "ConstraintMethod" -> {"Projection", 
         "Invariants" -> Restricciones}}}];
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] == \[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}];
Sol = NDSolve[{EqNewton, Restricciones, CondicionesIniciales}, 
   Vars, {t, 0, tiempoFinal}, AccuracyGoal -> 2, PrecisionGoal -> 2, 
   MaxStepSize -> 0.01, 
   Method -> {"IndexReduction" -> {True, 
       "ConstraintMethod" -> {"Projection", 
         "Invariants" -> Restricciones}}}];
2 deleted 1384 characters in body
source | link
1
source | link