The Golden Ratio and Point Vortices on the Halfplane¶

This demo is generated by Hanchun Wang base on the work with Prof. Boris Khesin

The paper is available at http://www.math.toronto.edu/khesin/papers/goldenhydroTMIN.pdf

The published version can be viewed at https://link.springer.com/article/10.1007/s00283-021-10099-1

Introduction¶

A single point vortex on the halfplane can be regarded as a vortex dipole that symmetric on the boundary $(y=0)$ in the full plane. Thus an $N$ point vortex system in the halfplane is equivalent to a $2N$ point vortex system in the full plane which are $N$ point vortices and their $N$ images.

The Green's function on the half-plane $\mathbb{R}_+^2$ is \begin{equation} G_{\mathbb{R}^2_+}\left( {z,z'} \right) = - \frac{1}{{2\pi }}\log ||z - z'||+\frac{1}{{2\pi }}\log ||z - z'^*|| \end{equation}

and the Hamiltonian is \begin{equation} H\left( {{z_1},...,{z_k}} \right) = \frac{1}{2}\sum\limits_{i \ne j}^{} {{\Gamma _i}{\Gamma _j}{G_{\mathbb{R}_ + ^2}}\left( {{z_i},{z_j}} \right) + \frac{1}{2}\sum\limits_i^{} {{\Gamma _i^2}{\gamma _{\mathbb{R}_ + ^2}}\left( {{z_i}} \right)} } \end{equation} The first term represents the interaction between different vortices with their images; and the second term with ${{\gamma _{R_ + ^2}}\left( {{z_i}} \right)} = \frac{1}{{2\pi }}\log ||z - z^*||$ represents the interaction between a vortex and its image.

For the two point vortices system on the half plane, the Hamiltonian is $$H= \frac{1}{{4\pi }}\log \left( {{{\left( {2{y_1}} \right)}^{\Gamma _1^2}}{{\left( {2{y_2}} \right)}^{\Gamma _2^2}}{{\left( {\frac{{{{\left( {{x_1} - {x_2}} \right)}^2} + {{\left( {{y_1} + {y_2}} \right)}^2}}}{{{{\left( {{x_1} - {x_2}} \right)}^2} + {{\left( {{y_1} - {y_2}} \right)}^2}}}} \right)}^{{\Gamma _1}{\Gamma _2}}}} \right)$$

The motion of equations are $${{\dot x}_1} = \frac{{{\Gamma _1}}}{{4\pi }}\frac{1}{{{y_1}}} + \frac{{{\Gamma _2}}}{{4\pi }}\left( {\frac{{2\left( {{y_1} + {y_2}} \right)}}{{{{\left( {{x_1} - {x_2}} \right)}^2} + {{\left( {{y_1} + {y_2}} \right)}^2}}} - \frac{{2\left( {{y_1} - {y_2}} \right)}}{{{{\left( {{x_1} - {x_2}} \right)}^2} + {{\left( {{y_1} - {y_2}} \right)}^2}}}} \right)$$ and $${{\dot y}_1} = \frac{{ - {\Gamma _2}}}{{4\pi }}\left( {\frac{{2\left( {{x_1} - {x_2}} \right)}}{{{{\left( {{x_1} - {x_2}} \right)}^2} + {{\left( {{y_1} + {y_2}} \right)}^2}}} - \frac{{2\left( {{x_1} - {x_2}} \right)}}{{{{\left( {{x_1} - {x_2}} \right)}^2} + {{\left( {{y_1} - {y_2}} \right)}^2}}}} \right)$$

To study vortex bifurcations we normalize their strengths by setting $\Gamma_1=\Gamma_2=1$ for a vortex pair and $\Gamma_1=-\Gamma_2=1$ for a dipole. Furthermore, introduce the following dimensionless parameter

$$W:=(P / \Gamma)^2 \exp \left(-4 \pi \mathcal{H} / \Gamma^2\right)$$

measuring the vortex interaction where $\Gamma:=\Gamma_1=\pm \Gamma_2$. As we will see below, the increase of $W$ corresponds to the weakening of the interaction between the point vortices.

We show that $$W=\phi, 1, \frac{1}{\phi}$$ are three bifurcation values.

Dependency¶

In [2]:
import numpy as np
import matplotlib.pyplot as plt

from numpy.testing import assert_almost_equal
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import solve_ivp
from matplotlib.animation import FuncAnimation, PillowWriter

plt.style.use('seaborn-poster')
%matplotlib inline

Utils¶

In [105]:
def calculate_distance_sq(x1, y1, x2, y2):
    return (x1 - x2) ** 2 + (y1 - y2) ** 2

def calculate_energy(x1, y1, x2, y2, s1, s2):
    t1 = np.square(s1) * np.log(2 * y1)
    t2 = np.square(s2) * np.log(2 * y2)
    t3 = np.divide((x1 - x2) ** 2 + (y1 + y2) ** 2, (x1 - x2) ** 2 + (y1 - y2) ** 2)
    t4 = s1 * s2 * np.log(t3)
    return 1 / (4 * np.pi) * (t1 + t2 + t4)

def calculate_p(pos):
    y_coord = pos[:, 1]
    s_lst = pos[:, 2]
    return y_coord @ s_lst

def calculate_lambd_vec(pos):
    x1, y1, s1 = pos[0]
    x2, y2, s2 = pos[1]
    return s2/s1

def calculate_c_vec(pos):
    x1, y1, s1 = pos[0]
    x2, y2, s2 = pos[1]
    e = calculate_energy(x1, y1, x2, y2, s1, s2)
    return np.exp(np.divide(-1 * 4 * np.pi * e, s1 ** 2))

