Regularity conditions for tensor components in polar coordinates¶
In [3]:
%load_ext autoreload
%autoreload 2
import numpy as np
import sympy as sp
import matplotlib.pyplot as plt
import plotting_utils as plottings
from IPython.display import display
The autoreload extension is already loaded. To reload it, use: %reload_ext autoreload
In [4]:
plt.rcParams['mathtext.fontset'] = 'cm'
Utilities for symbolic manipulation
In [5]:
s = sp.Symbol("s", positive=True)
p = sp.Symbol(r"\phi")
x = sp.Symbol("x")
y = sp.Symbol("y")
A = sp.Symbol("A")
A_x = sp.Symbol("A_x")
A_y = sp.Symbol("A_y")
A_s = sp.Symbol("A_s")
A_p = sp.Symbol(r"A_{\phi}")
A_xx = sp.Symbol("A_{xx}")
A_xy = sp.Symbol("A_{xy}")
A_yx = sp.Symbol("A_{yx}")
A_yy = sp.Symbol("A_{yy}")
A_ss = sp.Symbol("A_{ss}")
A_sp = sp.Symbol(r"A_{s\phi}")
A_ps = sp.Symbol(r"A_{\phi s}")
A_pp = sp.Symbol(r"A_{\phi\phi}")
cart2polar_map = {
x: s*sp.cos(p),
y: s*sp.sin(p)
}
polar2cart_map = {
s: sp.sqrt(x**2 + y**2),
sp.exp(+sp.I*p): (x + sp.I*y)/sp.sqrt(x**2 + y**2),
sp.exp(-sp.I*p): (x - sp.I*y)/sp.sqrt(x**2 + y**2)
}
In [6]:
R = sp.Matrix([[sp.cos(p), sp.sin(p)], [-sp.sin(p), sp.cos(p)]])
R
Out[6]:
$\displaystyle \left[\begin{matrix}\cos{\left(\phi \right)} & \sin{\left(\phi \right)}\\- \sin{\left(\phi \right)} & \cos{\left(\phi \right)}\end{matrix}\right]$
In [7]:
U = sp.Matrix([[1, 1], [-sp.I, sp.I]])
U
Out[7]:
$\displaystyle \left[\begin{matrix}1 & 1\\- i & i\end{matrix}\right]$
In [8]:
R @ U
Out[8]:
$\displaystyle \left[\begin{matrix}- i \sin{\left(\phi \right)} + \cos{\left(\phi \right)} & i \sin{\left(\phi \right)} + \cos{\left(\phi \right)}\\- \sin{\left(\phi \right)} - i \cos{\left(\phi \right)} & - \sin{\left(\phi \right)} + i \cos{\left(\phi \right)}\end{matrix}\right]$
In [9]:
np.linalg.svd(np.array([[1, 1j], [1, -1j]]))
Out[9]:
(array([[-0.70710678+0.j, -0.70710678+0.j],
[-0.70710678+0.j, 0.70710678+0.j]]),
array([1.41421356, 1.41421356]),
array([[-1.+0.j, 0.+0.j],
[ 0.+0.j, 0.-1.j]]))
In [24]:
A_val = s**3*sp.cos(p)
# A_val = s*sp.exp(sp.I*3*p)
display(sp.Eq(A, A_val))
display(sp.Eq(A, A_val.rewrite(sp.exp)))
$\displaystyle A = s^{3} \cos{\left(\phi \right)}$
$\displaystyle A = s^{3} \left(\frac{e^{i \phi}}{2} + \frac{e^{- i \phi}}{2}\right)$
does not have the leading order behaviour $s^1$ required.
As we convert it to Cartesian coordinates, it takes the following form in Cartesian coordinates
In [25]:
A_cart = A_val.rewrite(sp.exp).subs(polar2cart_map).simplify()
sp.Eq(A, A_cart)
Out[25]:
$\displaystyle A = x \left(x^{2} + y^{2}\right)$
In [9]:
sp.diff(A_cart, (y, 2)).simplify().factor()
Out[9]:
$\displaystyle - \frac{x \left(x^{2} - 2 y^{2}\right)}{\left(x^{2} + y^{2}\right)^{\frac{5}{2}}}$
The field is thus singular, even discontinous at the origin.
The visualization of the field is as follows
In [44]:
A_numerical = sp.lambdify((s, p), A_val)
fig, ax = plottings.polar_singularity_scalar(sfunc=A_numerical, shape='polar', num=160)
ax.set_title(r'$f(r, \phi)$', fontsize=18)
ax.set_aspect('equal', 'box')
# ax.set_title(r'$\widetilde{f}(x,y) = f(r(x,y), \phi(x,y))$', fontsize=16)
save_name = '../out/imgs/Singularity_scalar__cosp__sp'
# save_name = None
if save_name is not None:
plt.savefig(save_name + '.png', dpi=150, bbox_inches='tight')
plt.savefig(save_name + '.pdf', bbox_inches='tight')
plt.show()
In [101]:
A_fields = [sp.lambdify((s, p), A_val), ]
fig, axes = plt.subplots(nrows=1, ncols=len(A_fields), figsize=(5.5, 5))
axes = np.atleast_1d(axes)
X_mesh, Y_mesh, S_mesh, P_mesh = plottings.polar_mesh(shape='Cartesian', num=200)
A_field_vals = [A_tmp(S_mesh, P_mesh) for A_tmp in A_fields]
vmin = min([np.real(A_tmp_val).min() for A_tmp_val in A_field_vals])
vmax = max([np.real(A_tmp_val).max() for A_tmp_val in A_field_vals])
plot_kw = {'shading': 'gouraud', 'vmin': vmin, 'vmax': vmax, 'cmap': 'inferno'}
for i_f, ax in enumerate(axes):
im = ax.pcolormesh(X_mesh, Y_mesh, np.real(A_field_vals[i_f]), **plot_kw)
ax.set_aspect('equal', 'box')
ax.set_title(r'$\widetilde{f}(x,y) = f(r(x,y), \phi(x,y))$', fontsize=18)
ax.set_xlabel('$x$', fontsize=16)
ax.set_ylabel('$y$', fontsize=16)
fig.subplots_adjust(right=0.85)
cbar_ax = fig.add_axes([0.9, 0.15, 0.04, 0.7])
fig.colorbar(im, cax=cbar_ax)
save_name = '../out/imgs/Singularity_scalar__cosp__xy_vertbar'
# save_name = None
if save_name is not None:
plt.savefig(save_name + '.png', dpi=150, bbox_inches='tight')
plt.savefig(save_name + '.pdf', bbox_inches='tight')
plt.show()
Example: confirmation of the conditions¶
Now we simplify modify the field by adding a $s^1$ prefactor
In [49]:
A_val = s*sp.cos(p)
display(sp.Eq(A, A_val))
display(sp.Eq(A, A_val.rewrite(sp.exp).expand()))
$\displaystyle A = s \cos{\left(\phi \right)}$
$\displaystyle A = \frac{s e^{i \phi}}{2} + \frac{s e^{- i \phi}}{2}$
As we convert it to Cartesian coordinates, it takes the following form in Cartesian coordinates
In [7]:
sp.Eq(A, A_val.rewrite(sp.exp).subs(polar2cart_map).simplify())
Out[7]:
$\displaystyle A = x$
The field is thus regular:
In [8]:
A_numerical = sp.lambdify((s, p), A_val)
fig, ax = plottings.polar_singularity_scalar(sfunc=A_numerical)
plt.show()
Vector¶
The regularity condition for vector (Lewis and Bellan 1990): $$ \begin{aligned} A_s &= sp(s^2) + \sum_{m\neq 0} \left(\lambda_m s^{|m|-1} + s^{|m|+1}p(s^2)\right) e^{im\phi},\\ A_\phi &= sp(s^2) + \sum_{m\neq 0} \left(i\mathrm{sgn}(m)\lambda_m s^{|m|-1} + s^{|m|+1}p(s^2)\right) e^{im\phi} \end{aligned} $$
In [85]:
def rotate_polar2cart_vec(F_s, F_p):
return F_s*sp.cos(p) - F_p*sp.sin(p), F_s*sp.sin(p) + F_p*sp.cos(p)
def rotate_cart2polar_vec(F_x, F_y):
return F_x*sp.cos(p) + F_y*sp.sin(p), -F_x*sp.sin(p) + F_y*sp.cos(p)
(Counter)Example: Lewis-Bellan relations are not necessary for one component to be regular¶
Consider the following vector field
In [10]:
A_s_val = (1 - sp.cos(2*p))/s
A_p_val = sp.sin(2*p)/s
A_s_exp = A_s_val.rewrite(sp.exp).expand()
A_p_exp = A_p_val.rewrite(sp.exp).expand()
display(sp.Eq(A_s, A_s_val))
display(sp.Eq(A_p, A_p_val))
display(sp.Eq(A_s, A_s_exp))
display(sp.Eq(A_p, A_p_exp))
$\displaystyle A_{s} = \frac{1 - \cos{\left(2 \phi \right)}}{s}$
$\displaystyle A_{\phi} = \frac{\sin{\left(2 \phi \right)}}{s}$
$\displaystyle A_{s} = - \frac{e^{2 i \phi}}{2 s} + \frac{1}{s} - \frac{e^{- 2 i \phi}}{2 s}$
$\displaystyle A_{\phi} = - \frac{i e^{2 i \phi}}{2 s} + \frac{i e^{- 2 i \phi}}{2 s}$
which does not fulfill the Lewis-Bellan relations at all (even the Fourier coefficients are singular!).
Nevertheless, as we convert it to Cartesian coordinates, it takes the following form in Cartesian coordinates
In [11]:
A_x_val, A_y_val = rotate_polar2cart_vec(A_s_exp, A_p_exp)
display(sp.Eq(A_x, A_x_val.rewrite(sp.exp).expand().subs(polar2cart_map).simplify()))
display(sp.Eq(A_y, A_y_val.rewrite(sp.exp).expand().subs(polar2cart_map).simplify()))
$\displaystyle A_{x} = 0$
$\displaystyle A_{y} = \frac{2 y}{x^{2} + y^{2}}$
The field is of course singular, even blows up at the origin for $A_y$, but the $A_x$ component is however perfectly regular. Therefore, Lewis and Bellan shouldn't have jumped to their conclusion by deriving only the $A_x$ component (eq. 14 in their paper).
In [12]:
A_s_numerical = sp.lambdify((s, p), A_s_val)
A_p_numerical = sp.lambdify((s, p), A_p_val)
fig, ax = plottings.polar_singularity_vector(vfunc_s=A_s_numerical, vfunc_p=A_p_numerical)
plt.show()
Example: a singular vector field showing necessity of the conditions¶
We look at a field whose Fourier coefficients have the correct prefactor in $s$, but have the wrong configuration of coefficients
In [104]:
A_s_val = (1 + s**2)*sp.cos(p)
A_p_val = (-1 - s**2)*sp.sin(p)
A_s_exp = A_s_val.rewrite(sp.exp).expand()
A_p_exp = A_p_val.rewrite(sp.exp).expand()
display(sp.Eq(A_s, A_s_val))
display(sp.Eq(A_p, A_p_val))
display(sp.Eq(A_s, A_s_exp))
display(sp.Eq(A_p, A_p_exp))
$\displaystyle A_{s} = \left(s^{2} + 1\right) \cos{\left(\phi \right)}$
$\displaystyle A_{\phi} = \left(- s^{2} - 1\right) \sin{\left(\phi \right)}$
$\displaystyle A_{s} = \frac{s^{2} e^{i \phi}}{2} + \frac{s^{2} e^{- i \phi}}{2} + \frac{e^{i \phi}}{2} + \frac{e^{- i \phi}}{2}$
$\displaystyle A_{\phi} = \frac{i s^{2} e^{i \phi}}{2} - \frac{i s^{2} e^{- i \phi}}{2} + \frac{i e^{i \phi}}{2} - \frac{i e^{- i \phi}}{2}$
As we convert it to Cartesian coordinates, it takes the following form in Cartesian coordinates
In [105]:
A_x_val, A_y_val = rotate_polar2cart_vec(A_s_exp, A_p_exp)
A_x_xy = A_x_val.rewrite(sp.exp).expand().subs(polar2cart_map).simplify()
A_y_xy = A_y_val.rewrite(sp.exp).expand().subs(polar2cart_map).simplify()
display(sp.Eq(A_x, A_x_xy))
display(sp.Eq(A_y, A_y_xy))
$\displaystyle A_{x} = x^{2} + y^{2} + 1$
$\displaystyle A_{y} = 0$
In [106]:
dAx_dx = sp.diff(A_x_xy, x).subs(cart2polar_map).simplify()
dAx_dy = sp.diff(A_x_xy, y).subs(cart2polar_map).simplify()
dAy_dx = sp.diff(A_y_xy, x).subs(cart2polar_map).simplify()
dAy_dy = sp.diff(A_y_xy, y).subs(cart2polar_map).simplify()
display(dAx_dx, dAx_dy, dAy_dx, dAy_dy)
$\displaystyle 2 s \cos{\left(\phi \right)}$
$\displaystyle 2 s \sin{\left(\phi \right)}$
$\displaystyle 0$
$\displaystyle 0$
The field is singular and discontinous at the origin:
In [108]:
A_s_numerical = sp.lambdify((s, p), A_s_val)
A_p_numerical = sp.lambdify((s, p), A_p_val)
fig, axes = plottings.polar_singularity_vector(vfunc_s=A_s_numerical, vfunc_p=A_p_numerical, shape='Cartesian', components='Cartesian', num=200)
axes[0].set_title(r'$V_x$', fontsize=18)
axes[1].set_title(r'$V_y$', fontsize=18)
# axes[0].set_title(r'$V_r(r,\phi)$', fontsize=18)
axes[0].set_aspect('equal', 'box')
# axes[1].set_title(r'$V_\phi(r,\phi)$', fontsize=18)
axes[1].set_aspect('equal', 'box')
save_name = '../out/imgs/Regular_vector__1s^2cps__xy'
# save_name = None
if save_name is not None:
plt.savefig(save_name + '.png', dpi=150, bbox_inches='tight')
plt.savefig(save_name + '.pdf', bbox_inches='tight')
plt.show()
In [94]:
A_fields = [sp.lambdify((s, p), dAx_dx), sp.lambdify((s, p), dAy_dy)]
fig, axes = plt.subplots(nrows=1, ncols=len(A_fields), figsize=(12, 5.2))
X_mesh, Y_mesh, S_mesh, P_mesh = plottings.polar_mesh(shape='Cartesian', num=200)
A_field_vals = [A_tmp(S_mesh, P_mesh) for A_tmp in A_fields]
vmin = min([np.real(A_tmp_val).min() for A_tmp_val in A_field_vals])
vmax = max([np.real(A_tmp_val).max() for A_tmp_val in A_field_vals])
plot_kw = {'shading': 'gouraud', 'vmin': vmin, 'vmax': vmax, 'cmap': 'inferno'}
for i_f, ax in enumerate(axes):
im = ax.pcolormesh(X_mesh, Y_mesh, np.real(A_field_vals[i_f]), **plot_kw)
ax.set_aspect('equal', 'box')
axes[0].set_xlabel('$x$', fontsize=16)
axes[0].set_ylabel('$y$', fontsize=16)
axes[0].set_title(r'$dV_x/dx$', fontsize=18)
axes[1].set_xlabel('$x$', fontsize=16)
axes[1].set_title(r'$dV_x/dy$', fontsize=18)
fig.subplots_adjust(right=0.87)
cbar_ax = fig.add_axes([0.9, 0.15, 0.02, 0.7])
fig.colorbar(im, cax=cbar_ax)
save_name = '../out/imgs/Singularity_vector__scosp-ssinp__derivatives'
# save_name = None
if save_name is not None:
plt.savefig(save_name + '.png', dpi=150, bbox_inches='tight')
plt.savefig(save_name + '.pdf', bbox_inches='tight')
plt.show()
Example: confirmation of the conditions¶
Now we simply change the sign of $A_\phi$ so that the Fourier coefficients are indeed $i\mathrm{sgn}(m)$ of the $A_s$ coefficients at the lowest order
In [16]:
A_s_val = sp.cos(p)
A_p_val = -sp.sin(p)
A_s_exp = A_s_val.rewrite(sp.exp).expand()
A_p_exp = A_p_val.rewrite(sp.exp).expand()
display(sp.Eq(A_s, A_s_val))
display(sp.Eq(A_p, A_p_val))
display(sp.Eq(A_s, A_s_exp))
display(sp.Eq(A_p, A_p_exp))
$\displaystyle A_{s} = \cos{\left(\phi \right)}$
$\displaystyle A_{\phi} = - \sin{\left(\phi \right)}$
$\displaystyle A_{s} = \frac{e^{i \phi}}{2} + \frac{e^{- i \phi}}{2}$
$\displaystyle A_{\phi} = \frac{i e^{i \phi}}{2} - \frac{i e^{- i \phi}}{2}$
The field components are perfectly regular in Cartesian coordinates
In [17]:
A_x_val, A_y_val = rotate_polar2cart_vec(A_s_exp, A_p_exp)
display(sp.Eq(A_x, A_x_val.rewrite(sp.exp).expand().subs(polar2cart_map).simplify()))
display(sp.Eq(A_y, A_y_val.rewrite(sp.exp).expand().subs(polar2cart_map).simplify()))
$\displaystyle A_{x} = 1$
$\displaystyle A_{y} = 0$
In [18]:
A_s_numerical = sp.lambdify((s, p), A_s_val)
A_p_numerical = sp.lambdify((s, p), A_p_val)
fig, ax = plottings.polar_singularity_vector(vfunc_s=A_s_numerical, vfunc_p=A_p_numerical)
plt.show()
Rank-2 tensor¶
Holdenried-Chernoff (2021) derives the following, $$ A_{ij} = \left[A_{ij}^{00} + O(s^2)\right] + \sum_{|m|=1} \left[A_{ij}^{m0} s + O(s^3)\right] e^{im\phi} + \sum_{|m|\geq 2} \left[A_{ij}^{m0} s^{|m|-2} + A_{ij}^{m1} s^{|m|} + O(s^{|m|+2})\right] e^{im\phi} $$ where the following relations apply $$ m = 0: \left\{\begin{aligned} &A_{ss}^{00} = A_{\phi\phi}^{00} \\ &A_{s\phi}^{00} = -A_{\phi s}^{00} \end{aligned}\right.,\qquad |m| \geq 2: \left\{\begin{aligned} &A_{ss}^{m0} = - A_{\phi\phi}^{m0}\\ &A_{s\phi}^{m0} = A_{\phi s}^{m0} \\ &A_{s\phi}^{m0} = i \mathrm{sgn}(m) A_{ss}^{m0} \end{aligned}\right. $$
I show that these relations are not exhaustive, and two other conditions are required to make the conditions sufficient. The complete relations are $$ m = 0: \left\{\begin{aligned} &A_{ss}^{00} = A_{\phi\phi}^{00} \\ &A_{s\phi}^{00} = -A_{\phi s}^{00} \end{aligned}\right.,\qquad |m| = 1: A_{s\phi}^{m0} + A_{\phi s}^{m0} = i\mathrm{sgn}(m) \left(A_{ss}^{m0} - A_{\phi\phi}^{m0}\right),\qquad |m| \geq 2: \left\{\begin{aligned} &A_{ss}^{m0} = - A_{\phi\phi}^{m0}\\ &A_{s\phi}^{m0} = A_{\phi s}^{m0} \\ &A_{s\phi}^{m0} = i \mathrm{sgn}(m) A_{ss}^{m0} \\ &A_{s\phi}^{m1} + A_{\phi s}^{m1} = i\mathrm{sgn}(m)\left(A_{ss}^{m1} - A_{\phi\phi}^{m1}\right) \end{aligned}\right. $$
In [10]:
def rotate_polar2cart_tensor2(F_ss, F_sp, F_ps, F_pp):
F_xx = sp.cos(p)**2*F_ss + sp.sin(p)**2*F_pp - sp.sin(p)*sp.cos(p)*(F_sp + F_ps)
F_yy = sp.sin(p)**2*F_ss + sp.cos(p)**2*F_pp + sp.sin(p)*sp.cos(p)*(F_sp + F_ps)
F_xy = sp.sin(p)*sp.cos(p)*(F_ss - F_pp) + sp.cos(p)**2*F_sp - sp.sin(p)**2*F_ps
F_yx = sp.sin(p)*sp.cos(p)*(F_ss - F_pp) + sp.cos(p)**2*F_ps - sp.sin(p)**2*F_sp
return F_xx, F_xy, F_yx, F_yy
Example: a singular field showing necessity of the |m|=1 condition¶
The old conditions involve no additional requirements at $|m|=1$. To test if the new condition is really necessary, let us look at a field
In [11]:
A_ss_val = s*sp.sin(p)
A_pp_val = sp.S.Zero
A_sp_val = A_ps_val = s*sp.cos(p)
display(sp.Eq(A_ss, A_ss_val))
display(sp.Eq(A_sp, A_sp_val))
display(sp.Eq(A_ps, A_ps_val))
display(sp.Eq(A_pp, A_pp_val))
A_ss_exp = A_ss_val.rewrite(sp.exp).expand()
A_sp_exp = A_sp_val.rewrite(sp.exp).expand()
A_ps_exp = A_ps_val.rewrite(sp.exp).expand()
A_pp_exp = A_pp_val.rewrite(sp.exp).expand()
display(sp.Eq(A_ss, A_ss_exp))
display(sp.Eq(A_sp, A_sp_exp))
display(sp.Eq(A_ps, A_ps_exp))
display(sp.Eq(A_pp, A_pp_exp))
$\displaystyle A_{ss} = s \sin{\left(\phi \right)}$
$\displaystyle A_{s\phi} = s \cos{\left(\phi \right)}$
$\displaystyle A_{\phi s} = s \cos{\left(\phi \right)}$
$\displaystyle A_{\phi\phi} = 0$
$\displaystyle A_{ss} = - \frac{i s e^{i \phi}}{2} + \frac{i s e^{- i \phi}}{2}$
$\displaystyle A_{s\phi} = \frac{s e^{i \phi}}{2} + \frac{s e^{- i \phi}}{2}$
$\displaystyle A_{\phi s} = \frac{s e^{i \phi}}{2} + \frac{s e^{- i \phi}}{2}$
$\displaystyle A_{\phi\phi} = 0$
This field does have all the correct prefactors $s^1$ at $|m|=1$, but $A_{s\phi}^{m0} + A_{\phi s}^{m0} \neq i \mathrm{sgn}(m) (A_{ss}^{m0} - A_{\phi\phi}^{m0})$. As we convert it to Cartesian coordinates, it takes the following form in Cartesian coordinates
In [12]:
A_xx_val, A_xy_val, A_yx_val, A_yy_val = rotate_polar2cart_tensor2(A_ss_exp, A_sp_exp, A_ps_exp, A_pp_exp)
display(sp.Eq(A_xx, A_xx_val.rewrite(sp.exp).expand().subs(polar2cart_map).simplify()))
display(sp.Eq(A_xy, A_xy_val.rewrite(sp.exp).expand().subs(polar2cart_map).simplify()))
display(sp.Eq(A_yx, A_yx_val.rewrite(sp.exp).expand().subs(polar2cart_map).simplify()))
display(sp.Eq(A_yy, A_yy_val.rewrite(sp.exp).expand().subs(polar2cart_map).simplify()))
$\displaystyle A_{xx} = - \frac{x^{2} y}{x^{2} + y^{2}}$
$\displaystyle A_{xy} = \frac{x^{3}}{x^{2} + y^{2}}$
$\displaystyle A_{yx} = \frac{x^{3}}{x^{2} + y^{2}}$
$\displaystyle A_{yy} = \frac{y \left(2 x^{2} + y^{2}\right)}{x^{2} + y^{2}}$
Different from previous examples, this field is not only bounded, but also at least continuous. However, its derivative is undefined at the origin.
The whole function is not continuously differentiable, so it is still singular.
In [15]:
A_ss_numerical = sp.lambdify((s, p), A_ss_val)
A_sp_numerical = sp.lambdify((s, p), A_sp_val)
A_ps_numerical = sp.lambdify((s, p), A_ps_val)
A_pp_numerical = sp.lambdify((s, p), A_pp_val)
fig, ax = plottings.polar_singularity_rank2tensor(
tfunc_ss=A_ss_numerical, tfunc_sp=A_sp_numerical,
tfunc_ps=A_ps_numerical, tfunc_pp=A_pp_numerical
)
plt.show()
Example: confirmation of the |m|=1 condition¶
Now we change slightly the $A_{ss}$ component from $\cos$ to $\sin$:
In [16]:
A_ss_val = s*sp.sin(p)
A_pp_val = sp.S.Zero
A_sp_val = A_ps_val = s*sp.cos(p)/2
display(sp.Eq(A_ss, A_ss_val))
display(sp.Eq(A_sp, A_sp_val))
display(sp.Eq(A_ps, A_ps_val))
display(sp.Eq(A_pp, A_pp_val))
A_ss_exp = A_ss_val.rewrite(sp.exp).expand()
A_sp_exp = A_sp_val.rewrite(sp.exp).expand()
A_ps_exp = A_ps_val.rewrite(sp.exp).expand()
A_pp_exp = A_pp_val.rewrite(sp.exp).expand()
display(sp.Eq(A_ss, A_ss_exp))
display(sp.Eq(A_sp, A_sp_exp))
display(sp.Eq(A_ps, A_ps_exp))
display(sp.Eq(A_pp, A_pp_exp))
$\displaystyle A_{ss} = s \sin{\left(\phi \right)}$
$\displaystyle A_{s\phi} = \frac{s \cos{\left(\phi \right)}}{2}$
$\displaystyle A_{\phi s} = \frac{s \cos{\left(\phi \right)}}{2}$
$\displaystyle A_{\phi\phi} = 0$
$\displaystyle A_{ss} = - \frac{i s e^{i \phi}}{2} + \frac{i s e^{- i \phi}}{2}$
$\displaystyle A_{s\phi} = \frac{s e^{i \phi}}{4} + \frac{s e^{- i \phi}}{4}$
$\displaystyle A_{\phi s} = \frac{s e^{i \phi}}{4} + \frac{s e^{- i \phi}}{4}$
$\displaystyle A_{\phi\phi} = 0$
Now this field satisfies $A_{s\phi}^{m0} + A_{\phi s}^{m0} = i \mathrm{sgn}(m) (A_{ss}^{m0} - A_{\phi\phi}^{m0})$. As we convert it to Cartesian coordinates, it takes the form
In [17]:
A_xx_val, A_xy_val, A_yx_val, A_yy_val = rotate_polar2cart_tensor2(A_ss_exp, A_sp_exp, A_ps_exp, A_pp_exp)
display(sp.Eq(A_xx, A_xx_val.rewrite(sp.exp).expand().subs(polar2cart_map).simplify()))
display(sp.Eq(A_xy, A_xy_val.rewrite(sp.exp).expand().subs(polar2cart_map).simplify()))
display(sp.Eq(A_yx, A_yx_val.rewrite(sp.exp).expand().subs(polar2cart_map).simplify()))
display(sp.Eq(A_yy, A_yy_val.rewrite(sp.exp).expand().subs(polar2cart_map).simplify()))
$\displaystyle A_{xx} = 0$
$\displaystyle A_{xy} = \frac{x}{2}$
$\displaystyle A_{yx} = \frac{x}{2}$
$\displaystyle A_{yy} = y$
And the whole field becomes perfect regular.
In [25]:
A_ss_numerical = sp.lambdify((s, p), A_ss_val)
A_sp_numerical = sp.lambdify((s, p), A_sp_val)
A_ps_numerical = sp.lambdify((s, p), A_ps_val)
A_pp_numerical = sp.lambdify((s, p), A_pp_val)
fig, ax = plottings.polar_singularity_rank2tensor(
tfunc_ss=A_ss_numerical, tfunc_sp=A_sp_numerical,
tfunc_ps=A_ps_numerical, tfunc_pp=A_pp_numerical
)
plt.show()
Example: a singular field showing necessity of the 2nd-order condition for |m|>1¶
The old conditions involve no additional requirements at the second order for $|m|\geq 2$. To test if the new condition is really necessary, let us look at a field
In [27]:
A_ss_val = (+1 + s**2)*sp.sin(2*p)
A_pp_val = (-1 - s**2)*sp.sin(2*p)
A_sp_val = A_ps_val = (1 + s**2)*sp.cos(2*p)
display(sp.Eq(A_ss, A_ss_val))
display(sp.Eq(A_sp, A_sp_val))
display(sp.Eq(A_ps, A_ps_val))
display(sp.Eq(A_pp, A_pp_val))
A_ss_exp = A_ss_val.rewrite(sp.exp).expand()
A_sp_exp = A_sp_val.rewrite(sp.exp).expand()
A_ps_exp = A_ps_val.rewrite(sp.exp).expand()
A_pp_exp = A_pp_val.rewrite(sp.exp).expand()
display(sp.Eq(A_ss, A_ss_exp))
display(sp.Eq(A_sp, A_sp_exp))
display(sp.Eq(A_ps, A_ps_exp))
display(sp.Eq(A_pp, A_pp_exp))
$\displaystyle A_{ss} = \left(s^{2} + 1\right) \sin{\left(2 \phi \right)}$
$\displaystyle A_{s\phi} = \left(s^{2} + 1\right) \cos{\left(2 \phi \right)}$
$\displaystyle A_{\phi s} = \left(s^{2} + 1\right) \cos{\left(2 \phi \right)}$
$\displaystyle A_{\phi\phi} = \left(- s^{2} - 1\right) \sin{\left(2 \phi \right)}$
$\displaystyle A_{ss} = - \frac{i s^{2} e^{2 i \phi}}{2} + \frac{i s^{2} e^{- 2 i \phi}}{2} - \frac{i e^{2 i \phi}}{2} + \frac{i e^{- 2 i \phi}}{2}$
$\displaystyle A_{s\phi} = \frac{s^{2} e^{2 i \phi}}{2} + \frac{s^{2} e^{- 2 i \phi}}{2} + \frac{e^{2 i \phi}}{2} + \frac{e^{- 2 i \phi}}{2}$
$\displaystyle A_{\phi s} = \frac{s^{2} e^{2 i \phi}}{2} + \frac{s^{2} e^{- 2 i \phi}}{2} + \frac{e^{2 i \phi}}{2} + \frac{e^{- 2 i \phi}}{2}$
$\displaystyle A_{\phi\phi} = \frac{i s^{2} e^{2 i \phi}}{2} - \frac{i s^{2} e^{- 2 i \phi}}{2} + \frac{i e^{2 i \phi}}{2} - \frac{i e^{- 2 i \phi}}{2}$
This field have all the correct prefactors at $|m|=2$, and the lowest order coefficient in $s$ satisfy the conditions required by Holdenried-Chernoff (2021), but the second-order relation is not satisfied, i.e. $A_{s\phi}^{m1} + A_{\phi s}^{m1} \neq i \mathrm{sgn}(m) (A_{ss}^{m1} - A_{\phi\phi}^{m1})$. As we convert it to Cartesian coordinates, it takes the form
In [28]:
A_xx_val, A_xy_val, A_yx_val, A_yy_val = rotate_polar2cart_tensor2(A_ss_exp, A_sp_exp, A_ps_exp, A_pp_exp)
A_xx_xy = A_xx_val.rewrite(sp.exp).expand().subs(polar2cart_map).simplify()
A_xy_xy = A_xy_val.rewrite(sp.exp).expand().subs(polar2cart_map).simplify()
A_yx_xy = A_yx_val.rewrite(sp.exp).expand().subs(polar2cart_map).simplify()
A_yy_xy = A_yy_val.rewrite(sp.exp).expand().subs(polar2cart_map).simplify()
display(sp.Eq(A_xx, A_xx_xy))
display(sp.Eq(A_xy, A_xy_xy))
display(sp.Eq(A_yx, A_yx_xy))
display(sp.Eq(A_yy, A_yy_xy))
$\displaystyle A_{xx} = 0$
$\displaystyle A_{xy} = x^{2} + y^{2} + 1$
$\displaystyle A_{yx} = x^{2} + y^{2} + 1$
$\displaystyle A_{yy} = 0$
In [29]:
dAxx_dx2 = sp.diff(A_xx_xy, (x, 2)).subs(cart2polar_map).simplify()
dAxx_dy2 = sp.diff(A_xx_xy, (y, 2)).subs(cart2polar_map).simplify()
display(dAxx_dx2, dAxx_dy2)
dAxy_dx2 = sp.diff(A_xy_xy, (x, 2)).subs(cart2polar_map).simplify()
dAxy_dy2 = sp.diff(A_xy_xy, (y, 2)).subs(cart2polar_map).simplify()
display(dAxy_dx2, dAxy_dy2)
dAyy_dx2 = sp.diff(A_yy_xy, (x, 2)).subs(cart2polar_map).simplify()
dAyy_dy2 = sp.diff(A_yy_xy, (y, 2)).subs(cart2polar_map).simplify()
display(dAyy_dx2, dAyy_dy2)
$\displaystyle 0$
$\displaystyle 0$
$\displaystyle 2$
$\displaystyle 2$
$\displaystyle 0$
$\displaystyle 0$
The field is again bounded and continuous, but is not continuously differentiable, and hence still singular.
In [28]:
A_ss_numerical = sp.lambdify((s, p), A_ss_val)
A_sp_numerical = sp.lambdify((s, p), A_sp_val)
A_ps_numerical = sp.lambdify((s, p), A_ps_val)
A_pp_numerical = sp.lambdify((s, p), A_pp_val)
fig, ax = plottings.polar_singularity_rank2tensor(
tfunc_ss=A_ss_numerical, tfunc_sp=A_sp_numerical,
tfunc_ps=A_ps_numerical, tfunc_pp=A_pp_numerical
)
plt.show()
In [32]:
A_fields = [
# [sp.lambdify((s, p), A_xx_val), sp.lambdify((s, p), A_xy_val), sp.lambdify((s, p), A_yy_val)],
[sp.lambdify((s, p), dAxx_dx2), sp.lambdify((s, p), dAxy_dx2), sp.lambdify((s, p), dAyy_dx2)]
]
labels = [
# ['$T_{xx}$', '$T_{xy} = T_{yx}$', '$T_{yy}$']
['$d^2T_{xx}/dx^2$', '$d^2T_{xy}/dx^2$', '$d^2T_{yy}/dx^2$']
]
A_fields = np.array(A_fields)
fig, axes = plt.subplots(nrows=A_fields.shape[0], ncols=A_fields.shape[1], figsize=(18, 5.2))
axes = np.atleast_2d(axes)
X_mesh, Y_mesh, S_mesh, P_mesh = plottings.polar_mesh(shape='Cartesian', num=200)
A_field_vals = np.array([[A_tmp(S_mesh, P_mesh)*np.ones_like(S_mesh, dtype=np.float64) for A_tmp in A_field_r] for A_field_r in A_fields])
vmin = min([np.real(A_tmp_val).min() for A_tmp_val in A_field_vals])
vmax = max([np.real(A_tmp_val).max() for A_tmp_val in A_field_vals])
plot_kw = {'shading': 'gouraud', 'vmin': vmin, 'vmax': vmax, 'cmap': 'inferno'}
for i_r, ax_r in enumerate(axes):
for i_c, ax in enumerate(ax_r):
im = ax.pcolormesh(X_mesh, Y_mesh, np.real(A_field_vals[i_r, i_c]), **plot_kw)
ax.set_aspect('equal', 'box')
ax.set_title(labels[i_r][i_c], fontsize=20)
if i_c == 0:
ax.set_ylabel('$y$', fontsize=18)
if i_r == A_fields.shape[0] - 1:
ax.set_xlabel('$x$', fontsize=18)
ax.tick_params(axis='both', which='major', labelsize=14)
# axes[0].set_xlabel('$x$', fontsize=16)
# axes[0].set_ylabel('$y$', fontsize=16)
# axes[0].set_title(r'$dV_x/dx$', fontsize=18)
# axes[1].set_xlabel('$x$', fontsize=16)
# axes[1].set_title(r'$dV_x/dy$', fontsize=18)
fig.subplots_adjust(right=0.87)
cbar_ax = fig.add_axes([0.9, 0.15, 0.02, 0.7])
cbar_ax.tick_params(axis='both', which='major', labelsize=14)
fig.colorbar(im, cax=cbar_ax)
save_name = '../out/imgs/Regular_tensor2_m2_derivative'
# save_name = None
if save_name is not None:
plt.savefig(save_name + '.png', dpi=150, bbox_inches='tight')
plt.savefig(save_name + '.pdf', bbox_inches='tight')
plt.show()
Example: confirmation of the 2nd-order condition for |m|>1¶
We only need to flip the sign for the second order term in $A_{\phi\phi}$ to construct a field which completely abides by the newly derived regularity conditions, even at second order:
In [5]:
A_ss_val = (+1 + s**2)*sp.sin(2*p)
A_pp_val = (-1 - s**2)*sp.sin(2*p)
A_sp_val = A_ps_val = (1 + s**2)*sp.cos(2*p)
display(sp.Eq(A_ss, A_ss_val))
display(sp.Eq(A_sp, A_sp_val))
display(sp.Eq(A_ps, A_ps_val))
display(sp.Eq(A_pp, A_pp_val))
A_ss_exp = A_ss_val.rewrite(sp.exp).expand()
A_sp_exp = A_sp_val.rewrite(sp.exp).expand()
A_ps_exp = A_ps_val.rewrite(sp.exp).expand()
A_pp_exp = A_pp_val.rewrite(sp.exp).expand()
display(sp.Eq(A_ss, A_ss_exp))
display(sp.Eq(A_sp, A_sp_exp))
display(sp.Eq(A_ps, A_ps_exp))
display(sp.Eq(A_pp, A_pp_exp))
$\displaystyle A_{ss} = \left(s^{2} + 1\right) \sin{\left(2 \phi \right)}$
$\displaystyle A_{s\phi} = \left(s^{2} + 1\right) \cos{\left(2 \phi \right)}$
$\displaystyle A_{\phi s} = \left(s^{2} + 1\right) \cos{\left(2 \phi \right)}$
$\displaystyle A_{\phi\phi} = \left(- s^{2} - 1\right) \sin{\left(2 \phi \right)}$
$\displaystyle A_{ss} = - \frac{i s^{2} e^{2 i \phi}}{2} + \frac{i s^{2} e^{- 2 i \phi}}{2} - \frac{i e^{2 i \phi}}{2} + \frac{i e^{- 2 i \phi}}{2}$
$\displaystyle A_{s\phi} = \frac{s^{2} e^{2 i \phi}}{2} + \frac{s^{2} e^{- 2 i \phi}}{2} + \frac{e^{2 i \phi}}{2} + \frac{e^{- 2 i \phi}}{2}$
$\displaystyle A_{\phi s} = \frac{s^{2} e^{2 i \phi}}{2} + \frac{s^{2} e^{- 2 i \phi}}{2} + \frac{e^{2 i \phi}}{2} + \frac{e^{- 2 i \phi}}{2}$
$\displaystyle A_{\phi\phi} = \frac{i s^{2} e^{2 i \phi}}{2} - \frac{i s^{2} e^{- 2 i \phi}}{2} + \frac{i e^{2 i \phi}}{2} - \frac{i e^{- 2 i \phi}}{2}$
As we convert it to Cartesian coordinates, it takes the form
In [30]:
A_xx_val, A_xy_val, A_yx_val, A_yy_val = rotate_polar2cart_tensor2(A_ss_exp, A_sp_exp, A_ps_exp, A_pp_exp)
display(sp.Eq(A_xx, A_xx_val.rewrite(sp.exp).expand().subs(polar2cart_map).simplify()))
display(sp.Eq(A_xy, A_xy_val.rewrite(sp.exp).expand().subs(polar2cart_map).simplify()))
display(sp.Eq(A_yx, A_yx_val.rewrite(sp.exp).expand().subs(polar2cart_map).simplify()))
display(sp.Eq(A_yy, A_yy_val.rewrite(sp.exp).expand().subs(polar2cart_map).simplify()))
$\displaystyle A_{xx} = 0$
$\displaystyle A_{xy} = x^{2} + y^{2} + 1$
$\displaystyle A_{yx} = x^{2} + y^{2} + 1$
$\displaystyle A_{yy} = 0$
The field is now polynomial and undoubtedly regular throughout the domain
In [31]:
A_ss_numerical = sp.lambdify((s, p), A_ss_val)
A_sp_numerical = sp.lambdify((s, p), A_sp_val)
A_ps_numerical = sp.lambdify((s, p), A_ps_val)
A_pp_numerical = sp.lambdify((s, p), A_pp_val)
fig, ax = plottings.polar_singularity_rank2tensor(
tfunc_ss=A_ss_numerical, tfunc_sp=A_sp_numerical,
tfunc_ps=A_ps_numerical, tfunc_pp=A_pp_numerical
)
plt.show()
Spectral basis¶
In [60]:
r_array = np.linspace(0, 1, num=100)
phi_array = np.linspace(0., 2*np.pi, num=100)
X_mesh = np.outer(r_array, np.cos(phi_array))
Y_mesh = np.outer(r_array, np.sin(phi_array))
def plot_func(func, ax, vclip=None, title=None, **kwargs):
F_val = func(r_array, phi_array)
if vclip is None:
vclip = np.max(np.abs(F_val))
ax.pcolormesh(X_mesh, Y_mesh, F_val, vmin=-vclip, vmax=vclip, **kwargs)
if title is not None:
ax.set_title(title, fontsize=16)
ax.set_aspect('equal', 'box')
ax.axis('off')
In [123]:
from scipy.special import eval_jacobi
class ZernikePolynomial:
def __init__(self, n, m):
self.n = n
self.m = m
assert (self.n >= abs(self.m)) and (self.n - abs(self.m))%2 == 0
self.n_j = (n - abs(m)) // 2
def __call__(self, r_array, phi_array):
f_phi = np.cos(self.m*phi_array) if self.m >= 0 else np.sin(abs(self.m)*phi_array)
f_r = (-1)**self.n_j*(r_array**abs(self.m))*eval_jacobi(self.n_j, abs(self.m), 0, 1 - 2*r_array**2)
return np.outer(f_r, f_phi)
class OnesidedJacobiPoly:
def __init__(self, n, m, alpha, beta):
self.n, self.m = n, m
self.a, self.b = alpha, beta
assert (self.n >= abs(self.m)) and (self.n - abs(self.m))%2 == 0
self.n_j = (n - abs(m)) // 2
def __call__(self, r_array, phi_array):
f_phi = np.cos(self.m*phi_array) if self.m >= 0 else np.sin(abs(self.m)*phi_array)
f_r = (r_array**abs(self.m))*eval_jacobi(self.n_j, self.a, self.b, 1 - 2*r_array**2)
return np.outer(f_r, f_phi)
class SkewedOnesidedJacobiPoly:
def __init__(self, n, m, skew, alpha, beta):
self.n, self.m = n, m
self.rp = abs(m + skew)
self.a, self.b = alpha, beta
# assert (self.n >= self.rp) and (self.n - self.rp)%2 == 0
self.n_j = (n - abs(m)) // 2
def __call__(self, r_array, phi_array):
# f_phi = np.cos(self.m*phi_array) if self.m >= 0 else np.sin(abs(self.m)*phi_array)
f_phi = np.cos(self.m*phi_array)
f_r = (r_array**self.rp)*eval_jacobi(self.n_j, self.a, self.b, 1 - 2*r_array**2)
return np.outer(f_r, f_phi)
In [125]:
Nrows = 6
Ngrids = 10
Ninch = 2.5
fwidth = Ngrids*Nrows
mid_pt = fwidth//2
save_name = '../out/imgs/Skewed_Jacobi_n_N5M5'
# save_name = None
# vclip = 0
# for i_r in range(Nrows):
# for i_c in range(-i_r, i_r+1, 2):
# # F_val = ZernikePolynomial(i_r, i_c)(r_array, phi_array)
# F_val = OnesidedJacobiPoly(i_r, i_c, 3/2, 3/2)(r_array, phi_array)
# if np.max(np.abs(F_val)) > vclip:
# vclip = np.max(np.abs(F_val))
vclip = None
fig = plt.figure(figsize=(Ninch*Nrows, Ninch*Nrows))
gs = fig.add_gridspec(nrows=fwidth, ncols=fwidth, wspace=5, hspace=5)
plot_kw = {'shading': 'gouraud', 'cmap': 'RdBu_r'}
for i_r in range(Nrows):
row_range = slice(Ngrids*i_r, Ngrids*(i_r + 1))
for i_c in range(-i_r, i_r+1, 2):
col_range = slice(mid_pt + (i_c-1)*(Ngrids//2), mid_pt + (i_c+1)*(Ngrids//2))
ax = fig.add_subplot(gs[row_range, col_range])
# plot_func(ZernikePolynomial(i_r, i_c), ax=ax, title=r'$Z_{%d}^{%d}$' % (i_r, i_c), vclip=vclip, **plot_kw)
# plot_func(OnesidedJacobiPoly(i_r, i_c, abs(i_c)-1/2, -1/2), ax=ax, title=r'$\Psi_{%d}^{%d}$' % (i_r, i_c), vclip=vclip, **plot_kw)
# plot_func(OnesidedJacobiPoly(i_r, i_c, 3/2, 3/2), ax=ax, title=r'$r^{%d} P_{%d}^{(\frac{1}{2}, \frac{1}{2})}(1 - 2r^2) \cos {%d}\phi$' % (i_r, abs(i_c), abs(i_c)), vclip=vclip, **plot_kw)
plot_func(SkewedOnesidedJacobiPoly(i_r, i_c, -1, abs(i_c-1), 0), ax=ax, title=r'$\Psi_{-,%d}^{%d}$' % (i_r, i_c), vclip=vclip, **plot_kw)
if save_name is not None:
plt.savefig(save_name + '.pdf', bbox_inches='tight')
plt.savefig(save_name + '.png', dpi=150, bbox_inches='tight')
plt.show()
In [ ]: