I want to calculate trajectories of an n-pendulum. For the number n=1 it's just a pendulum, for n=2 a double pendulum , for n=3 a triple and so on. The Solver has to solve n equations for the angle. The phi's are the angles and ParametricNDSolve has to solve these. The parameters of ParametricNDSolve are the initial conditions - position and velocity at t=0 of each pendulum. So for n=3 pendulums I have 6 Parameters (position and velocity three times).

The Problem is that I cannot figure out how to define these parameters later and plot the function. I tried the solution from the documentation, but I cannot adapt it to my subscripted functions (I guess - maybe the problem lies somewhere else).

(The whole thing works if I use NDSolve and define the starting conditions beforehand.)

my code:

n = 1; (*number of pendulums chained to each other*)

g = 10; 
α =  0; 
tend = 50; 
(* g is gravity, alpha  is friction, tend is end-time of solve *)

m = Table[1, {i, n}]; (* masses of pendulums *)
l = Table[1 + 0 i, {i, n}]; (* lengths of pendulume *)

r[i_, t_] := Sum[l[[j]]*{Sin[ Subscript[ϕ, j][t]], -Cos[Subscript[ϕ, j][t]]}, {j, 1, i}]; (* positions of pendulums *)
T = Sum[m[[i]]*D[r[i, t], t].D[r[i, t], t]/2, {i, 1, n}];(* kinetic energy *)
U = Sum[m[[i]]*r[i, t].{0, g}, {i, 1, n}];(* potential energy *)
F = Sum[1/2*α*m[[i]]*D[r[i, t], t].D[r[i, t], t], {i, 1, n}]; (* friction term *)
L = T - U // Simplify // Evaluate (* lagrangian *)

sol = ParametricNDSolve[
  Table[{D[D[L, Subscript[ϕ, i]'[t]], t] - 
      D[L, Subscript[ϕ, i][t]] == -D[F, 
       Subscript[ϕ, i]'[t]], 
    Subscript[ϕ, i][0] == Subscript[ϕ0s, i], 
    Subscript[ϕ, i]'[0] == Subscript[ϕ0v, i]}, {i, n}], 
  Table[Subscript[ϕ, i][t], {i, 1, n}], {t, 0, tend + 1}, 
  Flatten[Table[{Subscript[ϕ0s, i], Subscript[ϕ0v, i]}, {i, 
     n}], 1]]

(* ϕ0s is start value, ϕ0v is start velocity *)

f10 = Subscript[ϕ, 1][1, 0][t] /. sol
Plot[f10[t], {t, 0, 1}]

Help would be much appreciated. For anyone interested, I'm trying to find cyclic solutions of the pendulum.


Example for n=3

n = 3;(*number of pendulums chained to each other*)g = 10;
\[Alpha] = 0;
tend = 50;
(*g is gravity,alpha is friction,tend is end-time of solve*)

m = Table[1, {i, n}];(*masses of pendulums*)l = 
 Table[1 + 0 i, {i, n}];(*lengths of pendulume*)
r[i_, t_] := 
     Subscript[\[Phi], j][t]], -Cos[Subscript[\[Phi], j][t]]}, {j, 1, 
   i}];(*positions of pendulums*)T = 
 Sum[m[[i]]*D[r[i, t], t].D[r[i, t], t]/2, {i, 1, 
   n}];(*kinetic energy*)U = 
 Sum[m[[i]]*r[i, t].{0, g}, {i, 1, n}];(*potential energy*)F = 
 Sum[1/2*\[Alpha]*m[[i]]*D[r[i, t], t].D[r[i, t], t], {i, 1, 
   n}];(*friction term*)L = 
 T - U // Simplify // Evaluate (*lagrangian*);

sol = ParametricNDSolve[
  Table[{D[D[L, Subscript[\[Phi], i]'[t]], t] - 
      D[L, Subscript[\[Phi], i][t]] == -D[F, 
       Subscript[\[Phi], i]'[t]], 
    Subscript[\[Phi], i][0] == Subscript[\[Phi]0s, i], 
    Subscript[\[Phi], i]'[0] == Subscript[\[Phi]0v, i]}, {i, n}], 
  Table[Subscript[\[Phi], i][t], {i, 1, n}], {t, 0, tend + 1}, 
   Table[{Subscript[\[Phi]0s, i], Subscript[\[Phi]0v, i]}, {i, n}], 1]]

(*\[Phi]0s is start value,\[Phi]0v is start velocity*)

f10 = Subscript[\[Phi], 1][t] /. sol
Plot[f10[1, .1, .1, 1, 0, .1], {t, 0, 1}, 
 PlotLabel -> Row[{"n = ", n}]]


To show several functions we use

Table[f[i] = Subscript[\[Phi], i][t] /. sol, {i, 1, n}]

Plot[Evaluate[Table[f[i][1, .1, .1, 1, 0, .1], {i, 1, n}]], {t, 0, 1},
  PlotLabel -> Row[{"n = ", n}]]


  • $\begingroup$ wow thank you very much! It seems the exact notation of mathematica was a little too complicated for me. I have one more question (since you seem to know mathematica quite well): Is there a way to make the number of parameters that I give to the parametrized solution arbitrary? What I mean is the [1, .1, .1, 1,0,.1]-part in your solution. I'd like to use a vector of length n*2. Right now I have to give each individual parameter to the function. I have to change the expression when I change n. $\endgroup$ – exocortex Feb 28 at 15:35
  • $\begingroup$ Parameters can be set separately, for example par = {1, .1, .1, 1, 0, .1}; Plot[Evaluate[Table[f[i] @@ par, {i, 1, n}]], {t, 0, 1}, PlotLabel -> Row[{"n = ", n}]] $\endgroup$ – Alex Trounev Feb 28 at 16:15
  • $\begingroup$ Thanks a lot! This is TRULY helpful! I cannot upvote your solution unfortunatly! I am quite new here. Do you think I should change my original question a little bit so it might help others more? $\endgroup$ – exocortex Feb 28 at 21:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.