# Finding the source of numerical instability in a electrostatic problem solved by conformal mapping

I'm using conformal mapping to solve a 2D electrostatic problem (calculating the potential $$u(x,y)$$ in the plane). Let $$C_1$$ and $$C_2$$ be two circles at an electric potential $$U_1$$ and $$U_2$$, respectively. They have radii $$R_1$$ and $$R_2$$, respectively, and their centers are separated by a distance $$h$$.

I noticed that, for some values of $$R_1$$, $$R_2$$, $$U_1$$, $$U_2$$ and $$h$$, the potential $$u(x,y)$$ explodes in certain points. I'm led to believe that this is due to numerical instability and I would like to pinpoint where is that instability coming from and how can I attenuate it (maybe by manipulating the formulas). More specifically: how do I perform a sensitivity analysis of a boundary value problem solved by conformal mapping?

# Algorithm to find the solution

My steps to solve that problems are as follows:

• find a conformal mapping $$T$$ that makes the two circles $$C_1$$ and $$C_2$$ concentric;
• solve the potential problem in the annulus (has a closed formula solution $$v(r)$$), see for example: Conformal Mapping and Bipolar Coordinate for Eccentric Laplace Problems;
• map the potential $$v(r)$$ in the annulus back into the original plane.

The conformal mapping $$T(z)$$ can be found by chaining Möbius Transformations, where $$z$$ is a point in the complex plane. The parameters $$p$$, $$\xi_1$$ and $$\xi_2$$ are found numerically (I'll omit the details to find them).

$$T(z) = \frac{R_1 (p + \xi_1 + 1) - z(p + \xi_1 - 1)}{R_1 (p + \xi_2 + 1) - z(p + \xi_2 - 1)}$$

The inverse transform $$T^{-1}(z)$$ is given by:

$$T^{-1}(z) = R_1 \frac{(p + \xi_1 + 1) - z (p + \xi_2 + 1)}{(p + \xi_1 - 1) - z (p + \xi_2 - 1)}$$

The potential in the annulus is given by, where $$\log \equiv \ln$$, $$r_1 = |T(R_1)|$$ and $$r_2 = |T(R_2 + h)|$$:

$$v(r) = T(U_1) + \frac{[T(U_1) - T(U_2)] \log(r/r_1)}{\log(r_1 / r_2)}$$

Finally, the potential in the original plane is given by:

$$u(z) = T^{-1}(v(\,|T(z)|\,))$$

# Example Results

See julia script further down below.

## Good results

By using the following parameters, I get a good result that looks like what I expect from this problem.

$$R_1 = 2.2 \\ R_2 = 1.1 \\ U_1 = 1/2 \\ U_2 = -1/2 \\ h = 8$$

Contour plot:

Potential $$u$$ along a straight line connecting the circles:

## Bad results

By using the previous parameters but $$U_1 = 12$$, the potential $$u$$ along a straight line connecting the circles is the following:

See that, by just changing $$U_1$$, the result explodes in some points and we get, for example, $$u(5.67, 0) \approx -18935.11$$.

## Julia script

using NLsolve
using Polynomials
using Plots
gr()

const ϵ0 = 8.8541878176e-12  # free space electric permittivity

function find_trans(r1, r2, h, n=10)
"""
Finds the transformation that maps the area outside two circles of radii
R1 and R2, respectively, into an annulus. The first circle is centered
at (0,0) and the second on at (h,0). It is assumed that R1 ≥ R2.

Returns
-------
T : the direct transform
Ti : the inverse transform
params : tuple of the parameters (p, ξ1, ξ2)
"""
if r1 < r2
throw(ArgumentError("r1 < r2"))
end
if h < r1 + r2
throw(ArgumentError("h < r1 + r2"))
end
θ = range(0, 2π, length=n)
z1 = @. r1*(cos(θ) + 1im*sin(θ))
z2 = @. r2*(cos(θ) + 1im*sin(θ)) + h
# bigger circle C1 into a line over the imaginary axis
f(z) = (r1 + z)/(r1 - z)
# centers the transformed circle C2 into the origin
w1 = f(z2[1])
w2 = f(z2[cld(n,2)])
function f!(F, x)
F[1] = abs(w1)^2 - 2*real(w1*conj(x[1])) + abs(x[1])^2 - x[2]^2
F[2] = abs(w2)^2 - 2*real(w2*conj(x[1])) + abs(x[1])^2 - x[2]^2
end
# sol.zero: center of the transformed circle and its radius
sol1 = nlsolve(f!, [0.0, 2r1])
p = sol1.zero[1]
r = sol1.zero[2]
# transforms the line into a circle centered at the origin,
# while preserving the center of the other circle
# x1 and x2 are points symmetric to both circle and line
x1, x2 = roots(Polynomial([r^2, 2p, 1]))
# Resulting Möbius Transform
T(z) = (r1*(p + x1 + 1) - z*(p + x1 - 1))/(r1*(p + x2 + 1) - z*(p + x2 - 1))
Ti(z) = r1*((p + x1 + 1) - z*(p + x2 + 1))/((p + x1 - 1) - z*(p + x2 - 1))
return T, Ti, (p, x1, x2)
end

function find_potential(R1, R2, h, U1, U2, n=100)
"""
Returns the potential function u and the points on the
boundaries zc1 and zc2.
"""
θ = range(0, 2π, length=n)
zc1 = @. R1*(cos(θ) + 1im*sin(θ))
zc2 = @. R2*(cos(θ) + 1im*sin(θ)) + h
T, Ti, params = find_trans(R1, R2, h)

r1 = abs(T(R1))
r2 = abs(T(R2 + h))
v1 = (T(U1))
v2 = (T(U2))
# potential in the annulus:
v(w) = v1 + (v1 - v2) * log(abs(w)/r1) / log(r1/r2)
# potential in the original plane:
u(z) = Ti( v( abs(T(z)) ) )
return u, zc1, zc2
end

## Good result
R1 = 2.2  # conductor 1 radius
R2 = 1.1  # conductor 2 radius
h = 8  # speration between conductors' centers
U1 = 1/2  # potential of the conductor 1
U2 = -1/2  # potential of the conductor 2
@assert R1 >= R2 && h > R1 + R2 && U1 > U2  # sanity test
n = 100  # "resolution" of the plots
x = range(R1 * 1.1, h - R2 * 1.1, length=n)
u, zc1, zc2 = find_potential(R1, R2, h, U1, U2, n)
p1 = plot(x, u.(x), xlabel="x", ylabel="u(x,0)", label="")
savefig("goodpot.png")

s = 1.5
x = range(-s*R1, h + s*R2, length=n)
y = range(-2s*R1, 2s*R1, length=n)
X = repeat(reshape(x, 1, :), length(y), 1)
Y = repeat(y, 1, length(x))
function umasked(x, y, mask=true)
z = x + 1im*y
pot = u(z)
if mask
if pot > U1
pot = U1
elseif pot < U2
pot = U2
end
#=if abs(z) <= R1 || abs(z - h) <= R2
pot = nothing
end=#
end
return pot
end
U = map(umasked, X, Y);

contour(x, y, U, fill=true, color=:haline, aspect_ratio=1)
plot!(zc1, color=:red, label="")
plot!(zc2, color=:red, label="")
savefig("goodpot_contour.png")

## Bad result
R1 = 2.2  # conductor 1 radius
R2 = 1.1  # conductor 2 radius
h = 8  # speration between conductors' centers
U1 = 12  # potential of the conductor 1
U2 = -1/2  # potential of the conductor 2
@assert R1 >= R2 && h > R1 + R2 && U1 > U2  # sanity test
n = 100  # "resolution" of the plots
x = range(R1 * 1.1, h - R2 * 1.1, length=n)
u, zc1, zc2 = find_potential(R1, R2, h, U1, U2, n)
plot(x, u.(x), xlabel="x", ylabel="u(x,0)", label="")
savefig("badpot.png")

s = 1.5
x = range(-s*R1, h + s*R2, length=n)
y = range(-2s*R1, 2s*R1, length=n)
X = repeat(reshape(x, 1, :), length(y), 1)
Y = repeat(y, 1, length(x))
function umasked(x, y, mask=true)
z = x + 1im*y
pot = u(z)
if mask
if pot > U1
pot = U1
elseif pot < U2
pot = U2
end
#=if abs(z) <= R1 || abs(z - h) <= R2
pot = nothing
end=#
end
return pot
end
U = map(umasked, X, Y);

contour(x, y, U, fill=true, color=:haline, aspect_ratio=1)
plot!(zc1, color=:red, label="")
plot!(zc2, color=:red, label="")
savefig("badpot_contour.png")
$$$$
`
• The solution in the annulus does not have any problem for U1=12? Only the solution transformed to the two-circle domain has the problem? – Maxim Umansky Jul 16 at 21:19
• The solution in the annulus is correct. The problem is definitely in the transformations T and Ti because of z in the denominator. – Pedro H. N. Vieira Jul 19 at 12:42
• But the transformation should be exactly the same for U1=12 and U1=1/2, isn't it? Can you just save the transformation (not recalculating it) from the good case, and use it again for the case with U1=12? – Maxim Umansky Jul 19 at 14:25
• The transformations do not depend on the contour values U1 or U2. It depends only on the geometry of the problem and what intermediate steps I take to find the mapping (parameters xi and p). For this example, T_i blows near -3.56 while T(12) = -9.7, T(-0.5) = -0.48 and T(0.5) = -0.06. The solution is to tweak parameters xi and p. That's why I want to make a sensitivity analysis: what region is stable for given parameters . – Pedro H. N. Vieira Jul 19 at 18:13
• I would suggest to use one good solution (for a given geometry) to obtain all other solutions for this geometry, using the linear properties of the Laplace equation that you are solving. – Maxim Umansky Jul 19 at 21:52

## 1 Answer

The cause of the numerical glitch has probably something to do with poor convergence of some series involved in the conformal map calculation for a particular location. However, here is how to use one good solution for this problem, for a particular geometry, to produce any other one.

Suppose we have a good solution $$\phi_0(x,y)$$ corresponding to $$U_1$$=0 and $$U_2$$=1. Then, from the linear properties of the Laplace equation, the solution for the same geometry and some other choice of $$U_1$$ and $$U_2$$ is given by

$$\phi(x,y) = U_1 + (U_2-U_1) \phi_0(x,y)$$.