# NDSolve for complex variables

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 = {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][];
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, Subscript[x, i]' == 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). 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.

• Your code works fine in V10.4.1. What version of Mathematica are you using? – Michael E2 Jun 28 '16 at 23:29
• The solutions I get are real-valued, not complex. – Michael E2 Jun 28 '16 at 23:47
• 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. – publius10 Jun 29 '16 at 0:09
• 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. – Michael E2 Jun 29 '16 at 0:35

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


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}] • 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. >> – publius10 Jun 29 '16 at 3:12