3次元固有値問題の基本解近似法 -

基本解近似法(a.k.a 代用電荷法)を使ってる文献を見ると、二次元の問題にしか適用されてない。考案されたのが1960年代で、まだ計算機の能力的に2次元問題が中心だったのかもしれないけど、3次元以上でも使えるのか、気になった。

基本解近似法によるクラドニ図形の計算
https://m-a-o.hatenablog.com/entry/2021/08/09/122806

で単位円板に対するDirichlet固有値問題の近似計算をやったので、3次元単位球で同じ計算をやってみる。ヘルムホルツ方程式の基本解が必要だけど、よく知られてる通り、3次元では、
$G_{k}(x) = \dfrac{e^{i k |x| }}{4 \pi |x|}$
と書ける。

厳密な固有関数は、3次元の場合、
$\dfrac{1}{\sqrt{r}}J_{l+\frac{1}{2}}(\sqrt{\lambda_{l,p}}r) Y_{lm}(\omega)$
という形で書ける。 $Y_{lm}$ は球面調和関数。固有値は $\lambda_{l,p}$ で、ベッセル関数 $J_{l+\frac{1}{2}}$ のゼロ点の2乗。

次数が半整数のBessel関数のゼロ点の近似値は、以下のような感じ。10以下の値は、丁度10個ある。

	p=1	p=2	p=3
1/2	3.141592653589793	6.283185307179586	9.42477796076938
3/2	4.4934094579090641	7.72525183693770716	10.904121659428898
5/2	5.76345919689454979	9.09501133047635515
7/2	6.9879320005005199
9/2	8.182561452571242
11/2	9.355812111042746

実装は、2次元の時と殆ど同じで、以下のような感じ。

import numpy as np
from scipy.spatial import distance_matrix
import matplotlib.pyplot as plt

def make_points_on_sphere(N):
   U1 = np.random.uniform(size = N)
   U2 = np.random.uniform(size = N)
   X = np.sin(2 * np.pi * U1) * np.cos(2 * np.pi * U2)
   Y = np.cos(2 * np.pi * U1) * np.cos(2 * np.pi * U2)
   Z = np.sin(2 * np.pi * U2)
   return np.stack([X,Y,Z]).T


def make_points_in_ball(N):
   U0 = np.random.uniform(size = N)
   U1 = np.random.uniform(size = N)
   U2 = np.random.uniform(size = N)
   X = U0 * np.sin(2 * np.pi * U1) * np.cos(2 * np.pi * U2)
   Y = U0 * np.cos(2 * np.pi * U1) * np.cos(2 * np.pi * U2)
   Z = U0 * np.sin(2 * np.pi * U2)
   return np.stack([X,Y,Z]).T


"""
Ntest : number of test points (distributed in the ball)
Nsrc : number of source points (exist on the outside of the ball)
Nbnd : number of boundary points (distributed on the sphere)
"""
def main(Ntest , Nsrc , Nbnd):
   R = 1.6
   Delta_k = 1.0
   p_ex = np.array([[5.0, 5.0, 5.0]])
   test_points = make_points_in_ball(Ntest)
   boundary_points = make_points_on_sphere(Nbnd)
   source_points = R * make_points_on_sphere(Nsrc)
   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((Nbnd,))
   D = distance_matrix(boundary_points , source_points)
   istart = 301
   iend = 800
   eps = 0.01
   ks = []
   values = []
   for i in range(istart , iend):
       k = i * eps
       ks.append( k )
       w = np.exp(1.0j * (k+Delta_k) * D_ex) / D_ex
       G = np.exp(1.0j * k * D) / D
       A = np.dot(G.T , G)
       B = np.dot(G.T , w)
       Q = np.linalg.solve(A , -B)
       w_test = np.exp(1.0j * (k+Delta_k) * D_test_ex) / D_test_ex
       G_test = np.exp(1.0j * k * D_test) / D_test
       F = np.linalg.norm( w_test + np.dot(G_test,Q) )
       values.append(F)
   plt.plot(np.array(ks) , np.log(np.array(values)))
   plt.show()


if __name__=="__main__":
   np.random.seed(12345)
   main(3000,1500,2000)

パラメータは、l=0,p=1,2,3の時の固有関数の値が再現できるかどうかで決めた。点の数が、二次元の時より、かなり多いけど、これくらいないと、厳密解の値を再現できなかった。多分、二次元の時と違って、点をランダムに取ってるのも影響してるっぽい

また、固有値が大きいほど、多くの点を必要とするっぽい。以下は、l=0,p=3で計算した固有関数の(300点における)実部の値を適当に定数倍したものと、厳密解を重ねたもの。l=0の時は、球対称なので、固有関数の値は、原典からの距離rのみに依存するはずで、横軸が、原点からの距離r

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