基本解近似法によるクラドニ図形の計算 -

基本解近似法は、日本語文献だと、代用電荷法の方が通りがいいっぽい。理由としては、日本語の唯一の成書
代用電荷法とその応用 : 境界値問題の半解析的近似解法
http://webcatplus.nii.ac.jp/webcatplus/details/book/ncid/BN0003658X.html
が、代用電荷法と呼んでるからだと思う。けど、別に、電気工学に限った方法というわけではないので、もうちょっと分野非依存的な名前の方が好ましい。英語文献では、method of fundamental solution(略して、MFS)の方が良く使われてるっぽくて、直訳すると、"基本解の方法"とかだけど、日本語としては、なんか違和感があるので、基本解近似法と呼んでおく。

名前の違いは、独立した考案者が存在したせいらしく、数学の文献では Kupradze and Aleksidze (1964)，電気工学の文献ではSteinbiegler (1969，学位論文)により提唱されたらしい。

元々は、ユークリッド空間内の(有界とは限らない)真部分領域上のポテンシャル問題の高速解法として提案されたっぽい。問題の解を、外部領域にソースがある基本解の線形和で表し(従って、内部領域と境界上では、微分作用素の作用によって消える)、線形結合の係数を境界条件をなるべくよく満たすように決定するという感じの方法。外部領域のどこにソースを置くかは、ユーザが適当に決めるが、普通、境界に近いとこに置く。境界条件は、Dirichlet境界条件、Neumann境界条件、混合型境界条件等どれでもいいけど、これも境界上に複数の代表点(これもユーザが適当に決める)を取って、例えば、Dirichlet境界条件なら、代表点で、境界値を与える関数と一致するという条件を解くというのが、オーソドックスな方法。

ソースを外部領域に配置する都合上、問題領域が、全ユークリッド空間だと使えない。また、基本解が解析的に分かる微分作用素というのは、さほど多いとは言えないというあたりが弱点。

一方、差分法や有限要素法では面倒な非有界領域でも使えるし、有限要素法や境界要素法のように領域をイイ感じに分割しなくてよくて、実装が割と簡単なこと(特殊関数の値からなる行列を作って、適当に解くだけ)が長所。

多分、計算量とかの点では、基本解近似法と境界要素法は似てる部分が多いと思うけど、基本解近似法は、境界要素の分割も必要ない。差分法や有限要素法、そして、SPHのような粒子法は、3次元より大きい次元での計算には、計算量以前にメモリ消費の観点から考えて、適しているとは言えない。粒子法に関しては、かなり粗い離散化も原理的には可能であるが、少なくとも、私の経験上では、有限要素法や差分法と同程度の精度を得るには、同程度のメモリを使う必要がある。

実際上は、3次元より大きい次元での計算が必要になることは通常ない(※)ので、これが凄くありがたいわけではないけれど。

※)例外の一つは、ボルツマン方程式などの輸送方程式で、6次元での計算が必要になる。別の例外は、量子化学で、電子数Nの系を考える場合、(原子核の自由度は分離できるとして)3N次元空間上の波動関数を考える必要がある。量子化学では通常の離散化手法は全く使えず、良い基底で解を展開し、多くの場合、最低固有値のみを計算する(頑張れば、励起状態の計算もできるけど、いずれにせよ、1つか2つの固有値を計算するのが精々)。多分、原子核物理も、同様の困難があるはずだけど、どうしているのかは知らない。これらの分野の計算で、基本解近似法が役に立つことは多分ないけど、いい基底を見つけて展開するというのが、離散化の本命なのかもしれない

基本解近似法をDirichlet固有値問題なんかに適用する場合には、固有関数をHelmholtz方程式の基本解の一次結合として書く。しかし、Helmholtz方程式の基本解は固有値に依存してるので、散乱問題のように連続スペクトルが出る場合はともかく、離散スペクトルを決定するには、境界要素法でDirichlet固有値問題を解く時と同様、Determinant Searchなどをやることになるらしい。