def calculate_w(pos):
    c = calculate_c_vec(pos)
    p = calculate_p(pos)
    lambd = calculate_lambd_vec(pos)
    s = pos[0][2]
    w = np.power(p/s, 1+lambd*lambd)*c
    return w
In [77]:
def rhs(t, S):
    s1,s2 = S[0], S[3]
    x1,x2 = S[1], S[4]
    y1,y2 = S[2], S[5]
    ds1, ds2 = 0,0
    dx1 = -1 * s2 * np.divide((y1 - y2), calculate_distance_sq(x1, y1, x2, y2)) \
        + (-1 * (-1 * s1) * np.divide(1, 2 * y1)) \
        + (-1 * (-1 * s2) * np.divide((y1 + y2), calculate_distance_sq(x1, y1, x2, -1 * y2)))

    dy1 = 1 * s2 * np.divide((x1 - x2), calculate_distance_sq(x1, y1, x2, y2)) \
            + ((-1 * s2) * np.divide((x1 - x2), calculate_distance_sq(x1, y1, x2, -1 * y2)))

    dx2 = -1 * s1 * np.divide((y2 - y1), calculate_distance_sq(x1, y1, x2, y2)) \
            + (-1 * (-1 * s2) * np.divide(1, 2 * y2)) \
            + (-1 * (-1 * s1) * np.divide((y2 + y1), calculate_distance_sq(x1, -1 * y1, x2, y2)))

    dy2 = 1 * s1 * np.divide((x2 - x1), calculate_distance_sq(x1, y1, x2, y2)) \
            + ((-1 * s1) * np.divide((x2 - x1), calculate_distance_sq(x1, -1 * y1, x2, y2)))
    return np.array([ds1, dx1, dy1,ds2,dx2,dy2])

Experiments¶

In [363]:
xmin=-1
xmax=1
ymax=1.5
Frames = 100
skip = 1200
def main(S0,time=4, fig=True, gif=False):
    w= calculate_w(np.array([[S0[1],S0[2],S0[0]],[S0[4],S0[5],S0[3]]]))
    print(f"W={w}")
    w = round(w,3)
    t_eval = np.arange(0, time, 0.0001)
    sol = solve_ivp(rhs, [0, time], S0, t_eval=t_eval, method='Radau')
    if fig:
        plt.figure(figsize = (12, 5))
        plt.plot(sol.y[1], sol.y[2], c='coral')
        plt.plot(sol.y[4], sol.y[5], c='teal')
        plt.text(xmin+(xmax-xmin)*0.75, ymax*0.8, f'W = {w}', fontsize = 26, style='italic')
        plt.ylim(bottom=0)
        plt.grid(alpha=0.5)
        plt.show()
    def animate(n):
        fig.clear()
        N = 8
        n = n*int((len(sol.y[1])-skip*N)/Frames)
        M=len(sol.y[1][n:n+skip*N:skip])
        size = np.linspace(2,10,M)**2
        # print(n)
        plt.scatter(sol.y[1][n:n+skip*N:skip], sol.y[2][n:n+skip*N:skip], s=size, color='coral',edgecolors='grey')
        plt.scatter(sol.y[4][n:n+skip*N:skip], sol.y[5][n:n+skip*N:skip], s=size, color='teal', edgecolors='grey')
        plt.grid(alpha=0.5)
        plt.xlim(xmin, xmax)
        plt.ylim(0,ymax)
        plt.text(xmin+(xmax-xmin)*0.75, ymax*0.8, f'W = {w}', fontsize = 26, style='italic')
    if gif:
        fig = plt.figure(figsize=(12,5))
        ani = FuncAnimation(fig, animate, interval=200, frames=Frames)
        ani.save(f"{w}.gif", dpi=200, writer=PillowWriter(fps=30))
In [ ]:
S0=[ 1.        , -1.70651513,  0.16676086, -1.        ,  1.41108151,    1.14676086]
main(S0, time=5, fig=1, gif=0)

In [ ]:
S0=[ 1.        , -2.30651513,  0.14658999, -1.        ,  0.81108151,    1.16676]
main(S0, time=3.5, fig=1, gif=0)

In [ ]:
S0=[ 1.        , -2.30651513,  0.12658999, -1.        ,  0.81108151,    1.16676]
main(S0, time=3, fig=1, gif=0)

In [ ]:
xmin, xmax, ymax= -15, 7, 8
S0=[ 1.        , -1.30651513,  0.21658999, -1.        ,  5.81108151,    1.16676]
main(S0, time=25, fig=1, gif=0)

Negative dipole¶

In [ ]:
xmin=-2
xmax=8
ymax=1.4
S0=[ 1.        , -1,  0.4545, 1.        ,  -1,    1.06676]
Frames = 120
_=main(S0, time=7, fig=1, gif=0)

In [ ]:
xmin=-2
xmax=14
ymax=1.4
S0=[ 1.        , -1.30651513,  0.27545, 1.        ,  -1.31108151,    1.16676]
main(S0, time=12, fig=1, gif=0)

In [ ]:
xmin, xmax, ymax= -2, 35, 1.4
S0=[ 1.        , -1.30651513,  0.2145, 1.        ,  -1.31108151,    1.16676]
main(S0, time=40, fig=1, gif=0)

In [ ]:
xmin, xmax, ymax= -2, 180, 1.4
skip=10000
S0=[ 1.        , -1.30651513,  0.1945, 1.        ,  -1.31108151,    1.16676]
main(S0, time=200, fig=1, gif=0)