5
$\begingroup$

I've come across a problem that requires me to solve a differential equation of the form $$ (x'(s)+iy'(s))^3 = F(x(s),y(s)) $$ for a given starting point $(x(0),y(0))=(x_0,y_0)$, where $F$ is a complex-valued function of the (complex) plane (not necessarily analytical), and I would like to know if there are good ways to convince NSolve to solve it.

Now: yes, I know, this is an ugly system, and yes, I know, there's multiple values of the derivative consistent with that equation at any given point, but it's still perfectly possible to ask for curves that are continuously differentiable and which obey it, and that rules out any kinks in the curve outside of the zeroes of $F(x,y)$. And yes, I know that, even then, there's still three valid starting values of the velocity at the initial point - I want all of them, either one by one or, ideally, in one go.

The first thing I thought to try with this one was to just take the cubic root of both sides and hope for the best, via e.g.

Block[{x0, y0, eqns, soln}, x0 = -1; y0 = 1; eqns = Join[ Thread[{x'[s], y'[s]} == ({Re[#], Im[#]}/Abs[#])&[((x[s] - I y[s])^2)^(1/3)]], Thread[{x[0], y[0]} == {x0, y0}] ]; soln = First@NDSolve[eqns, {x, y}, {s, 0, 3}]; ParametricPlot[ {x[s], y[s]} /. soln , {s, 0, 3} , Frame -> True ] ]

(where $F(x,y) = (x-iy)^2$ is quite representative, with more general $F(x,y)=(x\pm i y)^n$ also interesting), but this runs into trouble with the cube-root's branch cut:

Mathematica graphics

That is, when the curve reaches the $y$ axis, and therefore $F(x,y)$ reaches the negative real axis from below, the cube root $\sqrt[3]{F(x,y)}$ crosses its branch cut, swaps to another root, and the solution switches to another track that it shouldn't be on.

My main goal is to solve this branch-choosing problem in a way that is as quick, simple, clean, built-in, and efficient as possible. (I'll be doing a lot of these solves once I have this sorted.) Some approaches I've tried:

  • If pressed, I can build an Euler's-method integrator where the branch cut is explicitly chosen based on the current value of the velocity. However, I would rather not re-implement stuff that NDSolve is meant to be taking care of for me.
  • So far, NDSolve hasn't much liked it when I try and force the right-hand-side to have conditionals (or otherwise, e.g. powers of $x'+iy'$ inside the cube root cancelled by powers outside it) that depend on the velocity on the right-hand side. For all of the stuff I've tried, NDSolve just flat-out refuses to take it.
  • I also tried feeding the equation itself, $(x'(s)+iy'(s))^3 = F(x(s),y(s))$ with the cubed velocity and all, into NDSolve, but it wasn't particularly keen on that, either. It does do some of the solving, but it throws a rather wide variety of numerical errors and warnings, it is liable to just plain stop (or even worse, backtrack and re-track and backtrack and re-track over the same bit of trajectory), and it just doesn't seem to be a fan of the form.

So: how can I get Mathematica to integrate equations of this form with as minimal an amount of pain as possible?

Improve this question
$\endgroup$
1
  • $\begingroup$ Surd might be helpful. $\endgroup$ Commented Nov 15, 2017 at 1:33

1 Answer 1

Reset to default
4
$\begingroup$

See update below

You could try using a complex dependent value:

sol = NDSolve[{z'[s]^3 == Conjugate[z[s]]^2, z[0] == -1 + I}, z, {s, 0, 3}];

NDSolve::mxst: Maximum number of 180980 steps reached at the point s == 0.9165575857155464`.

This gives the following parametric plots:

ParametricPlot[ Evaluate @ ReIm @ Through @ (Values @ Flatten @ sol)[s], {s, 0, 3}, PlotRange->{{-8, 8}, {-8, 8}} ]

enter image description here

Update

The problem with the blue curve is the branch cut of the cube root along the negative $x$-axis. For instance, note:

(-1. + 10`*^-16I)^(1/3) (-1. - 10`*^-16I)^(1/3)

0.5 + 0.866025 I

0.5 - 0.866025 I

This discontinuity prevents NDSolve from extending the blue line below the $x$-axis. Probably this could be worked around by using an alternate cube root with a branch cut along the positive $x$-axis, but I haven't pursued this idea. Instead, we can control the branch cut being used by taking a cube root of the equation, and inserting cube roots of 1 in appropriate places. Here is this idea in action. First, the naive cube root:

sol1 = NDSolveValue[{z'[s] == (Conjugate[z[s]]^2)^(1/3), z[0] == -1 + I}, z, {s, 0, 3}]; ParametricPlot[ReIm @ sol1[s], {s, 0, 3}]

enter image description here

This reproduces the issue in the OP. To work around this, I will tweak things slightly:

sol1 = NDSolveValue[{z'[s] == Exp[2 Pi I/3] Conjugate[z[s]]^(2/3), z[0] == -1 + I}, z, {s, 0, 3}]; ParametricPlot[ReIm @ sol1[s], {s, 0, 3}]

enter image description here

That's much better. Now for the other two cube roots of 1:

sol2 = NDSolveValue[{z'[s] == Exp[2 Pi I/3] (Conjugate[z[s]]^2)^(1/3), z[0] == -1 + I}, z, {s, 0, 3}]; sol3 = NDSolveValue[{z'[s] == Exp[4 Pi I/3] (Conjugate[z[s]]^2)^(1/3), z[0] == -1 + I}, z, {s, 0, 3}]; ParametricPlot[Evaluate @ ReIm @ {sol1[s], sol2[s], sol3[s]}, {s, 0, 3}]

enter image description here

This looks pretty good!

Improve this answer
$\endgroup$
6
  • $\begingroup$ Was going to suggest this myself... (+1) $\endgroup$ Commented Nov 15, 2017 at 1:54
  • $\begingroup$ This has an equivalent problem with the kink when the curve reaches the negative real axis. This one should be smooth and stay in the third quadrant. $\endgroup$ Commented Nov 15, 2017 at 7:30
  • $\begingroup$ @Emilio, that might be related to the need to take the principal branch of Power[]. I'm guessing the solution you actually want corresponds to a different branch. $\endgroup$ Commented Nov 15, 2017 at 10:51
  • $\begingroup$ @J.M. Frankly, I don't understand why it's encountering any hiccups at that point - the RHS maps to 1, and that's smack in the safest zone of the principal branch of Power[]. The kink in Carl's plot seems to be extraneous, though - if I run that NDSolve, I get the same mxst error, and looking at the internals of the first InterpolatingFunction, it doesn't have any data at all in either the third or the fourth quadrants; the solution just stops at $s\approx 0.916\leftrightarrow y=0$. $\endgroup$ Commented Nov 15, 2017 at 11:32
  • 1
    $\begingroup$ @CarlWoll What version of Mathematica are you using? When I run your code with 11.2 on Windows 10, the blue curve does not appear at all, although the InterpolationFunction for the blue curve (in the second quadrant) certainly exists. $\endgroup$ Commented Nov 16, 2017 at 2:28

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.