Notice removed Draw attention by Community♦ occurred Dec 13 '13 at 22:22 Bounty Ended with no winning answer by Community♦ occurred Dec 13 '13 at 22:22 Notice added Draw attention by José D. occurred Dec 5 '13 at 20:40 Bounty Started worth 150 reputation by José D. occurred Dec 5 '13 at 20:40 Notice removed Draw attention by Community♦ occurred Dec 5 '13 at 19:55 Bounty Ended with no winning answer by Community♦ occurred Dec 5 '13 at 19:55 Notice added Draw attention by José D. occurred Nov 27 '13 at 17:59 Bounty Started worth 50 reputation by José D. occurred Nov 27 '13 at 17:59 5 Code make-up added. edit approved Nov 25 '13 at 19:01 Svend Tveskæg 35744 silver badges1313 bronze badges 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 edited Nov 25 '13 at 16:30 Dr. belisarius 108k1111 gold badges173173 silver badges390390 bronze badges 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 occurred Nov 25 '13 at 7:42 3 added 16 characters in body edited Nov 25 '13 at 5:49 José D. 41566 silver badges1919 bronze badges 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 edited Nov 25 '13 at 5:37 José D. 41566 silver badges1919 bronze badges 1 asked Nov 25 '13 at 4:24 José D. 41566 silver badges1919 bronze badges