簡単な例として、単位円板 $D=\{z \in \mathbf{C} | |z|=1 \}$ 上でのLaplacianに対するDirichlet固有値問題
$\Delta \phi = -\lambda \phi$
$\phi(z)=0 (z \in \partial D)$
を考える。よく知られていることとして、この問題の固有値は、第一種Bessel関数の零点の二乗として書ける。第一種Bessel関数の零点を、いくつか掲載すると、以下のようになってる

	n=0	n=1	n=2	n=3	n=4	n=5	n=6	n=7
p=1	2.4048	3.8317	5.1356	6.3802	7.5883	8.7715	9.9361	11.0864
p=2	5.5201	7.0156	8.4172	9.7610	11.0647	12.3386	13.5893	14.8213
p=3	8.6537	10.1735	11.6198	13.0152	14.3725	15.7002	17.0038	18.2876
p=4	11.7915	13.3237	14.7960	16.2235	17.6160	18.9801	20.3208	21.6415
p=5	14.9309	16.4706	17.9598	19.4094	20.8269	22.2178	23.5861	24.9439

nはBessel関数の次数で、pは、零点を小さい順に並べた時のランク。100以下の固有値(実際の固有値は、上の表の値を二乗したものであることに注意)は12個出ることが分かる。

基本解近似法に基づく計算を、以下に実装してみた。固有値かどうかの指標として、近似固有関数の(近似)ノルムを見ている。固有値じゃないスペクトルに対応する解は0しかないので、固有値から離れたところで近似解を探しても、ノルムの小さい解しか出ないと期待される。横軸は波数、縦軸は近似ノルムの値で、波数ごとに近似ノルム出してるだけで、具体的な固有値の数値を出す所までやってないけど、適当にピークピッキングするとかすればいいはず。

f:id:m-a-o:20210809122302p:plain

実装の詳細は、色々な方法がありえると思うけど、

The method of fundamental solutions for Helmholtz eigenvalue problems in simply and multiply connected domains
https://doi.org/10.1016/j.enganabound.2005.08.011

を参考にした。k-procedureと書かれてる方法は必須ではないけど、とりあえず使ってみた。論文では、固有値を検索するのに、Algorithm Aというのを書いてあるけど、これは面倒なので実装してない。単に、波数kを一定ステップで刻んで、近似固有関数のノルムを得ている。実装は、以下のような感じ

from scipy.spatial import distance_matrix
from scipy.special import hankel1
import numpy as np
from itertools import product
import matplotlib.pyplot as plt

Ntest=100
U1 = np.random.uniform(size = Ntest)
U2 = np.random.uniform(size = Ntest)
X = np.sqrt(U2) * np.cos(2 * np.pi * U1)
Y = np.sqrt(U2) * np.sin(2 * np.pi * U1)
test_points = np.stack([X,Y]).T

N = 30         #--number of source points
Nc = N*2
R = 1.3
Delta_k = 1.0   #-- used in k-procedure
p_ex = np.array([[10.0 , 10.0]])
boundary_points = np.stack([np.cos(2*np.pi*np.arange(Nc)/Nc) , np.sin(2*np.pi*np.arange(Nc)/Nc)]).T
source_points = np.stack([R*np.cos(2*np.pi*np.arange(N)/N) , R*np.sin(2*np.pi*np.arange(N)/N)]).T

D_test = distance_matrix(test_points , source_points)
D_test_ex = distance_matrix(test_points , p_ex).reshape((Ntest,))
D_ex = distance_matrix(boundary_points , p_ex).reshape((Nc,))
D = distance_matrix(boundary_points , source_points)

Niter = 10000
eps = 0.001
ks = eps*np.arange(1,Niter+1)
values = []
for i in range(1,Niter+1):
   k = i*eps
   w = (1.0j/4)*hankel1(0 , (k+Delta_k)*D_ex)
   G = (1.0j/4)*hankel1(0 , k*D)
   A = np.dot(G.T,G)
   B = np.dot(G.T,w)
   Q = np.linalg.solve(A,-B)
   w_test = (1.0j/4)*hankel1(0 , (k+Delta_k)*D_test_ex)
   G_test = (1.0j/4)*hankel1(0 , k*D_test)
   F = np.linalg.norm( w_test + np.dot(G_test,Q) )
   values.append(F)



plt.plot(ks,np.log(np.array(values)))
plt.show()

