I am trying to simulate a multiple pendulum system with damping as follows:

l = 0.6 (1 + 0.00 I); g = 9.8;
x0 = 0.01;
coord[n_] := Table[Subscript[x, i][t], {i, n}];
vel[n_] := D[coord[n], t];
mass[1] = {x0 Cos[\[Omega] t], 0};
mass[i_] := 
  mass[i - 1] + {l Sin[Subscript[x, i - 1][t]], -l Cos[
      Subscript[x, i - 1][t]]};
ke[i_] := With[{v = D[mass[i + 1], t]}, v.v/2];
pe[i_] := g mass[i + 1][[2]];
L[n_] := Sum[ke[i] - pe[i], {i, n}];
eq[n_] := 
  Table[D[D[L[n], vel[n][[m]]], t] - D[L[n], coord[n][[m]]] == 0, {m, 
ic[n_] := 
  Table[{Subscript[x, i][0] == 0, Subscript[x, i]'[0] == 0}, {i, n}] //

So that eq[n] is the Lagrangian equation of motion for n linked pendula. Now I can get a solution with

SOL[n_] := 
 NDSolve[{FullSimplify[eq[n]] /. \[Omega] -> 1, ic[n]} // Flatten, 
  Table[Subscript[x, i][t], {i, n}], {t, 0, 100}, 
  Method -> {"EquationSimplification" -> "Residual"}]

(try for instance f=SOL[4]). However if I add a little bit of damping by changing l=0.6(1+0.01I), the solver fails with the error

NDSolve::icfail: Unable to find initial conditions that satisfy the residual function within specified tolerances. Try giving initial conditions for both values and derivatives of the functions. >>

So my question is whether there is some method available to NDSolve that I can call to make this work for complex variables.

  • $\begingroup$ Your code works fine in V10.4.1. What version of Mathematica are you using? $\endgroup$
    – Michael E2
    Jun 28, 2016 at 23:29
  • $\begingroup$ The solutions I get are real-valued, not complex. $\endgroup$
    – Michael E2
    Jun 28, 2016 at 23:47
  • $\begingroup$ yes it works as is, but replace the first line with l=0.6(1+0.01I) to make l complex. currently trying your solution with that replacement, it's taking a while. $\endgroup$
    – publius10
    Jun 29, 2016 at 0:09
  • $\begingroup$ Try the Solve step with a symbolic l (Clear[l]). Then substitute l -> 0.6(1+0.01I). The approximate real numbers are bogging Solve down. $\endgroup$
    – Michael E2
    Jun 29, 2016 at 0:35

1 Answer 1


If you solve for the second derivatives, you won't have to use "EquationSimplification" -> "Residual" and things will work ok.

Solving for the second derivatives be faster if you start with exact coefficients. Also, if you clear l, solve for the derivatives, and then substitute a value for l, Solve won't choke on the algebra. The long time it takes is probably why NDSolve advises the use of the option:

NDSolve::ntdv: Cannot solve to find an explicit formula for the derivatives. Consider using the option Method->{"EquationSimplification"->"Residual"}. >>

(*l = 0.6 (1 + 0.00 I) // Rationalize;*)  (* old solution *)
g = 9.8 // Rationalize;
x0 = 0.01 // Rationalize;
(*... rest of OP's code...*)

(* This replaces eq[] in NDSolve; saves the result for ease of reuse *)
neweq[n_, ω0_] := neweq[n, ω0] = 
  Equal @@@ First@ Simplify@
     Solve[eq[n] /. ω -> ω0, Table[Subscript[x, i]''[t], {i, n}]]

(* Here is where we substitute the value for l *)
SOL[n_] := NDSolve[{neweq[n, 1] /. l -> 0.6 (1 + 0.01 I), ic[n]},
  Table[Subscript[x, i][t], {i, n}], {t, 0, 100}]

mysol = First@SOL[4]

Since the solutions are complex, one has to decide what to plot. I chose the magnitude. The plots nearly overlap, so I offset them:

 0.01 Range[n] + Table[Abs@ Subscript[x, i][t], {i, n}] /. mysol // Evaluate,
 {t, 0, 100}]

Mathematica graphics

  • $\begingroup$ thanks, that did it, although depending on the value of ω i use, i get numerical errors. for this example i used ω=1, but for 100 i get Infinity::indet: Indeterminate expression 0. ComplexInfinity encountered. >> $\endgroup$
    – publius10
    Jun 29, 2016 at 3:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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