二重周期的ポテンシャル問題に対する代用電荷法
緒方 秀教
∗ ∗電気通信大学大学院情報理工学研究科
情報・ネットワーク工学専攻
概要. 二重周期的二次元ポテンシャル問題に対する代用電荷法を，特にポテンシャル流 問題に対して提案する．この問題の解は二重周期関数を含むので，従来の代用電荷法で近 似解を求めるのは難しい．そこで，周期的ポテンシャルをWeierstrass楕円関数を用いて 構成し，その一次結合を用いて解を近似する．この方法は従来の代用電荷法の利点を受け 継ぐ上に，解が持つべき周期性を持つ．数値例により本方法の有効性が示される．Method of Fundamental Solutions for Doubly-Periodic
Potential Problems
Hidenori Ogata
∗∗
Department of Computer and Network Engineering,
Graduate School of Informatics and Engineering,
The University of Electro-Communications
Abstract. We propose a method of fundamental solutions for two-dimensional poten-tial problems, especially potenpoten-tial flow problems, with double periodicity. It is difficult to approximate the solution of these problems by the conventional method because the solution involves a periodic function. In our method, the solution is approximated by a linear combination of the periodic potentials, which is constructed using the Weierstrass elliptic functions. Our method inherits the advantanges of the conventional method, and, furthermore, it has the periodicity of the exact solution. Numerical examples show the effectiveness of our method.
1.
はじめに
代用電荷法[13]はポテンシャル問題に対する数値解法である．この方法は，ポテンシャ ル問題の近似解を問題領域外部に置かれた点電荷ポテンシャルの重ね合わせで近似するも のである．二次元領域D ⊂ R2 における二次元ポテンシャル問題*1 { △u = 0 in D u = f on ∂D *1 本論文では，整数全体の集合をZ，実数全体の集合をR，複素数全体の集合（複素平面）をCで表す．の場合について詳しく述べると，二次元Euclid平面R2 と複素平面Cを同一視して，解で あるポテンシャルuを次のように近似する*1： u(z) ≃ uN(z) = Q0 − 1 2π N ∑ j=1 Qjlog |z − ζj|, (1.1) ここで，ζ1, . . . , ζN は問題領域D 外部にとった点で電荷点とよばれ，Q0,Q1, . . . ,QN ∈ R は N ∑ j=1 Qj =0 (1.2) を満たす未知係数であり電荷とよばれる．近似解uN は問題領域D においてLaplace方程 式△uN =0を厳密に満たす．境界条件については，拘束条件を課して近似的に満たすよう にする．すなわち，境界∂D の境界上に拘束点とよばれる点z1, . . . ,zN をとって，条件 uN(zi) = f (zi), i = 1, . . . , N (1.3) を課す．方程式(1.3)を書き下すと Q0− 1 2π N ∑ j=1 Qjlog |zi− ζj| = f (zi), i = 1, . . . , N (1.4) となり，(1.2)と合わせて電荷Qj, j = 0, 1, . . . , N に対する連立一次方程式をなす．連立一 次方程式(1.2), (1.4)を解いて未知の電荷Qj を決定し，近似ポテンシャルuN を得る． このように，代用電荷法は素朴な発想に基づく方法であり，プログラミングは簡単で計 算量は多くないが，ある条件のもとでは指数関数的収束する[10]という優れた面を持つ． 代用電荷法は当初電気力学の問題に適用されたが[22]，その後，波動問題，弾性問題，重 調和問題[4]など，科学技術計算の様々な分野で用いられるようになった． 代用電荷法は複素解析関数の近似にも用いられる．ある領域D における複素解析関数 f (z)を然るべき境界条件のもとで求めたいとする．この時，実部Re f は領域D における 調和関数であるから，代用電荷法により Re f (z) ≃ Q0 − 1 2π N ∑ j=1 Qjlog |z − ζj| と近似できる．電荷 Qj は境界条件に対応する拘束条件を課すことにより決定される．す ると，虚部Im f は実部Re f の共役調和関数であるから， Im f (z) ≃ Q′ 0− 2π1 N ∑ j=1 Qjarg(z − ζj) *1 ここに示す近似解は，[14]の不変スキームによるものである．
の場合について詳しく述べると，二次元Euclid平面R2 と複素平面Cを同一視して，解で あるポテンシャルuを次のように近似する*1： u(z) ≃ uN(z) = Q0 − 1 2π N ∑ j=1 Qjlog |z − ζj|, (1.1) ここで，ζ1, . . . , ζN は問題領域D 外部にとった点で電荷点とよばれ，Q0,Q1, . . . ,QN ∈ R は N ∑ j=1 Qj =0 (1.2) を満たす未知係数であり電荷とよばれる．近似解uN は問題領域D においてLaplace方程 式△uN =0を厳密に満たす．境界条件については，拘束条件を課して近似的に満たすよう にする．すなわち，境界∂D の境界上に拘束点とよばれる点z1, . . . ,zN をとって，条件 uN(zi) = f (zi), i = 1, . . . , N (1.3) を課す．方程式(1.3)を書き下すと Q0− 1 2π N ∑ j=1 Qjlog |zi− ζj| = f (zi), i = 1, . . . , N (1.4) となり，(1.2)と合わせて電荷Qj, j = 0, 1, . . . , N に対する連立一次方程式をなす．連立一 次方程式(1.2), (1.4)を解いて未知の電荷 Qj を決定し，近似ポテンシャルuN を得る． このように，代用電荷法は素朴な発想に基づく方法であり，プログラミングは簡単で計 算量は多くないが，ある条件のもとでは指数関数的収束する[10]という優れた面を持つ． 代用電荷法は当初電気力学の問題に適用されたが[22]，その後，波動問題，弾性問題，重 調和問題[4]など，科学技術計算の様々な分野で用いられるようになった． 代用電荷法は複素解析関数の近似にも用いられる．ある領域D における複素解析関数 f (z)を然るべき境界条件のもとで求めたいとする．この時，実部Re f は領域D における 調和関数であるから，代用電荷法により Re f (z) ≃ Q0− 1 2π N ∑ j=1 Qjlog |z − ζj| と近似できる．電荷 Qj は境界条件に対応する拘束条件を課すことにより決定される．す ると，虚部Im f は実部Re f の共役調和関数であるから， Im f (z) ≃ Q′ 0− 2π1 N ∑ j=1 Qjarg(z − ζj) *1 ここに示す近似解は，[14]の不変スキームによるものである． と近似される．この2式を合わせて，複素解析関数 f (z)に対する代用電荷法近似 f (z) ≃ fN(z) = Q0+iQ′0− 1 2π N ∑ j=1 Qjlog(z − ζj) (1.5) を得る．この考えに基づき，天野らは代用電荷法による数値等角写像の方法を提案してい る[2, 3]． 本論文では，二重周期性をもつ二次元ポテンシャル問題を代用電荷法で解くことを考え る．具体的には，Fig.1のような二重周期性を持つ領域を流れる二次元ポテンシャル流を 求める問題を考える．一般に二次元ポテンシャル流は複素速度ポテンシャルという関数で 記述される．これは流れの領域における解析関数 f (z)で然るべき境界条件を満たすもの であり，流速場v = (u, v)を f′(z) = u − ivとその導関数により与える．したがって，ポテ ンシャル流を数値的に求めるにはこの複素速度ポテンシャルの近似関数を代用電荷法で求 めればよいのだが，ここでひとつ問題がある．それは，解が周期関数であり，従来の代用 電荷法では近似解を求めるのが難しいということである．二重周期的領域におけるポテン シャル流の流速場vは二重周期関数となり，したがってそれを与える複素速度ポテンシャ ル f (z)も当然二重周期関数を含む．ところが，そのような解析関数を従来の代用電荷法 による(1.5)の形では近似することは難しい． したがって，今考えている二重周期的ポテンシャル流の問題に対しては，別の形の代用 電荷法近似を考えなければならない．そこで，本論文では二重周期的ポテンシャルの一次 結合による近似解を考えることを提案する．従来の代用電荷法近似 (1.1)あるいは(1.5) では1点にのみソースを持つ点電荷ポテンシャルの一次結合で近似解を与えているが，代 わりに，二重周期的に並ぶ同じ強さの点電荷列によるポテンシャルを考え，それらの一次 結合により周期的ポテンシャル流の近似解を与えるのである．二重周期的ポテンシャルは 具体的には，Weierstrass楕円関数を用いて構成する．これにより，従来の代用電荷法の利 点を引き継ぎ，さらに，問題から要請される周期性を近似解が満たす代用電荷法を構成す ることができると期待される． 著者はすでに周期性をもつ問題に対し，同様の発想に基づく代用電荷法を提案してい る．一重周期的二次元ポテンシャル問題，具体的には周期的複素領域に対する数値等角 写像の問題に対し，周期的ポテンシャルの一次結合を用いた代用電荷法を提案してい る[21]．ここでは，一重周期的ポテンシャルを三角関数を用いて構成している．また，周 期的 Stokes流問題に対しても，著者は同様の発想による代用電荷法，すなわち，周期的 に並ぶ点荷重列により生じる流れ（周期的Stokeslet）の一次結合により近似解を与える方 法を提案している[16–19]．なお，流体力学の分野では周期性をもつ問題の研究は多く行 われており，例えば[5, 25]では周期的Stokes流問題に対する積分方程式による解法が研 究されている．応用的研究として，[8]では周期的Stokes流の解法の応用として，生物の 鞭毛による流れの研究が行われている． 本論文の構成は次のとおりである．第2節では，二次元ポテンシャル流の理論を大ま
かに述べて，本論文で扱う二重周期的二次元ポテンシャル流の問題を数学的に定式化す る．第 3 節では，前節で数学的に与えた二重周期的二次元ポテンシャル流問題に対し， Weierstrass 楕円関数を用いて構成した二重周期的ポテンシャルの一次結合による代用電 荷法を提案する．第4節では，数値例により本論文の代用電荷法の有効性を検証する．第 5節では本論文の総括を行い，今後の課題について触れる．
2.
二重周期的ポテンシャル流問題
はじめに二重周期的ポテンシャル流の流域は，数学的には次の複素領域D で与えられ る（Fig.1参照）． D = C\ ∪ m,n∈Z Dmn, (2.1) ここで，D00は有界単連結領域，ω1, ω2 ∈ C \ { 0 }をIm (ω2/ω1) > 0なる定数として Dmn = { z + mω1+nω2 | z ∈ D00} , m, n ∈ Z (2.2) とする．ω1, ω2 は二重周期的領域の基本周期を表し，Dmn, m, n ∈ Zは二重周期的に並ぶ 障害物列を表す．Fig. 1. The domain D with a doubly-periodic array of obstacles Dmn, m, n ∈ Z.
式(2.1)の二重周期的領域Dにおける二次元ポテンシャル流を求める．ところで，複素 領域D におけるポテンシャル流は，複素速度ポテンシャルという複素関数 f (z)によって 特徴づけられる（二次元ポテンシャル流の詳細については[7, 12]などを参照）．これはD における複素解析関数で，その実部（速度ポテンシャル）Φ = Re f および虚部（流れ関 数）Ψ =Im f により流速場v = (u, v)を u = ∂Φ ∂x = ∂Ψ ∂y, v = ∂Φ ∂y = − ∂Ψ ∂x で与える量である．この式は，(Φ, Ψ)に対するCauchy-Riemann関係式に他ならない．し
かに述べて，本論文で扱う二重周期的二次元ポテンシャル流の問題を数学的に定式化す る．第 3 節では，前節で数学的に与えた二重周期的二次元ポテンシャル流問題に対し， Weierstrass 楕円関数を用いて構成した二重周期的ポテンシャルの一次結合による代用電 荷法を提案する．第4節では，数値例により本論文の代用電荷法の有効性を検証する．第 5節では本論文の総括を行い，今後の課題について触れる．
2.
二重周期的ポテンシャル流問題
はじめに二重周期的ポテンシャル流の流域は，数学的には次の複素領域D で与えられ る（Fig.1参照）． D = C\ ∪ m,n∈Z Dmn, (2.1) ここで，D00は有界単連結領域，ω1, ω2 ∈ C \ { 0 }をIm (ω2/ω1) > 0なる定数として Dmn = { z + mω1+nω2 | z ∈ D00} , m, n ∈ Z (2.2) とする．ω1, ω2 は二重周期的領域の基本周期を表し，Dmn, m, n ∈ Zは二重周期的に並ぶ 障害物列を表す．Fig. 1. The domain D with a doubly-periodic array of obstacles Dmn, m, n ∈ Z.
式(2.1)の二重周期的領域D における二次元ポテンシャル流を求める．ところで，複素 領域D におけるポテンシャル流は，複素速度ポテンシャルという複素関数 f (z)によって 特徴づけられる（二次元ポテンシャル流の詳細については[7, 12]などを参照）．これはD における複素解析関数で，その実部（速度ポテンシャル）Φ = Re f および虚部（流れ関 数）Ψ =Im f により流速場v = (u, v)を u = ∂Φ ∂x = ∂Ψ ∂y, v = ∂Φ ∂y = − ∂Ψ ∂x で与える量である．この式は，(Φ, Ψ)に対するCauchy-Riemann関係式に他ならない．し たがって，速度場v = (u, v)は f (z)の導関数を用いて f′(z) = u − iv で与えられる．そして，流れ関数Ψの等高線に沿って， 0 = dΨ = ∂Ψ ∂xdx + ∂Ψ ∂ydy = −vdx + udy が成り立つので，Ψ =Im f の等高線は流れの流線を与える．流体は流域の境界∂D に沿っ て流れるので，その上でΨ =Im f は一定である： Im f = constant on ∂D. (2.3) ゆえに，二重周期的領域D におけるポテンシャル流を求めるには，D における複素解析 関数 f (z)で境界条件(2.3)を満たすものを求めればよい． そして，今考えている問題では，流れの領域D は周期的構造が二次元的に無限遠まで 続く平面領域なので，解を一意的に与えるための通常の境界条件を設定することができな い．そこで代わりに，流れの平均流速が⟨v⟩ = (U, 0)であるという条件を課すことにする （U > 0は平均流速の大きさを表す定数，平均流速の方向に複素平面の実軸を取ることに する）．すなわち，次の条件を課す： ⟨ f′⟩ = |D1 0| D0 f′(z)dxdy = U, (2.4) ここで，D0は D0 ={ aω1+bω2 | 0 ≤ a, b ≤ 1 } \ ∪ m,n∈Z Dmn, すなわち，周期平行四辺形から障害物列を除いた部分，|D0|はその面積である． 以上により，二重周期的ポテンシャル問題は，領域D における複素解析関数（複素速 度ポテンシャル）f (z)で条件(2.3), (2.4)を満たすものを求める問題に帰着された．
3.
二重周期的ポテンシャル流問題に対する代用電荷法
前節で述べたように，二重周期的領域D におけるポテンシャル流を求めるには，複素 速度ポテンシャルとよばれるD における複素解析関数 f (z)で条件(2.3), (2.4)を満たすも のを求めればよい．この複素速度ポテンシャル f (z)を近似的に求めるのに代用電荷法を 使うのであるが，今の問題の場合，従来の代用電荷法，すなわち，式(1.5)の形を用いて ポテンシャル f (z)を近似するのは難しい．なぜなら，今考えているポテンシャル流は明 らかに二重周期性を持つので，それを与える複素速度ポテンシャルは二重周期関数を含む はずであるが，従来の代用電荷法による近似式 (1.5)では二重周期関数を含む解析関数を 近似するのは難しいからである．そこで，従来の代用電荷法の代わりに，次のような複素速度ポテンシャルの近似を考 える． f (z) ≃ fN(z) = Uz − i 2π N ∑ j=1 Qj{log σ(z − ζj) − vjz}. (3.1) ここで，ζ1, . . . , ζN は障害物D00内部に与えられた点（電荷点），Q1, . . . ,QN ∈ Rは N ∑ j=1 Qj =0 (3.2) を満たす未知係数（電荷）である．σ(z)はWeierstrassのシグマ関数*1 σ(z) = z ∏ ω∈Λ\{ 0 } ( 1 − ωz)exp( z ω + z2 2ω2 ) , (3.3) ここでΛは周期格子 Λ ={ mω1 +nω2 | m, n ∈ Z } であり，vj, j = 1, . . . , NはWeierstrassのゼータ関数 ζ(z) = d dzlog σ(z) = 1 z + ∑ ω∈Λ\{ 0 } ( 1 z − ω + 1 ω + z ω2 ) を用いて vj = 1 |D0| D0 ζ(z − ζj)dxdy (3.4) で与えられる．式(3.1)右辺の各項に現れる−(1/(2π))Qjlog σ(z−ζj)は各格子点z = ζj+ω, ω ∈ Λ近傍では −2π1 Qjlog σ(z − ζj− ω) ∼ −2π1 Qjlog(z − ζj− ω) のように振る舞うから，格子点z = ζj+ ω, ω ∈ Λに電荷Qj の二重周期列による周期的ポ テンシャルを表している． この近似速度ポテンシャル fN(z)は，シグマ関数(3.3)が格子点ω ∈ Λに1位の零点を 持つ整関数であるから，D における解析関数である．そして，fN(z)は擬周期性 fN(z + ωk) = fN(z) + constant, k = 1, 2 (3.5) を満たす．実際，シグマ関数の擬周期性 σ(z + ωk) = −eηk(z+ωk/2)σ(z), ηk =2ζ( ω2k ) , k = 1, 2 (3.6) *1 楕円関数論の詳細については，例えば[24]を参照すること．
そこで，従来の代用電荷法の代わりに，次のような複素速度ポテンシャルの近似を考 える． f (z) ≃ fN(z) = Uz − i 2π N ∑ j=1 Qj{log σ(z − ζj) − vjz}. (3.1) ここで，ζ1, . . . , ζN は障害物D00内部に与えられた点（電荷点），Q1, . . . ,QN ∈ Rは N ∑ j=1 Qj =0 (3.2) を満たす未知係数（電荷）である．σ(z)はWeierstrassのシグマ関数*1 σ(z) = z ∏ ω∈Λ\{ 0 } ( 1 − ωz)exp( z ω + z2 2ω2 ) , (3.3) ここでΛは周期格子 Λ ={ mω1+nω2 | m, n ∈ Z } であり，vj, j = 1, . . . , NはWeierstrassのゼータ関数 ζ(z) = d dzlog σ(z) = 1 z + ∑ ω∈Λ\{ 0 } ( 1 z − ω + 1 ω + z ω2 ) を用いて vj = 1 |D0| D0 ζ(z − ζj)dxdy (3.4) で与えられる．式(3.1)右辺の各項に現れる−(1/(2π))Qjlog σ(z−ζj)は各格子点z = ζj+ω, ω ∈ Λ近傍では −2π1 Qjlog σ(z − ζj− ω) ∼ −2π1 Qjlog(z − ζj− ω) のように振る舞うから，格子点z = ζj+ ω, ω ∈ Λに電荷Qj の二重周期列による周期的ポ テンシャルを表している． この近似速度ポテンシャル fN(z)は，シグマ関数(3.3)が格子点ω ∈ Λに1位の零点を 持つ整関数であるから，D における解析関数である．そして，fN(z)は擬周期性 fN(z + ωk) = fN(z) + constant, k = 1, 2 (3.5) を満たす．実際，シグマ関数の擬周期性 σ(z + ωk) = −eηk(z+ωk/2)σ(z), ηk =2ζ( ω2k ) , k = 1, 2 (3.6) *1 楕円関数論の詳細については，例えば[24]を参照すること． を用いれば， fN(z + ωk) =U(z + ωk) − i 2π N ∑ j=1 Qj { log(−1) + ηk ( z − ζj+ ωk 2 ) +log σ(z − ζj) − vj(z + ωk) } = fN(z) + i 2π N ∑ j=1 Qj(ηkζj+vjωk) = fN(z) + constant, k = 1, 2 を得る．ただし，2番目の等号で(3.2)を用いた．したがって，近似速度場 f′ N(z) = uN−ivN は期待されたとおり二重周期性 f′ N(z + ωk) = fN′(z), k = 1, 2 を満たす，すなわち，楕円関数である．そして，fN(z)の定義(3.1)からわかる通り，近似 速度場 f′ N(z)は平均流速の仮定(2.4)，すなわち， ⟨ fN′⟩ = |D1 0| D0 f′ N(z)dxdy = U (3.7) を満たす．残った境界条件(2.3)については，fN(z)に拘束条件を課して近似的に満たすこ とにする．すなわち，拘束点z1, . . . ,zN ∈ ∂D00をとり， Im fN(zi) = C, i = 1, . . . , N (3.8) （C は未知の実定数）を満たすように未知の電荷 Qj を定める．拘束条件(3.8) を書き直 すと， −C − 2π1 N ∑ j=1 Qj{log |σ(zi − ζj)| − Re (vjzi)} =−U(Im zi), i = 1, . . . , N (3.9) となり，(3.2) と合わせてQj, Cに対する連立一次方程式をなす．この連立一次方程式を 解くことにより未知の電荷Qj が定まり，近似速度ポテンシャル fN(z)が得られる．∂D00 上で拘束条件を課せば，fN(z)の擬周期性(3.5)より他の障害物上でも Im fN(zi+mω1+nω2) = constant, i = 1, . . . , N; m, n ∈ Z が成り立ち，境界条件が近似的に満たされる． なお，シグマ関数σ(z)およびゼータ関数ζ(z)の数値計算は次のようにして行う．シグ マ関数はテータ関数を用いて次のように表される（[24]第3.13節を参照）． σ(z) = ω1 ϑ′1 exp ( η1z2 2ω1 ) ϑ1( z ω1 �� �� � τ ) , τ = ω2 ω1, (3.10)
ここでϑ1(u|τ)はテータ関数 ϑ1(v|τ) (3.11) =2 ∞ ∑ n=0 (−1)nq(n+1/2)2sin((2n + 1)πv), =2q1/4sin πv ∞ ∏ n=1 (1 − q2n)(1 − 2q2ncos 2πv + q4n), q = eiπτ ( Im τ > 0 ), ϑ′1 = ϑ′1(0|τ)である．ゼータ関数ζ(z)は(3.10)の対数微分を計算することにより得られ る．(3.6)中で定義されたη1を η1 = π 2 3ω1 1 − 24 ∞ ∑ n=1 nq2n 1 − q2n で計算すれば（[23]第37節を参照*1），Im τ = Im (ω2/ω1) > 0より|q| < 1であるから， 式(3.11)における無限和あるいは無限積は速く収束するので，ϑ1(v|τ)は容易に計算する ことができる．
4.
数値例
本節では，いくつかの数値例により本論文で提案した代用電荷法の有効性を検証する． すべての数値計算はC++でプログラミングして，倍精度演算により行った．4.1
円柱列を過ぎる流れ
第一の例は，円柱列を過ぎるポテンシャル流である．すなわち，(2.1), (2.2)において D00 ={ z ∈ C | |z| < R } （R > 0は定数）ととった場合である．この問題に対し，いくつかの基本周期の組(ω1, ω2) について本論文の代用電荷法を適用してポテンシャル流を求めた．ただし，代用電荷法の 拘束点zi，電荷点ζj は次のようにとった． zj =R exp ( i2π( j − 1) N ) , ζj =qR exp ( i2π( j − 1) N ) , j = 1, . . . , N; 0 < q < 1. (4.1) 式 (3.4) の係数 vj を求めるには D0 上の二重積分を計算する必要があるが，これは Mersenne twister [11]により生成された乱数を用いて得た100万個の分点を用いてモンテ カルロ法により計算した．Fig.2 にいくつかの基本周期の組(ω1, ω2)に対して計算したポ テンシャル流の流線を示す．なお，図に示した数値計算例では，分点数 N = 32，(4.1)で 電荷点ζj の配置に用いたパラメータqの値はq = 0.4ととっている． *1 [23]では，本論文でω1, ω2, η1, η2 と記している量をそれぞれ2ω1,2ω3,2η1,2η3と記していることに注 意．ここでϑ1(u|τ)はテータ関数 ϑ1(v|τ) (3.11) =2 ∞ ∑ n=0 (−1)nq(n+1/2)2 sin((2n + 1)πv), =2q1/4sin πv ∞ ∏ n=1 (1 − q2n)(1 − 2q2ncos 2πv + q4n), q = eiπτ ( Im τ > 0 ), ϑ′1 = ϑ′1(0|τ)である．ゼータ関数ζ(z)は(3.10)の対数微分を計算することにより得られ る．(3.6)中で定義されたη1を η1 = π 2 3ω1 1 − 24 ∞ ∑ n=1 nq2n 1 − q2n で計算すれば（[23]第37節を参照*1 ），Im τ = Im (ω2/ω1) > 0より|q| < 1であるから， 式(3.11)における無限和あるいは無限積は速く収束するので，ϑ1(v|τ)は容易に計算する ことができる．
4.
数値例
本節では，いくつかの数値例により本論文で提案した代用電荷法の有効性を検証する． すべての数値計算はC++でプログラミングして，倍精度演算により行った．4.1
円柱列を過ぎる流れ
第一の例は，円柱列を過ぎるポテンシャル流である．すなわち，(2.1), (2.2)において D00 ={ z ∈ C | |z| < R } （R > 0は定数）ととった場合である．この問題に対し，いくつかの基本周期の組(ω1, ω2) について本論文の代用電荷法を適用してポテンシャル流を求めた．ただし，代用電荷法の 拘束点zi，電荷点ζj は次のようにとった． zj = R exp ( i2π( j − 1) N ) , ζj =qR exp ( i2π( j − 1) N ) , j = 1, . . . , N; 0 < q < 1. (4.1) 式 (3.4) の係数 vj を求めるには D0 上の二重積分を計算する必要があるが，これは Mersenne twister [11]により生成された乱数を用いて得た100万個の分点を用いてモンテ カルロ法により計算した．Fig.2 にいくつかの基本周期の組(ω1, ω2)に対して計算したポ テンシャル流の流線を示す．なお，図に示した数値計算例では，分点数 N = 32，(4.1)で 電荷点ζj の配置に用いたパラメータqの値はq = 0.4ととっている． *1 [23]では，本論文でω1, ω2, η1, η2 と記している量をそれぞれ2ω1,2ω3,2η1,2η3と記していることに注 意． 0 2 4 6 8 0 2 4 6 8 Im(z) Re(z) 0 2 4 6 8 0 2 4 6 8 Im(z) Re(z)(ω1, ω2) = (4R, 4Ri) (ω1, ω2) = (4R, 4Reiπ/3)
0 2 4 6 8 0 2 4 6 8 Im(z) Re(z) 0 2 4 6 8 0 2 4 6 8 Im(z) Re(z)
(ω1, ω2) = (4Reiπ/6,4Ri) (ω1, ω2) = (4R, 4Reiπ/4)
Fig. 2. The streamlines of the potential flows past a doubly-periodic array of cylinders.
本論文の方法の精度を調べるため，近似ポテンシャル fN(z)の境界条件(2.3)に関する 誤差 ϵN = UR1 max z∈∂D00|Im fN(z) − C| (4.2) を求めた*1．ここで，C は連立一次方程式 (3.2), (3.9)を解いて得られる定数C である． Fig.2の数値例に対し誤差ϵN を求め，それらの分点数N に対する変化をFig.3に示した． 図より，いずれの例でも，分点数N を大きくするにつれ誤差ϵN は指数関数的減衰するこ とがわかる．二次元円板内部における通常のポテンシャル問題に対し，拘束点・電荷点を 同心円上に等間隔にとった従来の代用電荷法を適用すると，誤差が指数関数的減衰するこ とが理論的に示されているが[10]，周期的円板列外部ポテンシャル問題に対する本論文 の代用電荷法でも同様の結果が成立することになる．各数値例について誤差ϵN の減衰率 を，電荷点 ζj を式(4.1) で与えるときのパラメータq の値を変えながら，ソフトウェア *1 ∂D00上の最大値は実際には，境界∂D00上に等間隔に1024個の点をとり，それらの点における最大値を 求めている．
gnuplotのfit コマンドを用いて最小二乗法により求めた．その結果をTable 1に示す． なお，ここで求めた誤差の減衰はFig.3中のグラフにも破線で示した．表より，基本周期 ω1, ω2を変えても誤差ϵN の減衰率はさほど変わらず，大まかに ϵN = { O(0.5N) q = 0.4 O(qN) q = 0.5, 0.6, 0.7 に従うことがわかる． 1x10-16 1x10-14 1x10-12 1x10-10 1x10-8 1x10-6 0.0001 0.01 1 0 10 20 30 40 50 60 70 error N q=0.4 q=0.5 q=0.6 q=0.7 1x10-16 1x10-14 1x10-12 1x10-10 1x10-8 1x10-6 0.0001 0.01 1 0 10 20 30 40 50 60 70 error N q=0.4 q=0.5 q=0.6 q=0.7
(ω1, ω2) = (4R, 4Ri) (ω1, ω2) = (4R, 4Reiπ/3)
1x10-16 1x10-14 1x10-12 1x10-10 1x10-8 1x10-6 0.0001 0.01 1 0 10 20 30 40 50 60 70 error N q=0.4 q=0.5 q=0.6 q=0.7 1x10-16 1x10-14 1x10-12 1x10-10 1x10-8 1x10-6 0.0001 0.01 1 0 10 20 30 40 50 60 70 error N q=0.4 q=0.5 q=0.6 q=0.7
(ω1, ω2) = (4Reiπ/6,4Ri) (ω1, ω2) = (4R, 4Reiπ/4)
Fig. 3. The error estimates ϵN of the method of fundamental solutions applied to the
examples shown in Fig.2. In the figures, “error” means ϵN, and q is the parameter q used
for positioning the charge points ζjby (4.1).
なお，流速平均の条件(3.7) は，近似ポテンシャル fN の式(3.1) においてvj を正確に 計算すれば厳密に満たされるが，D0 を与える式(3.4)においてD0 上の二重積分をモン テカルロ法により計算しているので，厳密には満たされない．そこで，実際の平均流速 ⟨vN⟩ = (⟨uN⟩, ⟨vN⟩), ⟨uN⟩ − i⟨vN⟩ = ⟨ fN′⟩ = |D1 0| D0
gnuplotのfit コマンドを用いて最小二乗法により求めた．その結果をTable 1に示す． なお，ここで求めた誤差の減衰はFig.3中のグラフにも破線で示した．表より，基本周期 ω1, ω2を変えても誤差ϵN の減衰率はさほど変わらず，大まかに ϵN = { O(0.5N) q = 0.4 O(qN) q = 0.5, 0.6, 0.7 に従うことがわかる． 1x10-16 1x10-14 1x10-12 1x10-10 1x10-8 1x10-6 0.0001 0.01 1 0 10 20 30 40 50 60 70 error N q=0.4 q=0.5 q=0.6 q=0.7 1x10-16 1x10-14 1x10-12 1x10-10 1x10-8 1x10-6 0.0001 0.01 1 0 10 20 30 40 50 60 70 error N q=0.4 q=0.5 q=0.6 q=0.7
(ω1, ω2) = (4R, 4Ri) (ω1, ω2) = (4R, 4Reiπ/3)
1x10-16 1x10-14 1x10-12 1x10-10 1x10-8 1x10-6 0.0001 0.01 1 0 10 20 30 40 50 60 70 error N q=0.4 q=0.5 q=0.6 q=0.7 1x10-16 1x10-14 1x10-12 1x10-10 1x10-8 1x10-6 0.0001 0.01 1 0 10 20 30 40 50 60 70 error N q=0.4 q=0.5 q=0.6 q=0.7
(ω1, ω2) = (4Reiπ/6,4Ri) (ω1, ω2) = (4R, 4Reiπ/4)
Fig. 3. The error estimates ϵN of the method of fundamental solutions applied to the
examples shown in Fig.2. In the figures, “error” means ϵN, and q is the parameter q used
for positioning the charge points ζjby (4.1).
なお，流速平均の条件(3.7)は，近似ポテンシャル fN の式(3.1) においてvj を正確に 計算すれば厳密に満たされるが，D0 を与える式(3.4)においてD0 上の二重積分をモン テカルロ法により計算しているので，厳密には満たされない．そこで，実際の平均流速 ⟨vN⟩ = (⟨uN⟩, ⟨vN⟩), ⟨uN⟩ − i⟨vN⟩ = ⟨ fN′⟩ = |D1 0| D0
f′(z)dxdy ≃ ⟨u − iv⟩
Table 1. The decay rates of the errors ϵN of the presented method applied to the examples
shown in Fig.2, where q is the parameter used for positioning the charge points ζjby (4.1).
(ω1, ω2) q
0.4 0.5 0.6 0.7
(4R, 4Ri) O(0.51N) O(0.51N) O(0.58N) O(0.67N) (4R, 4Reiπ/3) O(0.50N) O(0.50N) O(0.58N) O(0.67N) (4Reiπ/6,4Ri) O(0.49N) O(0.50N) O(0.58N) O(0.67N) (4R, 4Reiπ/4) O(0.54N) O(0.56N) O(0.60N) O(0.67N)
をFig.2の数値例についてモンテカルロ法積分により計算した．その結果は ⟨vN⟩ U = (1.00024, −5.4 × 10−5) (ω1, ω2) = (4R, 4Ri) (0.99981, −1.1 × 10−3) (ω1, ω2) = (4R, 4Reiπ/3)
(1.00030, 9.8 × 10−5) (ω1, ω2) = (4Reiπ/6,4Ri)
(0.99991, −1.4 × 10−4) (ω1, ω2) = (4R, 4Reiπ/4) である．
4.2
楕円柱列を過ぎる流れ
第2の例は，楕円柱列を過ぎるポテンシャル流である．すなわち，(2.1), (2.2)において D00= { x + iy������ x2 A2 + y2 B2 < 1 } . (4.3) ととった場合である．ここで，楕円の長軸の長さ2Aおよび短軸の長さ2Bは，ρ >R > 0 なる定数R, ρをとり， A = 1 2 ( ρ + R 2 ρ ) , B = 1 2 ( ρ− R 2 ρ ) (4.4) により与えている．この問題に対し，いくつかの基本周期の組(ω1, ω2)について本論文の 代用電荷法を適用してポテンシャル流を求めた．ただし，ρ =2Rととり，代用電荷法の拘 束点zi および電荷点ζj は次のようにとった． zj = J ( ρexp ( i2π( j − 1) N )) , ζj = J ( qρ exp ( i2π( j − 1) N )) , j = 1, . . . , N;R ρ < q < 1, (4.5) ここで，J(w)はJoukowski変換 J(w) = 1 2 ( w + R2 w ) (4.6)である．vj を与える二重積分は，第1の例と同様モンテカルロ法により計算した．Fig.4 にいくつかの基本周期の組(ω1, ω2)に対して計算したポテンシャル流の流線を示す．た だし，分点数は N = 32，電荷点ζj を(4.5)で配置するときのパラメータqはq = 0.6と とった． 0 2 4 6 8 0 2 4 6 8 Im(z) Re(z) 0 2 4 6 8 0 2 4 6 8 Im(z) Re(z)
(ω1, ω2) = (4R, 4Ri) (ω1, ω2) = (4R, 4Reiπ/3)
0 2 4 6 8 0 2 4 6 8 Im(z) Re(z) 0 2 4 6 8 0 2 4 6 8 Im(z) Re(z)
(ω1, ω2) = (4Reiπ/6,4Ri) (ω1, ω2) = (4R, 4Reiπ/4)
Fig. 4. The streamlines of the potential flows past a doubly-periodic array of elliptic obstacles.
この問題に対しても，(4.2)で与えられる境界条件に対する誤差評価ϵN を，電荷点ζj を 与える式(4.5)におけるパラメータqの値を変えながら，見積もった．ただし，qは制約 R/ρ = 0.5 < q < 1を受けるので，qはその範囲の値q = 0.6, 0.7, 0.8をとった．この結果 をFig.5に示す．図より，誤差評価ϵN は分点数N に対し指数関数的減衰することがわか る．二次元楕円内部領域における通常のポテンシャル問題に対し，拘束点・電荷点として 同心円上の等間隔点をJoukowski変換により写したものをとった従来の代用電荷法を適用 すると，誤差が指数関数的減衰することが理論的に示されているが[15]，周期的楕円列外 部ポテンシャル問題に対する本論文の代用電荷法でも同様の結果が成立することになる． 誤差の減衰率を最小二乗法により調べ，Table 2に示した．なお，ここで求めた誤差の減 衰はFig.5中のグラフにも破線で示した．表より，基本周期ω1, ω2 のとり方を変えても誤
である．vj を与える二重積分は，第1の例と同様モンテカルロ法により計算した．Fig.4 にいくつかの基本周期の組(ω1, ω2)に対して計算したポテンシャル流の流線を示す．た だし，分点数は N = 32，電荷点ζj を(4.5) で配置するときのパラメータqはq = 0.6と とった． 0 2 4 6 8 0 2 4 6 8 Im(z) Re(z) 0 2 4 6 8 0 2 4 6 8 Im(z) Re(z)
(ω1, ω2) = (4R, 4Ri) (ω1, ω2) = (4R, 4Reiπ/3)
0 2 4 6 8 0 2 4 6 8 Im(z) Re(z) 0 2 4 6 8 0 2 4 6 8 Im(z) Re(z)
(ω1, ω2) = (4Reiπ/6,4Ri) (ω1, ω2) = (4R, 4Reiπ/4)
Fig. 4. The streamlines of the potential flows past a doubly-periodic array of elliptic obstacles.
この問題に対しても，(4.2)で与えられる境界条件に対する誤差評価ϵN を，電荷点ζj を 与える式(4.5)におけるパラメータqの値を変えながら，見積もった．ただし，qは制約 R/ρ = 0.5 < q < 1を受けるので，qはその範囲の値q = 0.6, 0.7, 0.8をとった．この結果 をFig.5に示す．図より，誤差評価ϵN は分点数N に対し指数関数的減衰することがわか る．二次元楕円内部領域における通常のポテンシャル問題に対し，拘束点・電荷点として 同心円上の等間隔点をJoukowski変換により写したものをとった従来の代用電荷法を適用 すると，誤差が指数関数的減衰することが理論的に示されているが[15]，周期的楕円列外 部ポテンシャル問題に対する本論文の代用電荷法でも同様の結果が成立することになる． 誤差の減衰率を最小二乗法により調べ，Table 2に示した．なお，ここで求めた誤差の減 衰はFig.5中のグラフにも破線で示した．表より，基本周期ω1, ω2のとり方を変えても誤 差評価ϵN はあまり変わらず，大まかに ϵN =O(qN) に従うことがわかる． 1x10-16 1x10-14 1x10-12 1x10-10 1x10-8 1x10-6 0.0001 0.01 1 0 10 20 30 40 50 60 70 error N q=0.6 q=0.7 q=0.8 1x10-16 1x10-14 1x10-12 1x10-10 1x10-8 1x10-6 0.0001 0.01 1 0 10 20 30 40 50 60 70 error N q=0.6 q=0.7 q=0.8
(ω1, ω2) = (4R, 4Ri) (ω1, ω2) = (4R, 4Reiπ/3)
1×10−16 1×10−14 1×10−12 1×10−10 1×10−8 1×10−6 0.0001 0.01 1 0 10 20 30 40 50 60 70 error N q=0.6 q=0.7 q=0.8 1×10−16 1×10−14 1×10−12 1×10−10 1×10−8 1×10−6 0.0001 0.01 1 0 10 20 30 40 50 60 70 error N q=0.6 q=0.7 q=0.8
(ω1, ω2) = (4Reiπ/6,4Ri) (ω1, ω2) = (4R, 4Reiπ/4)
Fig. 5. The error estimates ϵN of the presented method applied to the examples shown in
Fig.4. In the figures, “error” means ϵN, and “q” is the parameter q used for positioning
the charge points ζjby (4.5).
なお，Fig.4の数値例について，モンテカルロ法積分により計算した平均流速は， ⟨vN⟩ U = (0.99652, 1.2 × 10−4) (ω1, ω2) = (4R, 4Ri) (1.00053, 2.5 × 10−4) (ω1, ω2) = (4R, 4Reiπ/3)
(1.00074, −3.3 × 10−4) (ω1, ω2) = (4Reiπ/6,4Ri)
(1.00003, −1.3 × 10−4) (ω1, ω2) = (4R, 4Reiπ/4)
Table 2. The decay rates of the errors ϵN of the presented method applied to the examples
shown in Fig.4, where q is the parameter used for positioning the charge points ζjby (4.5).
(ω1, ω2) q
0.6 0.7 0.8
(4R, 4Ri) O(0.58N) O(0.67N) O(0.77N) (4R, 4Reiπ/3) O(0.58N) O(0.67N) O(0.77N) (4Reiπ/6,4Ri) O(0.58N) O(0.67N) O(0.77N) (4R, 4Reiπ/4) O(0.57N) O(0.67N) O(0.77N)
4.3 Cassini
橙形列を過ぎる流れ
第3の例は，Cassini橙形列を過ぎる流れである．すなわち， ∂D00= x + iy ∈ C �� �� �� � ( xR0 +1 )2 +y2 ( xR0 − 1 )2 +y2 = a4 （R0 >0, a (1 < a ≤ √2)は定数）を境界に持つ障害物列の流れである．(1 < a < √2)の場 合，障害物Dmn は凸でなくなる．境界線は極座標(r, θ)を用いると r = r(θ) = R0 [ cos 2θ +(cos22θ + a4− 1)1/2]1/2, 0 ≤ θ ≤ 2π (4.7) と表される．この問題に対し，いくつかの基本周期の組(ω1, ω2)について本論文の代用電 荷法を適用してポテンシャル流を求めた．ただし，R0 =0.7R（R > 0は定数），a = 21/4と し，代用電荷法の拘束点zi および電荷点ζjは，[1]におけるとり方を若干変えて，次のよ うにとった． zi = r(θi)eiθi, θi = 2π(i − 1)N , i = 1, . . . , N, (4.8) ζj =zj+iq(zj+1− zj−1), j = 1, . . . , N, (4.9) ただし，z0 = zN, zN+1 = z1 である．この例でも，vj を与える二重積分はモンテカルロ法 積分により計算した．Fig.6にいくつかの基本周期の組(ω1, ω2)に対して計算したポテン シャル流の流線を示す．ただし，分点数はN = 64，電荷点ζj を(4.9)で配置するときの パラメータqはq = 1ととった． この問題における境界条件の誤差ϵN をFig.7に示す．図中，“q”は電荷点ζj を(4.9)で 与えるときのパラメータqの値である．図より，誤差ϵN はN を大きくするにつれて指数 関数的減衰することがわかる．代用電荷法の精度について，通常の二次元内部単連結領域 のポテンシャル問題の場合，境界が解析的な曲線であり，拘束点・電荷点として同心円上Table 2. The decay rates of the errors ϵN of the presented method applied to the examples
shown in Fig.4, where q is the parameter used for positioning the charge points ζjby (4.5).
(ω1, ω2) q
0.6 0.7 0.8
(4R, 4Ri) O(0.58N) O(0.67N) O(0.77N) (4R, 4Reiπ/3) O(0.58N) O(0.67N) O(0.77N) (4Reiπ/6,4Ri) O(0.58N) O(0.67N) O(0.77N) (4R, 4Reiπ/4) O(0.57N) O(0.67N) O(0.77N)
4.3 Cassini
橙形列を過ぎる流れ
第3の例は，Cassini橙形列を過ぎる流れである．すなわち， ∂D00 = x + iy ∈ C �� �� �� � ( xR0 +1 )2 +y2 ( xR0 − 1 )2 +y2 = a4 （R0 >0, a (1 < a ≤ √2)は定数）を境界に持つ障害物列の流れである．(1 < a < √2)の場 合，障害物Dmn は凸でなくなる．境界線は極座標(r, θ)を用いると r = r(θ) = R0 [ cos 2θ +(cos22θ + a4− 1)1/2]1/2, 0 ≤ θ ≤ 2π (4.7) と表される．この問題に対し，いくつかの基本周期の組(ω1, ω2)について本論文の代用電 荷法を適用してポテンシャル流を求めた．ただし，R0 = 0.7R（R > 0は定数），a = 21/4と し，代用電荷法の拘束点ziおよび電荷点ζjは，[1]におけるとり方を若干変えて，次のよ うにとった． zi =r(θi)eiθi, θi = 2π(i − 1)N , i = 1, . . . , N, (4.8) ζj =zj +iq(zj+1− zj−1), j = 1, . . . , N, (4.9) ただし，z0 = zN, zN+1 = z1 である．この例でも，vj を与える二重積分はモンテカルロ法 積分により計算した．Fig.6にいくつかの基本周期の組(ω1, ω2)に対して計算したポテン シャル流の流線を示す．ただし，分点数は N = 64，電荷点ζj を(4.9)で配置するときの パラメータqはq = 1ととった． この問題における境界条件の誤差ϵN をFig.7に示す．図中，“q”は電荷点ζj を(4.9)で 与えるときのパラメータqの値である．図より，誤差ϵN はN を大きくするにつれて指数 関数的減衰することがわかる．代用電荷法の精度について，通常の二次元内部単連結領域 のポテンシャル問題の場合，境界が解析的な曲線であり，拘束点・電荷点として同心円上 0 2 4 6 8 0 2 4 6 8 Im(z) Re(z) 0 2 4 6 8 0 2 4 6 8 Im(z) Re(z)(ω1, ω2) = (4R, 4Ri) (ω1, ω2) = (4R, 4Reiπ/3)
0 2 4 6 8 0 2 4 6 8 Im(z) Re(z) 0 2 4 6 8 0 2 4 6 8 Im(z) Re(z)
(ω1, ω2) = (4Reiπ/6,4Ri) (ω1, ω2) = (4R, 4Reiπ/4)
Fig. 6. The streamlines of the potential flows past a doubly-periodic array of Cassini’s ovals. の等間隔点を等角写像で写した点を用いた場合，誤差は指数関数的減衰することが理論 的に示されている[9, 20]．しかし，本数値例の場合，拘束点・電荷点はそのようなとり方 をしていない．それでも誤差が指数関数的減衰するのは意外な結果であった．誤差ϵN の 減衰率を最小二乗法により求めたのをTable 3に示す．なお，ここで求めた誤差の減衰は Fig.7 中のグラフにも破線で示した．基本周期ω1, ω2 のとり方を変えても誤差ϵN はあま り変わらないことがわかる．ただし，この例の場合，誤差ϵN の減衰率と電荷点ζj 配置を 決めるパラメータqとの間の定量的な関係はわからない． なお，Fig.6の数値例について，モンテカルロ法積分により計算した平均流速は， ⟨vN⟩ U = (0.99982, 2.9 × 10−4) (ω1, ω2) = (4R, 4Ri) (1.00021, 1.8 × 10−5) (ω1, ω2) = (4R, 4Reiπ/3)
(0.99957, −1.7 × 10−4) (ω1, ω2) = (4Reiπ/6,4Ri)
(0.99983, −3.3 × 10−4) (ω1, ω2) = (4R, 4Reiπ/4)
1×10−10 1×10−8 1×10−6 0.0001 0.01 1 0 10 20 30 40 50 60 70 error N q=1.5 q=1.0 q=0.5 1×10−10 1×10−8 1×10−6 0.0001 0.01 1 0 10 20 30 40 50 60 70 error N q=1.5 q=1.0 q=0.5
(ω1, ω2) = (4R, 4Ri) (ω1, ω2) = (4R, 4Reiπ/3)
1×10−10 1×10−8 1×10−6 0.0001 0.01 1 0 10 20 30 40 50 60 70 error N q=1.5 q=1.0 q=0.5 1×10−10 1×10−8 1×10−6 0.0001 0.01 1 0 10 20 30 40 50 60 70 error N q=1.5 q=1.0 q=0.5
(ω1, ω2) = (4Reiπ/6,4Ri) (ω1, ω2) = (4R, 4Reiπ/4)
Fig. 7. The error estimates ϵN of the presented method applied to the examples shown in
Fig.6. In the figures, “error” means ϵN, and “q” is the parameter q used for positioning
the charge points ζjby (4.9).
Table 3. The decay rates of the errors ϵN of the presented method applied to the examples
shown in Fig.6, where q is the parameter used for positioning the charge points ζjby (4.9).
(ω1, ω2) q
1.5 1.0 0.5
(4R, 4Ri) O(0.73N) O(0.78N) O(0.87N) (4R, 4Reiπ/3) O(0.72N) O(0.78N) O(0.87N) (4Reiπ/6,4Ri) O(0.72N) O(0.78N) O(0.87N) (4R, 4Reiπ/4) O(0.72N) O(0.78N) O(0.87N)
1×10−10 1×10−8 1×10−6 0.0001 0.01 1 0 10 20 30 40 50 60 70 error N q=1.5 q=1.0 q=0.5 1×10−10 1×10−8 1×10−6 0.0001 0.01 1 0 10 20 30 40 50 60 70 error N q=1.5 q=1.0 q=0.5
(ω1, ω2) = (4R, 4Ri) (ω1, ω2) = (4R, 4Reiπ/3)
1×10−10 1×10−8 1×10−6 0.0001 0.01 1 0 10 20 30 40 50 60 70 error N q=1.5 q=1.0 q=0.5 1×10−10 1×10−8 1×10−6 0.0001 0.01 1 0 10 20 30 40 50 60 70 error N q=1.5 q=1.0 q=0.5
(ω1, ω2) = (4Reiπ/6,4Ri) (ω1, ω2) = (4R, 4Reiπ/4)
Fig. 7. The error estimates ϵN of the presented method applied to the examples shown in
Fig.6. In the figures, “error” means ϵN, and “q” is the parameter q used for positioning
the charge points ζjby (4.9).
Table 3. The decay rates of the errors ϵN of the presented method applied to the examples
shown in Fig.6, where q is the parameter used for positioning the charge points ζjby (4.9).
(ω1, ω2) q
1.5 1.0 0.5
(4R, 4Ri) O(0.73N) O(0.78N) O(0.87N) (4R, 4Reiπ/3) O(0.72N) O(0.78N) O(0.87N) (4Reiπ/6,4Ri) O(0.72N) O(0.78N) O(0.87N) (4R, 4Reiπ/4) O(0.72N) O(0.78N) O(0.87N)
5.
まとめと今後の課題
本論文では二重周期的ポテンシャル問題に対する代用電荷法を，とくにポテンシャル流 問題に対して提案した．ポテンシャル流は流域における複素速度ポテンシャルという複素 解析関数で記述されるので，この問題は二重周期的領域における解析関数を然るべき境界 条件のもとで近似する問題に帰着される．この問題を代用電荷法で解く場合，近似すべき 解析関数が二重周期関数を含むので，従来の代用電荷法で解くのは難しい．そこで，周期 ポテンシャルをWeierstrass楕円関数を用いて構成し，それらの一次結合で解を近似する 方法を提案した．数値例により，本論文の有効性が示され，指数関数的収束という高い精 度を達成することが示された． 今後の課題として次の問題が挙げられる． 本論文ではポテンシャル流問題について議論し代用電荷法を提案したが，それ以外のポ テンシャル問題，例えば，二重周期性を持つ電磁場の問題についても，同様の代用電荷法 を構成することは簡単にできる． 本論文では二重周期的ポテンシャル問題に対し，Weierstrass楕円関数を用いて二重周 期関数を含むポテンシャル解を近似する代用電荷法を提案したが，テータ関数を用いた代 用電荷法を構成することもできる．そしてこれは，Weierstrassのシグマ関数のテータ関 数による表示(3.10)を用いれば，本論文で提案した方法と数学的には同値であることが示 される．これについては別稿で報告したい． 本論文の代用電荷法について，数値例より指数関数的収束という高い精度を与えること が示されたが，これに対する理論的根拠がほしい．本問題に限らず，代用電荷法に対する 理論解析は，二次元円板領域，楕円領域などごく限られた特別な場合についてしか行われ ていないのが現状であり，研究が急がれる． また，本論文で提案した方法を二重周期的ポテンシャル流から二重周期的 Stokes流 問題へ拡張することも，課題に挙げられる．これについては，著者らが [19] で周期的 Stokesletの重ね合わせによる代用電荷法をすでに提案しているが，ここで用いた周期的StokesletはFourier級数により与えたものである．一方，Weierstrass楕円関数を用いた周 期的Stokesletが[6]により与えられている．これを用いて，二重周期的Stokes流問題に
対し楕円関数に基づく代用電荷法を構成し，[19]の方法と結果に差異が見られるか，興味
が持たれる．
謝辞
本論文について貴重なご助言を頂いた査読者に感謝する．また，数値実験にあたり参考文献
[1] 天野要，代用電荷法に基づく外部等角写像の数値計算法，情報処理学会論文誌，
29(1988)，62–72．
[2] Amano, K., A charge simulation method for the numerical conformal mapping of inte-rior, exterior and doubly-connected domains, J. Comput. Appl. Math., 53 (1994), 353– 370.
[3] Amano, K., A charge simulation method for numerical conformal mapping onto circular and radial slit domains, SIAM J. Sci. Comput.,19 (1998), 1169–1187.
[4] Fairweather, G. and Karageorghis, A., The method of fundamental solutions for elliptic boundary value problems, Adv. Comp. Math.,9 (1998), 69–95.
[5] Greengard, L. and Kropinski, M. C., Integral equation methods for Stokes flow in doubly-periodic domains, J. Eng. Math.,48 (2004), 157–170.
[6] Hasimoto, H., Periodic fundamental solution of the two-dimensional Stokes equations, J. Phys. Soc. Japan,78 (2009), 074401.
[7] 今井功，流体力学（前編），裳華房，東京，1973．
[8] Liron, N., Fluid transport by cilia between parallel plates, J. Fluid. Mech., 86 (1978), 705–726.
[9] Katsurada, M., Asymptotic error analysis of the charge simulation method in a Jordan region with an analytic boundary, J. Fac. Sci. Univ. Tokyo, Sect. IA Math., 37 (1990), 635–657.
[10] Katsurada, M. and Okamoto, H., A mathematical study of the charge simulation method I, J. Fac. Sci. Univ. Tokyo, Sect. IA Math.,35 (1988), 507–518.
[11] Matsumoto, M. and Nishimura, T., Mersenne twister: a 623-dimensionally equidis-tributed uniform pseudo-random number generator, ACM Trans. Model. Comput. Simul.,8 (1998), 3–30.
[12] Milne-Thomson, L. M., Theoretical Hydrodynamics, Dover, New York, 2011.
[13] 村島定行，代用電荷法とその応用–境界値問題の半解析的数値解法，森北出版，東京， 1983． [14] 室田一雄，代用電荷法におけるスキームの「不変性」について，情報処理学会論文 誌，34 (1993)，533–535． [15] 西田詩，2次元楕円領域における代用電荷法の数学的及び数値的考察，日本応用数理 学会論文誌，5 (1995)，185–198．
参考文献
[1] 天野要，代用電荷法に基づく外部等角写像の数値計算法，情報処理学会論文誌，
29(1988)，62–72．
[2] Amano, K., A charge simulation method for the numerical conformal mapping of inte-rior, exterior and doubly-connected domains, J. Comput. Appl. Math.,53 (1994), 353– 370.
[3] Amano, K., A charge simulation method for numerical conformal mapping onto circular and radial slit domains, SIAM J. Sci. Comput.,19 (1998), 1169–1187.
[4] Fairweather, G. and Karageorghis, A., The method of fundamental solutions for elliptic boundary value problems, Adv. Comp. Math.,9 (1998), 69–95.
[5] Greengard, L. and Kropinski, M. C., Integral equation methods for Stokes flow in doubly-periodic domains, J. Eng. Math.,48 (2004), 157–170.
[6] Hasimoto, H., Periodic fundamental solution of the two-dimensional Stokes equations, J. Phys. Soc. Japan,78 (2009), 074401.
[7] 今井功，流体力学（前編），裳華房，東京，1973．
[8] Liron, N., Fluid transport by cilia between parallel plates, J. Fluid. Mech., 86 (1978), 705–726.
[9] Katsurada, M., Asymptotic error analysis of the charge simulation method in a Jordan region with an analytic boundary, J. Fac. Sci. Univ. Tokyo, Sect. IA Math.,37 (1990), 635–657.
[10] Katsurada, M. and Okamoto, H., A mathematical study of the charge simulation method I, J. Fac. Sci. Univ. Tokyo, Sect. IA Math.,35 (1988), 507–518.
[11] Matsumoto, M. and Nishimura, T., Mersenne twister: a 623-dimensionally equidis-tributed uniform pseudo-random number generator, ACM Trans. Model. Comput. Simul.,8 (1998), 3–30.
[12] Milne-Thomson, L. M., Theoretical Hydrodynamics, Dover, New York, 2011.
[13] 村島定行，代用電荷法とその応用–境界値問題の半解析的数値解法，森北出版，東京， 1983． [14] 室田一雄，代用電荷法におけるスキームの「不変性」について，情報処理学会論文 誌，34 (1993)，533–535． [15] 西田詩，2次元楕円領域における代用電荷法の数学的及び数値的考察，日本応用数理 学会論文誌，5 (1995)，185–198．
[16] Ogata, H., A fundamental solution method for three-dimensional Stokes flow problems with obstacles in a planar periodic array, J. Comput. Appl. Math.,189 (2006), 622–634. [17] Ogata, H. and Amano, K., A fundamental solution method for three-dimensional vis-cous flow problems with obstacles in a periodic array, J. Comput. Appl. Math., 193 (2006), 302–318.
[18] Ogata, H. and Amano, K., Fundamental solution method for two-dimensional Stokes flow problems with one-dimensional periodicity, Japan J. Indust. Appl. Math., 27 (2010), 191–215.
[19] Ogata, H., Amano, K., Sugihara, M. and Okano, D., A fundamental solution method for viscous flow problems with obstacles in a periodic array, J. Comput. Appl. Math., 152 (2003), 411–425.
[20] Ogata, H. and Katsurada, M., Convergence of the invariant scheme of the method of fundamental solutions for two-dimensional potental problems in a Jordan region, Japan J. Indust. Appl. Math.,31 (2014), 231–262.
[21] Ogata, H., Okano, D. and Amano, K., Numerical conformal mapping of periodic struc-ture domains, Japan J. Indust. Appl. Math.,19 (2002), 257–275.
[22] Singer, H., Steinbigler, H. and Weiss, P., A charge simulation method for the calculation of high voltage fields, IEEE Trans. Power Appar. Syst.,PAS-93 (1974), 1660–1668.
[23] 竹内端三，楕圓凾數論，岩波書店，東京，1936．
[24] 梅村浩，楕円関数論 増補新装版：楕円曲線の解析学，東京大学出版会，東京，2020．
[25] Zick, A. A. and Homsy, G. M., Stokes flow through periodic array of spheres, J. Fluid. Mech.,115 (1982), 13–26. 緒方 秀教(正会員) 〒182-8585 東京都調布市調布ケ丘1-5-1 1992年東京大学大学院工学系研究科修士課程修了，博士（工学）．現在，電気通信大学 大学院情報理工学研究科情報・ネットワーク工学専攻教授．日本応用数理学会，日本数学 会，情報処理学会，日本計算工学会会員．数値解析，とくに関数論・佐藤超函数論に基づ く数値計算手法に興味を持つ． （受付日 2020 年 8 月 9 日） （受理日 2021 年 1 月 26 日）