厳密な固有関数は
$\psi_{np}(r \cos \theta , r\sin \theta) = j_n(\sqrt{\lambda_{n,p}} r) e^{\pm i \theta}$
という形であるが、n=0の時は、回転不変で、原点からの距離のみに依存する。例えば、n=0で、k=5.5201の場合、計算された固有関数の値の実部を、原点からの距離の関数としてプロットして、厳密な固有関数のそれ(Bessel関数 $J_{0}(kr)$ )と比較したグラフは以下のようになった。赤いほうが、厳密な固有関数だけど、ほぼ完全に一致していると言える。

f:id:m-a-o:20210809122421p:plain

同じ著者による似た内容の論文で、
The method of fundamental solutions for eigenproblems with Laplace andbiharmonic operators
https://scholar.google.com/citations?view_op=view_citation&hl=fr&user=TwTH3u8AAAAJ&citation_for_view=TwTH3u8AAAAJ:d1gkVwhDpl0C
では、重調和方程式も試している。

論文では自由端境界条件は扱われてないが、重調和方程式といえば、クラドニ図形なので(?)、クラドニ図形の問題を、基本解近似法で解いてみることにする。

クラドニ図形のことは、以前に、
コンピュータ以前の数値計算(2)20世紀初頭の計算需要
https://m-a-o.hatenablog.com/entry/2020/06/14/170525
の一番最後にも、書いた。方程式の導出は、Sophie GermainとLagrangeが19世紀初頭に行った。

【重調和方程式の導出について】基本的には、板の微小振動の問題なので、現代的な導出としては、弾性体の波動方程式(特に決まった名前が付いてないけど、Wikipediaなどに、Navier-Cauchy方程式という名前があり、歴史的にも妥当な命名らしい)から始めるべきだろう。Sophie Germainの時には、この方程式は、まだ知られてなかった。この世のマクロな物体は、弦、膜、梁、板のいずれも、本来は、3次元物体なので、どういう近似の下で、方程式を二次元空間や一次元空間上の方程式に簡約できるか明確なのが望ましいと思う。一般的な板理論として、Mindlin理論があり、更に近似を導入した板理論として、Kirchhoff–Love板理論がある。これの特別な場合として、重調和方程式が出る。Mindlin理論から、beam theory(梁理論)を出すこともでき、Timoshenko beam theoryや更に近似を導入すると、Euler-Bernoulli beam theoryを導出できる。これらの導出は、適当な工学の教科書を見れば書いてあるが、計算の見通しは、それほど良くない。19世紀〜20世紀初頭にかけて、これらの一次元や二次元の方程式は、三次元の方程式に比べれば扱いやすいという理由で、詳しく調べられた。

正方形板で単純支持境界条件の場合は、Navierが1823年に解を構成したようだ(未確認)が、方程式が書かれてから約100年近く経った1908~1909年頃、正方形板で自由端境界条件を課した場合の数値解を、Ritzが、いくつか計算した。以前は、このRitzの計算をトレースした。Ritzは、49個の基底で張られる関数空間の中で、近似固有関数を探した。当時は、電子計算機もないので、49x49の行列の計算は、手に負えるものではなかったが、対称性を利用して、小さな行列の計算に帰着させ、数値解を求めたと思われる。

基本解近似法で必要な基本解だけど、上の論文では、 $\Delta^2-k^4=(\Delta+k^2)(\Delta-k^2)$ と可換な2つの微分作用素の積に分解するので(?)、Helmholtz operatorの基本解とmodified Helmholtz operatorの基本解の線形結合を考えている(従って、ソース点の個数Nに対して、決めるべき係数は、2N個ある)。

Ritzの基底は、境界条件を(厳密に)満たすように作られているのだが、基本解近似法では、逆に、重調和方程式は(厳密に)満たすが、境界条件は満たさないので、境界条件を近似的に満たすような解を探すのが仕事になる。

必要なのは自由端境界条件だけど、導出は意外と面倒で、正しい境界条件は、19世紀半ば頃、Kirchhoffが完成させた。円板のように尖った部分がなければKirchhoffの境界条件でいいけど、正方形板の場合は、四隅の角の境界条件がある(これは、Lambが得たそうだ)。まぁ、そんなわけで、境界条件も色々あるけど、天下り的に受け入れることにする。

計算領域を $[-1,+1] \times [-1,+1]$ として、

