一変数正則関数の零点と臨界点を可視化

はじめに

この記事では、前回投稿した「Pythonで連続関数の零点を網羅的に求める」で紹介した連続関数の零点を求める関数SearchZeroPtByDepthの一つの応用例を記載する。

目的

複素平面内の円板で定義された非定数の一変数正則関数

f(z)\space\space (z \in D_R)\\

\left(
\begin{eqnarray}

D_R &:=& \{|z| < R\}\space\space (R>0)\\
z &=& x + iy\\
&=& re^{i\theta}\space\space (0 < r < R)\\
\end{eqnarray}\right)

に対して、 その絶対値の２乗 $$A(z) := |f(z)|^2 = f(z)\overline{f(z)}$$ の臨界点を網羅的に求めるという問題設定をする。ここでは先ず数値計算によって、臨界点がどのように分布しているのかを観察してみる。

動機

はじめに、今後使用すると思われる、領域に関する記法についてまとめておく。

\begin{eqnarray}

r&>&0のとき\\

D_r &:=& \{|z| < r\}\space\space\space開円板\\
\overline{D_r} &:=& \{|z| \leqq r\}\space\space\space閉円板\\
\partial{D_r} &:=& \{|z| = r\}\space\space\space円板の境界（円）\\
\end{eqnarray}

ここで $r<R$ ならば、$z \in \overline{D_r}$で$f(z)$ は定数でない正則関数なので、最大絶対値の原理より、$|f(z)|$は$D_r$の内部では最大値をとらず、境界 $\partial{D_r}$で最大値をとる。これは、$A(z)=|f(z)|^2$なので、$A(z)$についても同様である。したがって、

M_r=max\{A(z)|z \in \partial{D_r}\}

とおくと、

r_1 < r_2 \Longrightarrow
M_{r_1} < M_{r_2} 

が成り立つ。
このことから始めは、$r$を固定したときの最大値になる点$(z,A(z))$が、$r$の変化に伴いどのような軌跡を描きながら増大していくのか興味を持った。そして、最大値だけでなく、より一般に臨界値になる点$(z,A(z))$の軌跡、最終的には、$\theta$を固定したときの臨界値になる点$(z,A(z))$が、$\theta$の変化に伴い描く軌跡も合わせて考えるようになった。

定義

ここで、用語について以下に定義しておく。
ただし、偏向臨界点（偏向臨界集合）は造語であり、ここだけで通用するものであることを断っておく。

定義1-1．臨界点、臨界値

$z$が正則関数$f$の臨界点とは、$z$が次を満たすこと

\frac{df}{dz}(z)=0\\

数式1によれば、この定義は写像

\phi:(x,y)\longmapsto(u(x,y),v(x,y))

の臨界点の定義と一致する。

定義1-2．A(z)の極座標臨界点、極座標臨界値

$z=re^{i\theta}\space\space(r\neq 0)$が$A(z)$の極座標臨界点とは、$z$が次を満たすこと

