NDSolve for complex variables

Question

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, 
    n}];
ic[n_] := 
  Table[{Subscript[x, i][0] == 0, Subscript[x, i]'[0] == 0}, {i, n}] //
    Flatten;

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.

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 *)
Clear[l];
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 *)
Clear[neweq];
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 *)
Clear[SOL];
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:

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