$x=\pm 1$ では $u_{xx} + \mu \cdot u_{yy} = u_{xxx} + (2-\mu)u_{xyy} = 0$

$y=\pm 1$ では $u_{yy} + \mu \cdot u_{xx} = u_{yyy} + (2-\mu)u_{xxy} - 0$

コーナーでは $u_{xy}=0$

が境界条件。

数値計算された固有値は、

Chladni Figures and the Tacoma Bridge: Motivating PDE Eigenvalue Problems via Vibrating Plates
https://doi.org/10.1137/10081931X

のFigure4.1の数値を信じることにする(固有値は、波数kの4乗であるのを忘れないこと)。

残念ながら、基本解近似法では、全く正しい答えが出なかった。

f:id:m-a-o:20210809121355p:plain

うまくいかない理由は、よく分からない(というか、そもそも一般的に、基本解近似法がうまくいく理由も、よく分からないが)。実装ミスってことは、ないと思うけど、上のDirichlet固有値問題と比べると、こっちは、そもそも境界条件の近似が、あんまりよくないように見える。そもそも、上の例と違って、 $k \to 0$ に近づくにつれて、近似ノルムが発散するっぽいし、ピークも全然シャープじゃない

一応、実装を貼っておく

import sympy
from scipy.special import hankel1,kn
from scipy.spatial import distance_matrix
from numpy import sqrt
import numpy as np
import matplotlib.pyplot as plt

besselk = kn
x,y,x0,y0,kk=sympy.symbols("x,y,x0,y0,kk")
F1 = sympy.hankel1(0 , kk*sympy.sqrt((x-x0)**2 + (y-y0)**2))
F2 = sympy.besselk(0 , kk*sympy.sqrt((x-x0)**2 + (y-y0)**2))
F1_xx = sympy.simplify(F1.diff(x).diff(x))
F1_xxx = sympy.simplify(F1_xx.diff(x))
F1_yy = sympy.simplify(F1.diff(y).diff(y))
F1_yyy = sympy.simplify(F1_yy.diff(y))
F1_xxy = sympy.simplify(F1_xx.diff(y))
F1_xy = sympy.simplify(F1.diff(x).diff(y))
F1_xyy = sympy.simplify(F1_xy.diff(y))

F2_xx = sympy.simplify(F2.diff(x).diff(x))
F2_xxx = sympy.simplify(F2_xx.diff(x))
F2_yy = sympy.simplify(F2.diff(y).diff(y))
F2_yyy = sympy.simplify(F2_yy.diff(y))
F2_xxy = sympy.simplify(F2_xx.diff(y))
F2_xy = sympy.simplify(F2.diff(x).diff(y))
F2_xyy = sympy.simplify(F2_xy.diff(y))

#-- evalf of sympy is too slow
f1_xx = eval("lambda x,y,x0,y0,kk:" + str(F1_xx))
f1_xxx = eval("lambda x,y,x0,y0,kk:" + str(F1_xxx))
f1_yy = eval("lambda x,y,x0,y0,kk:" + str(F1_yy))
f1_yyy = eval("lambda x,y,x0,y0,kk:" + str(F1_yyy))
f1_xxy = eval("lambda x,y,x0,y0,kk:" + str(F1_xxy))
f1_xyy = eval("lambda x,y,x0,y0,kk:" + str(F1_xyy))
f1_xy = eval("lambda x,y,x0,y0,kk:" + str(F1_xy))

f2_xx = eval("lambda x,y,x0,y0,kk:" + str(F2_xx))
f2_xxx = eval("lambda x,y,x0,y0,kk:" + str(F2_xxx))
f2_yy = eval("lambda x,y,x0,y0,kk:" + str(F2_yy))
f2_yyy = eval("lambda x,y,x0,y0,kk:" + str(F2_yyy))
f2_xxy = eval("lambda x,y,x0,y0,kk:" + str(F2_xxy))
f2_xyy = eval("lambda x,y,x0,y0,kk:" + str(F2_xyy))
f2_xy = eval("lambda x,y,x0,y0,kk:" + str(F2_xy))

Ntest=200
U1 = 2*(np.random.uniform(size = Ntest) - 0.5)
U2 = 2*(np.random.uniform(size = Ntest) - 0.5)
test_points = np.stack([U1,U2]).T

