NDSolve for Complex Algebraic-Differential Equation

Question

Let consider the following complex equations: $$ \frac{i x(s)}{2\pi} - \frac{\log(1+ e^{-y(s)}) - \log(1 + e^{y(s)})}{\dot x(s)} $$ which I will cal eq1 and $$ \frac{2 \left(\text{Li}_2\left(-e^{-y(s)}\right)+\text{Li}_2\left(-e^{y(s)}\right)-\frac{3 \pi ^2}{2}\right) x''(s)}{x'(s)^3}-\frac{y'(s) \log \left(e^{-y(s)}+1\right)-y'(s) \log \left(e^{y(s)}+1\right)}{x'(s)^2}+\frac{i y(s)}{2 \pi } $$ which I will call eq2. (btw, notice that since we are in the complex domain I didn't use on purpose sum/subraction formula for polylogs, since in principle I do not know in which branch I'm in) More specifically

eq1 = -((I*x[s])/(2*Pi)) - (2*(Log[1 - E^(-y[s])] - Log[1 + E^y[s]]))/Derivative[1][x][s]
eq2 = -((I*y[s])/(2*Pi)) + (2*(-Log[1 - E^(-y[s])] + Log[1 + E^y[s]])*Derivative[1][y][s])/Derivative[1][x][s]^2 + (4*(-(Pi^2/3) + PolyLog[2,E^(-y[s])] +PolyLog[2, -E^y[s]])*Derivative[2][x][s])/Derivative[1][x][s]^3

I would like to solve numerically them in the interval $s\in (-1/2,1/2)$ with some initial condition, say $x(0)=0$, $y(0)=0$ and $\dot x(0) = 1$. I try to input them in

NDSolve[{eq1 == 0, eq2 == 0, x[0] == 0, y[0] == 0, x'[0] == 1}, {x, y}, {s, -1/2, 1/2}]

I get the following warnings:

NDSolve::ntdvdae: Cannot solve to find an explicit formula for the derivatives. NDSolve will try solving the system as differential-algebraic equations.
NDSolve::mconly: For the method IDA, only machine real code is available. Unable to continue with complex values or beyond floating-point exceptions.
NDSolve::mconly: For the method IDA, only machine real code is available. Unable to continue with complex values or beyond floating-point exceptions.

I think that the first warns me about the method it will use to try to solve the equations, and the second and the third tell me that in order to use IDA method I have to use real valued functions. In order to deal with this I defined

eq1r = (eq1 /. {x -> (xr[#] + I xi[#] &), y -> (yr[#] + I yi[#] &)}) //Re // ComplexExpand;
eq1i = (eq1 /. {x -> (xr[#] + I xi[#] &), y -> (yr[#] + I yi[#] &)}) //Im // ComplexExpand;
eq2r = (eq2 /. {x -> (xr[#] + I xi[#] &), y -> (yr[#] + I yi[#] &)}) //Re // ComplexExpand;
eq2i = (eq2 /. {x -> (xr[#] + I xi[#] &), y -> (yr[#] + I yi[#] &)}) //Im // ComplexExpand;

and then, when I try to execute

NDSolveValue[{eq1r == 0, eq1i == 0, eq2r == 0, eq2i == 0, xr[0] == 0, xi[0] == 0, yr[0] == 0, yi[0] == 0, xr'[0] == 1, xi'[0] == 0}, {xr, xi, yr, yi}, {s, -1/2, 1/2}]

apart for the usual warning I get some other messages like

NDSolveValue::nrnum1: The function value 2. -3.10849*10^-8 I is not a real number when the arguments are {-1.95313*10^-7,-7.51414*10^-33+0. I,3.84724*10^-26+0. I,-1.95313*10^-7+0. I,1. +0. I,3.10849*10^-8+0. I,-1.19591*10^-33+0. I,1,1}.

which I really cannot understand since I explicitly use Re and Im. I think that maybe I have been too naive.

QUESTION

Do you have any suggestion to numerically solve these equations with Mathematica?

Answer 1

The first equation can be solved with respect to y[s]

Solve[eq1 == 0, y[s]]
(*{{y[s] -> 
   Log[1/2 (-1 + E^((I x[s] Derivative[1][x][s])/(4 \[Pi])) - Sqrt[
       1 - 6 E^((I x[s] Derivative[1][x][s])/(4 \[Pi])) + E^((
        I x[s] Derivative[1][x][s])/(2 \[Pi]))])]}, {y[s] -> 
   Log[1/2 (-1 + E^((I x[s] Derivative[1][x][s])/(4 \[Pi])) + Sqrt[
       1 - 6 E^((I x[s] Derivative[1][x][s])/(4 \[Pi])) + E^((
        I x[s] Derivative[1][x][s])/(2 \[Pi]))])]}}*)

We define two equations on two branches

eq3 = eq2 /. {y[s] -> 
    Log[1/2 (-1 + E^((I x[s] Derivative[1][x][s])/(4 \[Pi])) - Sqrt[
        1 - 6 E^((I x[s] Derivative[1][x][s])/(4 \[Pi])) + E^((
         I x[s] Derivative[1][x][s])/(2 \[Pi]))])], 
   y'[s] -> 
    D[Log[1/2 (-1 + E^((I x[s] Derivative[1][x][s])/(4 \[Pi])) - Sqrt[
         1 - 6 E^((I x[s] Derivative[1][x][s])/(4 \[Pi])) + E^((
          I x[s] Derivative[1][x][s])/(2 \[Pi]))]), s]]};
eq4 = eq2 /. {y[s] -> 
    Log[1/2 (-1 + E^((I x[s] Derivative[1][x][s])/(4 \[Pi])) + Sqrt[
        1 - 6 E^((I x[s] Derivative[1][x][s])/(4 \[Pi])) + E^((
         I x[s] Derivative[1][x][s])/(2 \[Pi]))])], 
   y'[s] -> 
    D[Log[1/2 (-1 + E^((I x[s] Derivative[1][x][s])/(4 \[Pi])) + Sqrt[
         1 - 6 E^((I x[s] Derivative[1][x][s])/(4 \[Pi])) + E^((
          I x[s] Derivative[1][x][s])/(2 \[Pi]))]), s]]};

We exclude the point s = 0 from the solution. We solve the equations eq,e4 on the segments {s,-1/2,-x0}, {s,x0,1/2}, where x0 =10^-10. For example

x0 = 10^-10; X2 = 
 NDSolve[{eq3 == 0, x[x0] == 0 + I 0, x'[x0] == 1 + I 0}, 
  x, {s, x0, 1/2}]

{Plot[Evaluate[ReIm[x[s] /. X2]], {s, x0, 1/2}], 
 Plot[Evaluate[
   ReIm[Log[
      1/2 (-1 + E^((I x[s] Derivative[1][x][s])/(4 \[Pi])) + Sqrt[
         1 - 6 E^((I x[s] Derivative[1][x][s])/(4 \[Pi])) + E^((
          I x[s] Derivative[1][x][s])/(2 \[Pi]))])]] /. X2], {s, x0, 
   1/2}]}