Computing gradient for conformal map between unit disk and unit square

Question

I am currently exploring a Schwarz-Christoffel transformation between the unit square and the unit disk and its inverse. I am following the publication of Fong (2015) which provides the equations $w=G(z)$ from the unit square ($z=x+iy$) to the unit disk ($w=u+iv$), and $z=G^{-1}(w)$ from the unit disk to the unit square:

$$G(z)=\frac{1-i}{\sqrt{2}}cn\left( K_e \frac{1+i}{2}z-K_e,\frac{1}{\sqrt{2}} \right)$$

$$G^{-1}(w)=\frac{1-i}{-K_e}F\left(cos^{-1} \left(\frac{1+i}{\sqrt{2}},\frac{1}{\sqrt{2}} \right) \right)+1-i$$

where $cn$ is a Jacobi elliptic function, $F$ is the incomplete Legendre ellipctic integral of the 1st kind, and $K_e \approx1.854$. I can use the map $G(z)$ to map points on the unit square to the unit disk:

However, a strange thing happens when I want to get the derivative of the transformation from the unit disk to the unit square,

$$\frac{\partial z}{\partial w}=\frac{\partial G^{-1}(w)}{\partial w}=\frac{\sqrt{2}}{K_e\sqrt{1-iw^2}\sqrt{1-\frac{1-iw^2}{\sqrt{2}}}}$$ I gott this derivative from Wolfram Alpha with the following entry:

Derivative of (1-i)/(-K)*EllipticF(acos((1+i)/sqrt(2)*w),1/sqrt(2))+1-i wrt w

It is pretty clear what the vector field should look like: it should follow the color gradient I have plotted above. Unfortunately, when I create a vector field / quiver plot of (say) the real component of the gradient, I encounter some artifacts along the axis from southwest to northeast - the rest of the vector field looks correct. (Please disregard the long, messy arrows in the northwestern and southwestern corners: these result from the singularities in the corners)

Did I make an error somewhere? I have attached a minimal Python code example below

import numpy as np
import matplotlib.pyplot as plt
from mpmath import ellipfun, ellipf

# Create a mesh on the unit square
X,Y = np.meshgrid(
    np.linspace(-1,1,51),
    np.linspace(-1,1,51))
X   = np.ndarray.flatten(X)
Y   = np.ndarray.flatten(Y)
XY = np.column_stack((X,Y))

# Convert the points to complex numbers
z   = XY[:,0] + 1j*XY[:,1]

# That's an approximate value
Ke  = 1.854

# Schwarz-Christoffel transformation from unit square to unit disk
cn = ellipfun('cn')
w   = [(1-1j)/np.sqrt(2)*cn(
        u   = Ke*(1+1j)/2*i-Ke,
        k   = 1/np.sqrt(2)) for i in z]
w   = np.asarray([np.real(i) + 1j*np.imag(i) for i in w],dtype=np.complex)

# Plot the current state
plt.figure()
plt.subplot(1,2,1)
plt.title('unit square')
plt.scatter(np.real(z),np.imag(z),c=np.real(z),s=0.5)
plt.axis('equal')
plt.xlabel('x')
plt.ylabel('y')

plt.subplot(1,2,2)
plt.title('unit disk')
plt.scatter(np.real(w),np.imag(w),c=np.real(z),s=0.5)
plt.axis('equal')
plt.xlabel('u')
plt.ylabel('v')


plt.figure()
# Calculate gradient from unit disk to unit square, using the inverse transformation
grad   = np.sqrt(2)/(1.854*np.sqrt(1-1j*w**2)*np.sqrt(1-(1-1j*w**2)/np.sqrt(2)))

dudx = np.real(grad)
dudy = -np.imag(grad)

plt.quiver(
    np.real(w),
    np.imag(w),
    dudx,
    dudy,
    scale = 100)

plt.title('derivative $\partial z / \partial w$')
plt.xlabel('u')
plt.ylabel('v')

plt.axis('equal')

Answer 1

hope I'm not too late. Thank you for this interesting question/topic.

I think I found a couple of problems with your inverse formula. It should be as follows:

$G^{-1}(w)= \frac{1-i}{-K_e} F\left(\text{cos}^{-1}\left(\frac{1+i}{\sqrt 2} w\right), \frac12\right) + 1-i$

This depends on how the function $F$ is defined (if is second argument should take the square root or not basically). I hope this goes in the same direction of the intentions of the authors of you article. If you input the above formula into wolfram alpha, you'll get:

$\frac{\partial G^{-1}}{\partial w}(w) = \frac{2}{K_e\sqrt{w^4+1}}$

this looks correct to me.

So, with this definition of the inverse function you can check if it is correct by trying to get to the square from the circle (which you could not do using the formula you provided):

# Does the inverse: given a circle, transforms into a square
# should be equal to z
zz = [
        ((1-1j)/(-Ke))*
        ellipf(
            np.arccos((1+1j)/np.sqrt(2)*u), 1/2
        ) for u in w
    ]
zz  = np.asarray([np.real(i)+1 + 1j*(np.imag(i)-1) for i in zz], dtype=np.complex128)

# this line would fail with the inverse definition that you'd used
assert np.allclose(zz, z), "Inverse transform is not working !!!"

# Plot the current state
plt.figure()
plt.subplot(1,2,1)
plt.title('unit disk')
plt.scatter(np.real(w),np.imag(w),c=np.real(z),s=0.5)
plt.axis('equal')
plt.xlabel('u')
plt.ylabel('v')

plt.subplot(1,2,2)
plt.title('unit square')
plt.scatter(np.real(zz),np.imag(zz), c=np.real(z),s=0.5)
plt.axis('equal')
plt.xlabel('x')
plt.ylabel('y')

plt.show()

And then you can finally test the formula of the gradient:

# Calculate gradient from unit disk to unit square, using the inverse transformation
grad   = 2/(1.854*np.sqrt(w**4+1))

dudx = np.real(grad)
dudy = -np.imag(grad)

plt.figure()
plt.quiver(np.real(w), np.imag(w), dudx, dudy, scale=100)

plt.title('derivative $\partial z / \partial w$')
plt.xlabel('u')
plt.ylabel('v')

plt.axis('equal')
plt.show()

I would post the vector field but I'm a newbie here, so I don't have the rights. But you can plot it, it looks good !

Hope this helps :)