1
$\begingroup$

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:

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:

good result contour

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

Good result

Bad results

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

Bad result

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")
```
$\endgroup$
5
  • $\begingroup$ 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? $\endgroup$ Jul 16 at 21:19
  • $\begingroup$ The solution in the annulus is correct. The problem is definitely in the transformations T and Ti because of z in the denominator. $\endgroup$ Jul 19 at 12:42
  • $\begingroup$ 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? $\endgroup$ Jul 19 at 14:25
  • $\begingroup$ 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 . $\endgroup$ Jul 19 at 18:13
  • $\begingroup$ 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. $\endgroup$ Jul 19 at 21:52
2
$\begingroup$

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)$.

$\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.