2
$\begingroup$

I have a system of equations I need to solve involving square roots. I don't always have a good guess what the solutions are, and so I am fortunate that NSolve can handle these equations. I want to be able to try moving the branch cut of the square root function to the positive imaginary axis since I want continuous functions as I cross the positive and negative real axis. I can easily accomplish this as follows:

Define my new square root function as

MySqrt[z_] := If[Re[z] < 0 && Im[z] >= 0, -Sqrt[z], +Sqrt[z]]

Now for some equations I pull out of a hat:

EQNS = {x + MySqrt[y] - 2/MySqrt[z] + z - 1, 
         y - x + MySqrt[x] + 1/MySqrt[z] - 2, 
           MySqrt[x] + 1/MySqrt[y] - 2 y + z};
NSolve[EQNS == 0, {x, y, z}] // Timing

Yields:

  Out[990]={0.92901, 
    {{x -> 0.698929 + 3.41798 I, y -> 1.03321 + 0.455899 I, z -> -0.301342 - 0.075144 I},
     {x -> 2.72138, y -> 1.42833, z -> 0.37027}}}

Now to check this against the standard definition:

EQNS = {x + Sqrt[y] - 2/Sqrt[z] + z - 1, 
       y - x + Sqrt[x] + 1/Sqrt[z] - 2, 
       Sqrt[x] + 1/Sqrt[y] - 2 y + z};
NSolve[EQNS == 0, {x, y, z}] // Timing

Yields:

Out[988]={0.279246, {
{x -> 0.698929 + 3.41798 I, y -> 1.03321 + 0.455899 I, z -> -0.301342 - 0.075144 I}, 
{x -> 0.698929 - 3.41798 I, y -> 1.03321 - 0.455899 I, z -> -0.301342 + 0.075144 I}, 
{x -> 0.698929 + 3.41798 I, y -> 1.03321 + 0.455899 I, z -> -0.301342 - 0.075144 I}}}

Which brings me to my question: Why does my definition for the square root make it so that it takes NSolve so much longer? It is about a factor of 4 here, which isn't a huge deal, but I am solving even more complicated equations with more variables and I am finding the increase to be a factor of 10 or greater. Can anyone suggest a better way to define a new (faster?) square root function? Or make any other suggestions to make this process faster?

$\endgroup$
4
$\begingroup$

Multiplying a number by $e^{\pi i/2} = i$ is equivalent to rotating the phase by $90^\circ$. So, if you pre-multiply your number by $i$ and take the square root, you will be using the same branch cut as the default. Then, to compensate, you will need to multiply by $e^{-\pi i/4}$ to rotate the phase backwards by $45^\circ$. Here is a function that does this:

isqrt[d_] := -Sqrt[I d](-1)^(3/4)

First, let's compare isqrt with MySqrt:

sample = RandomComplex[{-5 - 5 I, 5 + 5 I}, 100];
MinMax @ Abs[isqrt[sample] - MySqrt /@ sample]

{0., 6.28037*10^-16}

So, basically the same (I check the norm of the differences to account for differences in numerics). Next, let's check the timing:

EQNS={x+MySqrt[y]-2/MySqrt[z]+z-1,y-x+MySqrt[x]+1/MySqrt[z]-2,MySqrt[x]+1/MySqrt[y]-2 y+z};
NSolve[EQNS==0,{x,y,z}] //Timing

EQNS={x+isqrt[y]-2/isqrt[z]+z-1,y-x+isqrt[x]+1/isqrt[z]-2,isqrt[x]+1/isqrt[y]-2 y+z};
NSolve[EQNS==0,{x,y,z}] //Timing

{0.783682, {{x -> 0.698929 + 3.41798 I, y -> 1.03321 + 0.455899 I, z -> -0.301342 - 0.075144 I}, {x -> 2.72138, y -> 1.42833, z -> 0.37027}}}

{0.228611, {{x -> 0.698929 + 3.41798 I, y -> 1.03321 + 0.455899 I, z -> -0.301342 - 0.075144 I}, {x -> 2.72138, y -> 1.42833, z -> 0.37027}}}

The reason that MySqrt is so slow is because MySqrt[y] contains an "unevaluated" If object, while isqrt[y] is a simple arithmetic object:

MySqrt[y]
isqrt[y]

If[Re[y] < 0 && Im[y] >= 0, -Sqrt[y], +Sqrt[y]]

-(-1)^(3/4) Sqrt[I y]

NSolve probably has to use more complicated algebraic methods when presented with these If objects.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.