mu = 0.225  #--Poisson raio
N = 50      #-- number of source points
Naxis = 31
R = 1.45
Delta_k = 1.0
p_ex = np.array([[10.0 , 10.0]])
source_points = np.stack([R*np.cos(2*np.pi*np.arange(N)/N) , R*np.sin(2*np.pi*np.arange(N)/N)]).T


boundary_points1 = np.stack([np.repeat(1,Naxis-1) , np.linspace(+1 , -1,Naxis+1)[1:-1]]).T  #-- x=1 ,-1<y<1
boundary_points2 = np.stack([np.repeat(-1,Naxis-1) , np.linspace(-1 , +1,Naxis+1)[1:-1]]).T  #-- x=-1 ,-1<y<1
boundary_points3 = np.stack([np.linspace(-1,1 ,Naxis+1)[1:-1] , np.repeat(1,Naxis-1)]).T  #-- y=1 -1<x<1
boundary_points4 = np.stack([np.linspace(1,-1 ,Naxis+1)[1:-1] , np.repeat(-1,Naxis-1)]).T  #-- y=-1 -1<x<1
corners = np.array([[+1,+1],[+1,-1],[-1,+1],[-1,-1]])
boundary_points = np.concatenate([boundary_points1 , boundary_points2 , boundary_points3 , boundary_points4 , corners])
Nx = boundary_points1.shape[0] + boundary_points2.shape[0]
bnd_x = boundary_points[:Nx , :]
bnd_y = boundary_points[Nx:2*Nx , :]
bnd_c = boundary_points[-4: , :]


D_test = distance_matrix(test_points , source_points)
D_test_ex = distance_matrix(test_points , p_ex).reshape((Ntest,))

multi_source_points = np.concatenate([source_points,source_points])
multi_boundary_points = np.concatenate([bnd_x,bnd_x,bnd_y,bnd_y,bnd_c,bnd_c,bnd_c])
Istart = 100
Iend = 400
eps = 0.01
ks = []
values = []

