# Differential Equation in Complex Plane

I asked this question 8 years ago, but the answer of @Jens was not satisfactory and it didn't reproduce the plot exactly. I'm wondering if new versions of Mathematica can handle it.

• V 12.3.1 can NDSolve it up to t=0.42 as is. Not to t=10 as you have,. !Mathematica graphics You can try it online if you do not have current version installed? online version is free to try and uses latest version. Commented Nov 13, 2021 at 0:02
• The figure ignores the singularity in the vector field. Mathematica won't, and probably shouldn't, it. Commented Nov 13, 2021 at 1:01
• @Samane This equation can be solved exactly in terms of elliptic functions. If I find time in 24 hours I'll provide another answer, which I believe will be more clarifying. Commented Nov 13, 2021 at 4:05
• See my answer how to do it numerically with little effort. Commented Dec 11, 2021 at 21:13

As noted in a comment by eldo on the earlier question, the ODE can be solved symbolically

DSolve[x'[t] == (1 + (I x[t])^3)^0.5, x[t], t][[1, 1]]

(* x[t] -> InverseFunction[Hypergeometric2F1[1/3, 1/2, 4/3, I #1^3] #1 &][t + C[1]] *)


The inverse is

t + C[1] == Hypergeometric2F1[1/3, 1/2, 4/3, I z^3] z


where I have replaced x[t] by z. The real part of this equation simply defines t + Re[C[1]], which is of no interest. The imaginary part defines the desired curves (parameterized by Im[C[1]]), which can be plotted as follows.

Begin by noting that Hypergeometric2F1[1/3, 1/2, 4/3, I z^3] z has branch cuts at I z^3 == w, where 1 < w < Infinity. The branch cuts are at

bc = Solve[I z^3 == w, z] // Flatten
(* {z -> I w^(1/3), z -> -(-1)^(1/6) w^(1/3), z -> -(-1)^(5/6) w^(1/3)} *)


and the corresponding branch points at w = 1. When a curve crosses a branch cut, its parameter Im[C[1]] changes value. With z -> zr + I zi, the relationship between zr and zi on the branch cuts is

ComplexExpand[Im[I (zr + I zi)^3/zr]] == 0;
Solve[%, zr] // Flatten
(* {zr -> -Sqrt[3] zi, zr -> Sqrt[3] zi} *)


Based on this relationship, the following function computes the change in the curve parameters as they cross the second or third branch cut. (The first branch cut plays no role in the desired figure.)

f[b_] := Module[{t}, t = FindRoot[
Im[Hypergeometric2F1[1/3, 1/2, 4/3, I ((-Sqrt[3] + .9999 I) zi)^3]
(-Sqrt[3] + .9999 I) zi] + b, {zi, -1}];
Im[Hypergeometric2F1[1/3, 1/2, 4/3, I ((-Sqrt[3] + 1.0001 I) zi)^3]
(-Sqrt[3] + 1.0001 I) zi] /. t]


If, for instance, parameters {-1/2, 0, 1/2} are chosen for curves above the branch cuts, then the corresponding parameters below the branch cuts are

f /@ {1/2, 0, -1/2}
(* {-0.902182, -1.40218, -1.90218} *)


If the curve connecting the two conjugate branch points identified above also is desired, its parameter is

a = Im[Hypergeometric2F1[1/3, 1/2, 4/3, I z^3] z /. bc[[3]]] /. w -> 1
(* -((Sqrt[Pi] Gamma[4/3])/(2 Gamma[5/6])) *)


Combining these parameters then gives the desired plot.

ctr = Join[{a}, {1/2, 0, -(1/2)}, f /@ {1/2, 0, -(1/2)}];
ComplexContourPlot[Im[Hypergeometric2F1[1/3, 1/2, 4/3, I z^3] z],
{z, -4 - 6 I, 4 + 2 I}, Contours -> ctr, ContourShading -> None,
PlotPoints -> 500]


Incidentally, I just noticed by comparing numerical values that the change in curve parameters as they cross the branch cuts is equal to 2 a. Undoubtedly, as symbolic proof exists, but I do not have time to search for it.

• @Samane I just reduced the gap significantly.. Aligning the curves better should not be hard. Commented Nov 13, 2021 at 2:42

Artes mentions in a comment that the ODE $$\mathrm dx/\mathrm dt=\sqrt{1 + (i x)^3}$$ is solvable in terms of elliptic functions.

In particular, letting $$x(t)=4i g(t)$$ gives the equivalent nonlinear ODE

$$g^\prime(t)^2=4g(t)^3-1/16$$

which is readily recognized to have the same form as the Weierstrass ODE, and is thus solvable in terms of the Weierstrass elliptic function, built-in as WeierstrassP[]. Skipping the algebra (I have to leave some work for you to do, after all :P), the solution with the initial condition $$x(0)=1-5i$$ can be expressed as

4 I WeierstrassP[t - InverseWeierstrassP[(1 - 5 I)/(4 I), {0, 1/16}], {0, 1/16}]


Using the additional result (cf. the equianharmonic case):

WeierstrassHalfPeriods[{0, 1/16}]
{Gamma[1/3]^3/(2 2^(1/3) Pi), (E^((I Pi)/3) Gamma[1/3]^3)/(2 2^(1/3) Pi)}


allows us to make the following plot:

ParametricPlot[ReIm[4 I WeierstrassP[t - InverseWeierstrassP[(1 - 5 I)/(4 I), {0, 1/16}],
{0, 1/16}]], {t, 0, Gamma[1/3]^3/(2^(1/3) Pi)}]


You can do the whole thing numerically with little effort. Therefore differentiate equation and eliminate the square root term.

d0 = Derivative[1][x][t] - Sqrt[1 - I x[t]^3];

d1 = D[Derivative[1][x][t] - Sqrt[1 - I x[t]^3], t] // Together //
Numerator // Simplify;

eli = Eliminate[{d0 == 0, d1 == 0} /. Sqrt[1 - I x[t]^3] -> aa, aa]

Derivative[1][x][t] == Sqrt[1 - I x[t]^3] /. t -> 0

tab = Table[
xsol[j] =
x /. First@
NDSolve[{eli, Derivative[1][x][0] == Sqrt[1 - I x[0]^3],
x[0] == 1/10 + j I}, x, {t, 0, 10},
MaxSteps -> 30000], {j, {-5, -2, -1, -2/3, -1/3}}];

ParametricPlot[
Evaluate[Table[{Re[xsol[j][t]],
Im[xsol[j][t]]}, {j, {-5, -2, -1, -2/3, -1/3}}]], {t, 0, 10},
PlotRange -> All]


Separating the real and imaginary components.

d[z_, t_] := D[z, t] - Sqrt[1 + (I z)^3]
ode = d[x[t] + I y[t], t] // ComplexExpand
arg = 1 - I (x[t] + I y[t])^3;
arg0 = arg // ComplexExpand
ode0 = ode /. {Arg[arg] -> ArcTan[Re[arg0], Im[arg0]]}
ode1 = ode0 /. {I -> -I};
oder = (ode0 + ode1)/2;
odei = (ode0 - ode1)/(2 I);
tmax = 1.9;
solxy = NDSolve[{oder == 0, odei == 0, x[0] == 0, y[0] == 1/2}, {x, y}, {t, 0, tmax}]
ParametricPlot[Evaluate[{x[t], y[t]} /. solxy], {t, 0, tmax}, AspectRatio -> 1/2]