\begin{equation}
    \left\{
    \begin{aligned}
\frac{\partial A}{\partial r}(z) = 0\\
\frac{\partial A}{\partial \theta}(z) = 0
    \end{aligned}
    \right.
    \quad
\end{equation}

$z$が極座標臨界点のとき、$A(z)$を極座標臨界値と呼ぶ。
このとき、$(z,A(z))$も極座標臨界点と記述することがある。

注意

$r=0$においては数式2から、座標変換 $(r,\theta)\longmapsto(x,y)$のヤコビ行列式が$0$になるため、$r=0$において微分同相にならず

\begin{equation}
    \begin{aligned}
\frac{\partial A}{\partial r}(z) = 
\frac{\partial A}{\partial \theta}(z) = 0
    \end{aligned}
    \quad
\end{equation}

であっても、

\begin{equation}
    \begin{aligned}
\frac{\partial A}{\partial x}(z) = 
\frac{\partial A}{\partial y}(z) = 0
    \end{aligned}
    \quad
\end{equation}

とは限らない。そのような訳で上の定義で$r\neq 0$の条件を付けてある。
したがって、$r=0$つまり$z=0$においては、別途

\begin{equation}
\frac{\partial A}{\partial z}(0) = 0\\
\end{equation}

が成り立つかどうか確認する必要がある。

臨界点と極座標臨界点の関係

数式3によれば、$z\neq0$のとき、
$f(z)$の臨界点と零点を合わせた集合は、$A(z)$の極座標臨界点全体と一致する。
また、数式3と同じだが、微分演算子を使った表現を数式4に記載しておく。後で何かの役に立つかもしれない。

定義2-1．A(z)の偏角方向の偏向臨界点（偏向臨界集合）

$z=re^{i\theta}\space\space(r\neq 0)$が$A(z)$の偏角方向の偏向臨界点とは、$z$が次を満たすこと

\frac{\partial A}{\partial \theta}(z) = 0

このとき、$(z,A(z))$も偏角方向の偏向臨界点と記述することがある。
また、偏角方向の偏向臨界点全体を偏角方向の偏向臨界点集合と記述する。

定義2-2．A(z)の半径方向の偏向臨界点（偏向臨界集合）

$z=re^{i\theta}\space\space(r\neq 0)$が$A(z)$の半径方向の偏向臨界点とは、$z$が次を満たすこと

\frac{\partial A}{\partial r}(z) = 0

このとき、$(z,A(z))$も半径方向の偏向臨界点と記述することがある。
また、半径方向の偏向臨界点全体を半径方向の偏向臨界点集合と記述する。

臨界点の例

2次多項式

$f(z) = z^2 + z - 1$　の場合

左図：$r=0.5$における$A(0.5,\theta)$のグラフ（水色の線）と偏角方向の偏向臨界点（青点）
右図：$\theta=0$における$A(r,0)$のグラフ（紫色の線）と半径方向の偏向臨界点（赤点）

左図：$r$を変化させたときの、偏角方向の偏向臨界点（青点）の軌跡
右図：$\theta$を変化させたときの、半径方向の偏向臨界点（赤点）の軌跡
左図の偏角方向の偏向臨界点の内に、$(re^{i\theta},M_r)$の軌跡が含まれている。

偏角方向の偏向臨界点集合（青点）、半径方向の偏向臨界点集合（赤点）、$(z,A(z))$のグラフを重ねた3D表示
（ただし、零点や臨界点を強調するため、$A(z)$が$2.0$以上のときは$2.0$としてプロット）
偏角方向の偏向臨界点集合（青点）、半径方向の偏向臨界点集合（赤点）の交点が臨界点になる。
例えて言うならば、偏向臨界集合が骨格で、臨界点が関節のように見える。

Pythonプログラム

上記のグラフを描くPythonプログラムのコードは以下のとおり。

プログラム全文(ここをクリックすると表示します)

# -*- coding: utf-8 -*-

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import seaborn as snsni
from DepthMethod import SearchZeroPtByDepth,numerical_diff
import sys

class PartialValue:
    def __init__(self,f):
        self.func = f
        self.fixx = None
        self.fixy = None
        
    def fix_x(self,x):
        self.fixx=x
        
    def fix_y(self,y):
        self.fixy=y
        
    def func_x(self):
        return lambda x:self.func(x , self.fixy)
    
    def func_y(self):
        return lambda y:self.func(self.fixx , y)

#polar to real Cartesian
x=lambda r,t:r*np.cos(t)
y=lambda r,t:r*np.sin(t)

#real to complex
z=lambda x,y:x+y*1j

#definiton of holomorphic function
func,funcName = lambda z:z**2+z-1,'z**2+z-1'

diffdelt = 1e-6
zLim = 10.0 #max value of graph
abs2min = lambda a,zLim:np.minimum(a,zLim)
absfunc = lambda w:(abs(func(w))).real
absfunc2 = lambda w:(func(w)*func(w).conjugate()).real
abs2 = lambda r,t:absfunc2(x(r,t)+y(r,t)*1j).real
absfunc2_min = lambda w,zLim:abs2min(absfunc2(w),zLim)

print('-------------------------------------------')
print('partial critical points by polar coodinate')
print('-------------------------------------------')

"""
偏角方向の偏向臨界点
"""
print('Shift the radius,plot theta critical points')
#Shift the radius
rary=np.linspace(-1.0,1.0,200)
#Shift the theata
tary=np.linspace(0.0,2*np.pi,200)

Len = rary.max()
    
re_t=[]
im_t=[]
fv_t=[]
z_t=[]
re_t_c=[]
im_t_c=[]
fv_t_c=[]

#ra=[0.5]
ra=np.arange(0.05,1.0,0.05)

start = 0.0
ended = 2.5*np.pi

re_t.append(0.0)
im_t.append(0.0) 
fv_t.append(abs2(0,0))
z_t.append(0.0)

for r in ra:
    for t in tary:
        w = r*np.cos(t) + 1j*r*np.sin(t)
        re_t.append(w.real)
        im_t.append(w.imag) 
        #fv_t.append(abs2(r,t))
        fv_t.append(absfunc2_min(w,zLim))
        z_t.append(0.0)

    a = lambda r,t:absfunc2(x(r,t)+y(r,t)*1j)
    sv = PartialValue(a)
    sv.fix_x(r)
    fst = sv.func_y()
    dfst = lambda t:numerical_diff(fst,t,diffdelt)
    print('Radius = ',r)
    daZeros,rdGoodNum=SearchZeroPtByDepth(dfst,'ShiftRadius',start,ended,InitRandNum=50)
    daZeroSel = daZeros[0]
    if (daZeroSel.size > 0):
        #daZero=daZeroSel['EstimateZeroPt'][abs(daZeroSel['FunctionValue']) < np.sqrt(diffdelt)]
        daZero=daZeroSel['EstimateZeroPt']
        w=z(x(r,daZero),y(r,daZero))
        for i in w.index:
            if ((abs(w[i].real)<=Len) and (abs(w[i].imag) <= Len )):
                if (absfunc2(w[i]) < zLim):
                    re_t_c.append(w[i].real)
                    im_t_c.append(w[i].imag)
                    fv_t_c.append(absfunc2(w[i]))
##    add origine into critical points
#    re_t_c.append(0.0)
#    im_t_c.append(0.0)
#    fv_t_c.append(absfunc2(0.0))

"""
半径方向の偏向臨界点
"""
print('Shift the theta,plot radius critical points')

re_r=[]
im_r=[]
fv_r=[]
z_r=[]
re_r_c=[]
im_r_c=[]
fv_r_c=[]

#ta=[0]
ta=np.pi*1/30*np.arange(0,30,1)

start = -Len
ended = Len
for t in ta:
    for r in rary:
        w = r*np.cos(t) + 1j*r*np.sin(t)
        re_r.append(w.real)
        im_r.append(w.imag) 
        #fv_r.append(abs2(r,t))
        fv_r.append(absfunc2_min(w,zLim))
        z_r.append(0.0)       

    a = lambda r,t:absfunc2(x(r,t)+y(r,t)*1j)
    sv = PartialValue(a)
    sv.fix_y(t)
    fst = sv.func_x()
    dfst = lambda r:numerical_diff(fst,r,diffdelt)
    print('Theta = ',t)
    daZeros,rdGoodNum=SearchZeroPtByDepth(dfst,'ShiftAngle',start,ended,InitRandNum=50)
    daZeroSel = daZeros[0]
    if (daZeroSel.size > 0):
        #daZero=daZeroSel['EstimateZeroPt'][abs(daZeroSel['FunctionValue']) < np.sqrt(diffdelt)]
        daZero=daZeroSel['EstimateZeroPt']
        w=z(x(daZero,t),y(daZero,t))
        for i in w.index:
            if ((abs(w[i].real)<=Len) and (abs(w[i].imag) <= Len )):
                if (absfunc2(w[i]) < zLim):
                    re_r_c.append(w[i].real)
                    im_r_c.append(w[i].imag)
                    fv_r_c.append(absfunc2(w[i]))
    
sns.set_style("whitegrid", {'grid.linestyle': '--'})
fig = plt.figure(figsize=plt.figaspect(1.0))
fig.canvas.set_window_title('f(z) = '+funcName)

ax1 = fig.add_subplot(1, 2, 1, projection='3d')
ax1.set_xlabel("real part")
ax1.set_ylabel("imaginary part")
ax1.set_zlabel("absolute value")
#ax1.set_zscale('log',basex=10)
ax1.plot(re_t,im_t,fv_t,marker=".",color="cyan",linestyle='None',label='on cosntant radius',markersize=2)
ax1.plot(re_t,im_t,z_t,marker=".",color="darkgray",linestyle='None',markersize=2)
ax1.plot(re_t_c,im_t_c,fv_t_c,marker=".",color="blue",linestyle='None',label='theta critical points')

plt.xlim(-rary.max(),rary.max())
plt.ylim(-rary.max(),rary.max())
plt.legend()

ax2 = fig.add_subplot(1, 2, 2, projection='3d')
ax2.set_xlabel("real part")
ax2.set_ylabel("imaginary part")
ax2.set_zlabel("absolute value")
#ax2.set_zscale('log',basex=10)
ax2.plot(re_r,im_r,fv_r,marker=".",color="violet",linestyle='None',label='on cosntant theta',markersize=2)
ax2.plot(re_r,im_r,z_r,marker=".",color="darkgray",linestyle='None',markersize=2)
ax2.plot(re_r_c,im_r_c,fv_r_c,marker=".",color="red",linestyle='None',label='radius critical points')
plt.xlim(-rary.max(),rary.max())
plt.ylim(-rary.max(),rary.max())
plt.legend()

plt.show()

"""
偏向臨界集合とAの3D表示
"""
zLim2=2.0
sns.set_style("whitegrid", {'grid.linestyle': '--'})
fig = plt.figure(figsize=plt.figaspect(1.0))
fig.canvas.set_window_title('f(z) = '+funcName)

#for reference
dLen = Len/100
X, Y = np.mgrid[-Len:Len:dLen, -Len:Len:dLen]
#Z = absfunc(z(X, Y))
Z = absfunc2_min(z(X, Y),zLim2)
ax1 = fig.add_subplot(1, 1, 1, projection='3d')
ax1.set_xlabel("real part")
ax1.set_ylabel("imaginary part")
ax1.set_zlabel("absolute value")
plt.xlim(-Len,Len)
plt.ylim(-Len,Len)
ax1.set_zlim(0.0,zLim2)

re_t_c_rest=np.array(re_t_c)
im_t_c_rest=np.array(im_t_c)
fv_t_c_rest=np.array(fv_t_c)
idx=zLim2 > fv_t_c_rest
re_t_c_rest=re_t_c_rest[idx]
im_t_c_rest=im_t_c_rest[idx]
fv_t_c_rest=fv_t_c_rest[idx]

re_r_c_rest=np.array(re_r_c)
im_r_c_rest=np.array(im_r_c)
fv_r_c_rest=np.array(fv_r_c)
idx=zLim2 > fv_r_c_rest
re_r_c_rest=re_r_c_rest[idx]
im_r_c_rest=im_r_c_rest[idx]
fv_r_c_rest=fv_r_c_rest[idx]

ax1.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap='hsv_r', linewidth=0.1)
ax1.plot(re_t_c_rest,im_t_c_rest,fv_t_c_rest,marker=".",color="blue",linestyle='None',label='theta critical points')
ax1.plot(re_r_c_rest,im_r_c_rest,fv_r_c_rest,marker=".",color="red",linestyle='None',label='radius critical points')
plt.legend()