for i in range(Istart , Iend):
   k = eps*i
   D1_1_x = lambda a:f1_xx(a[0],a[1],a[2],a[3],a[4]) + mu*f1_yy(a[0],a[1],a[2],a[3],a[4])
   D2_1_x = lambda a:f2_xx(a[0],a[1],a[2],a[3],a[4]) + mu*f2_yy(a[0],a[1],a[2],a[3],a[4])
   D1_2_x = lambda a:f1_xxx(a[0],a[1],a[2],a[3],a[4]) + (2-mu)*f1_xyy(a[0],a[1],a[2],a[3],a[4])
   D2_2_x = lambda a:f2_xxx(a[0],a[1],a[2],a[3],a[4]) + (2-mu)*f2_xyy(a[0],a[1],a[2],a[3],a[4])
   D1_1_y = lambda a:f1_yy(a[0],a[1],a[2],a[3],a[4]) + mu*f1_xx(a[0],a[1],a[2],a[3],a[4])
   D2_1_y = lambda a:f2_yy(a[0],a[1],a[2],a[3],a[4]) + mu*f2_xx(a[0],a[1],a[2],a[3],a[4])
   D1_2_y = lambda a:f1_yyy(a[0],a[1],a[2],a[3],a[4]) + (2-mu)*f1_xxy(a[0],a[1],a[2],a[3],a[4])
   D2_2_y = lambda a:f2_yyy(a[0],a[1],a[2],a[3],a[4]) + (2-mu)*f2_xxy(a[0],a[1],a[2],a[3],a[4])
   D1_1_c = lambda a:f1_xx(a[0],a[1],a[2],a[3],a[4])
   D2_1_c = lambda a:f2_xx(a[0],a[1],a[2],a[3],a[4])
   D1_2_c = lambda a:f1_yy(a[0],a[1],a[2],a[3],a[4])
   D2_2_c = lambda a:f2_yy(a[0],a[1],a[2],a[3],a[4])
   D1_3_c = lambda a:f1_xy(a[0],a[1],a[2],a[3],a[4])
   D2_3_c = lambda a:f2_xy(a[0],a[1],a[2],a[3],a[4])
   conds = [(D1_1_x,bnd_x),(D1_2_x,bnd_x),(D1_1_y,bnd_y),(D1_2_y,bnd_y),(D1_1_c,bnd_c),(D1_2_c,bnd_c),(D1_3_c,bnd_c),
            (D2_1_x,None),(D2_2_x,None),(D2_1_y,None),(D2_2_y,None),(D2_1_c,None),(D2_2_c,None),(D2_3_c,None)]
   indices = [np.array(np.meshgrid(list(range(0,Nx)) , list(range(0,N)),indexing='ij')).T.reshape((N*Nx,2)) , 
              np.array(np.meshgrid(list(range(Nx,2*Nx)) , list(range(0,N)),indexing='ij')).T.reshape((N*Nx,2)) , 
              np.array(np.meshgrid(list(range(2*Nx,3*Nx)) , list(range(0,N)),indexing='ij')).T.reshape((N*Nx,2)) , 
              np.array(np.meshgrid(list(range(3*Nx,4*Nx)) , list(range(0,N)),indexing='ij')).T.reshape((N*Nx,2)) ,
              np.array(np.meshgrid(list(range(4*Nx,4*Nx+4)) , list(range(0,N)),indexing='ij')).T.reshape((N*4,2)) , 
              np.array(np.meshgrid(list(range(4*Nx+4,4*Nx+8)) , list(range(0,N)),indexing='ij')).T.reshape((N*4,2)),
              np.array(np.meshgrid(list(range(4*Nx+8,4*Nx+12)) , list(range(0,N)),indexing='ij')).T.reshape((N*4,2)),
              np.array(np.meshgrid(list(range(0,Nx)) , list(range(N,2*N)) ,indexing='ij')).T.reshape((N*Nx,2)) , 
              np.array(np.meshgrid(list(range(Nx,2*Nx)), list(range(N,2*N)),indexing='ij')).T.reshape((N*Nx,2)) , 
              np.array(np.meshgrid(list(range(2*Nx,3*Nx)), list(range(N,2*N)),indexing='ij')).T.reshape((N*Nx,2)) , 
              np.array(np.meshgrid(list(range(3*Nx,4*Nx)), list(range(N,2*N)),indexing='ij')).T.reshape((N*Nx,2)),
              np.array(np.meshgrid(list(range(4*Nx,4*Nx+4)), list(range(N,2*N)),indexing='ij')).T.reshape((N*4,2)) ,
              np.array(np.meshgrid(list(range(4*Nx+4,4*Nx+8)), list(range(N,2*N)),indexing='ij')).T.reshape((N*4,2)),
              np.array(np.meshgrid(list(range(4*Nx+8,4*Nx+12)), list(range(N,2*N)),indexing='ij')).T.reshape((N*4,2))]
   Gs = []
   ws = []
   for n,(fun , bnd_points) in enumerate(conds):
       Narg = indices[n].shape[0]
       args = np.concatenate([multi_boundary_points[indices[n][:,0]] ,multi_source_points[indices[n][:,1]] ,np.repeat(k , Narg).reshape((Narg,1))],axis=1)
       Gs.append( np.apply_along_axis(fun , 1 , args) )
       if bnd_points is not None:ws.append( np.apply_along_axis(lambda a:fun( np.array([a[0],a[1],p_ex[0,0],p_ex[0,1] ,k+Delta_k]) ) , 1 , bnd_points) )
   Nb = multi_boundary_points.shape[0]
   Ns = multi_source_points.shape[0]
   G = np.zeros( (Nb , Ns) ,dtype=np.complex )
   I = np.concatenate(indices)
   G[I[:,0],I[:,1]] = np.concatenate(Gs)[:]
   w = np.concatenate(ws)
   A = np.dot(G.T,G)
   B = np.dot(G.T,w)
   Q = np.linalg.solve(A,-B)
   w_test = hankel1(0 , (k+Delta_k)*D_test_ex)
   H1_test = hankel1(0 , k*D_test)
   H2_test = besselk(0 , k*D_test)
   F = np.linalg.norm( w_test + np.dot(H1_test, Q[:N]) + np.dot(H2_test , Q[N:]) )
   ks.append(k)
   values.append(F)



plt.plot(ks , np.array(values))
plt.show()