plt.show()

他の正則関数を対象にするには、

func,funcName = lambda z:z**2+z-1,'z**2+z-1'

において、λ式を変更すればよい。
また、zLim、zLim2を変更することで、臨界点付近の変化を強調することができる。

数式

数式1

\begin{eqnarray}
f(z) &:=& u(x,y) + iv(x,y)\\
J_f(z) &:=& 
\begin{vmatrix}
u_x & u_y \\
v_x & v_y \\
\end{vmatrix}\\
&=&u_xv_y - u_yv_x\\
&=&u_x^2 + v_x^2\\
&=&(u_x+iv_x)\overline{(u_x+iv_x)}\\
&=&\frac{df(z)}{dz}\overline{\frac{df(z)}{dz}}\\
\end{eqnarray}\\
\space\space\\
\therefore \space\space J_f(z)=0 \Leftrightarrow \frac{df(z)}{dz}=0

数式2

\begin{eqnarray}

\begin{pmatrix}
dx\\
dy
\end{pmatrix}
&=&
\begin{pmatrix}
cos\theta & -rsin\theta \\
sin\theta & rcos\theta \\
\end{pmatrix}
\begin{pmatrix}
dr\\
d\theta
\end{pmatrix}\\

J &:=& 
\begin{vmatrix}
cos\theta & -rsin\theta \\
sin\theta & rcos\theta \\
\end{vmatrix}
=r
\end{eqnarray}

数式3

$f(z)$が正則関数のとき、下記の式が成り立つ。部分的には、数式3-1～数式3-3を使う。
$\Re$、$\Im$は、それぞれ複素数の実部、虚部を表す。

\begin{eqnarray}

\frac{\partial A}{\partial r}&=&\frac{\partial A}{\partial x}\frac{\partial x}{\partial r}+\frac{\partial A}{\partial y}\frac{\partial y}{\partial r}\\
&=&(\frac{\partial A}{\partial z}\frac{\partial z}{\partial x}+\frac{\partial A}{\partial \overline{z}}\frac{\partial \overline{z}}{\partial x})\frac{\partial x}{\partial r}+(\frac{\partial A}{\partial z}\frac{\partial z}{\partial y}+\frac{\partial A}{\partial \overline{z}}\frac{\partial \overline{z}}{\partial y})\frac{\partial y}{\partial r}\\&=&2\Re\left(e^{i\theta}.\frac{\partial f(z)}{\partial z}.\overline{f(z)}\right)\\
\\
\frac{\partial A}{\partial \theta}&=&\frac{\partial A}{\partial x}\frac{\partial x}{\partial \theta}+\frac{\partial A}{\partial y}\frac{\partial y}{\partial \theta}\\
&=&(\frac{\partial A}{\partial z}\frac{\partial z}{\partial x}+\frac{\partial A}{\partial \overline{z}}\frac{\partial \overline{z}}{\partial x})\frac{\partial x}{\partial \theta}+(\frac{\partial A}{\partial z}\frac{\partial z}{\partial y}+\frac{\partial A}{\partial \overline{z}}\frac{\partial \overline{z}}{\partial y})\frac{\partial y}{\partial \theta}\\
&=&-2r\Im\left(e^{i\theta}.\frac{\partial f(z)}{\partial z}.\overline{f(z)}\right)\\

\end{eqnarray}\\

数式3-1

\begin{eqnarray}

\frac{\partial z}{\partial x}=1,\frac{\partial \overline{z}}{\partial x}=1,
\frac{\partial z}{\partial y}=i,\frac{\partial \overline{z}}{\partial y}=-i\\
\end{eqnarray}\\

数式3-2

\begin{eqnarray}
\frac{\partial x}{\partial r} = cos\theta, \frac{\partial y}{\partial r} = sin\theta,\frac{\partial x}{\partial \theta} =-rsin\theta, \frac{\partial y}{\partial \theta}=rcos\theta\\
\end{eqnarray}

数式3-3

\begin{eqnarray}

\frac{\partial A}{\partial z}&=&\frac{\partial}{\partial z}(f(z)\overline{f(z)})\\
&=&\frac{\partial f(z)}{\partial z}.\overline{f(z)}+f(z).\frac{\partial\overline{f(z)}}{\partial z}\\
&=&\frac{\partial f(z)}{\partial z}.\overline{f(z)}\\
\\
\frac{\partial A}{\partial \overline{z}}&=&\frac{\partial}{\partial \overline{z}}(f(z)\overline{f(z)})\\
&=&\frac{\partial f(z)}{\partial \overline{z}}.\overline{f(z)}+f(z).\frac{\partial\overline{f(z)}}{\partial \overline{z}}\\
&=&f(z).\frac{\partial\overline{f(z)}}{\partial \overline{z}}\\
&=&f(z).\overline{\frac{\partial f(z)}{\partial z}}\\
&=&\overline{\frac{\partial f(z)}{\partial z}.\overline{f(z)}}\\
\end{eqnarray}
\\
上では、fが正則関数なので、\frac{\partial f(z)}{\partial \overline{z}}=\frac{\partial\overline{f(z)}}{\partial z}=0を使っている。

数式4

\frac{1}{2}\left(\frac{\partial}{\partial r} - \frac{i}{r}\frac{\partial}{\partial \theta}\right)A = e^{i\theta}.\frac{\partial f(z)}{\partial z}.\overline{f(z)}

参考

１の応用としてこの記事を書いた。
3Dグラフの描画については、2を参考にさせていただいた。

１．Pythonで連続関数の零点を網羅的に求める
２．[pythonで3次元の散布図]
(https://qiita.com/pegasusBS15/items/3445c5ad1a70f61c85d4)