数理解析研究所講究録
1539 巻 2007 年 1-14

1

時間領域音響計算に用いる compact 差分と多段階積分法の最適化 (Optimization of Linear
Multi- step Methods for Acoustic Wave
Compu- tation)
岩津 玲磨 (Reima Iwatsu)”
鶴 秀生 (Hideo Tsuru)2

Denki University
2Nittobo Acoustic Engineering Corporation Ltd.
$*\mathrm{l}\mathrm{T}\mathrm{o}\mathrm{k}\mathrm{y}\mathrm{o}$

1

はじめに

音響のシミュレーションには時間領域または周波数領域において, 有限差分法, 有限要素法, 境界要素法などの解析手法が用いられる. そのうち線形な音波の伝播のシミュレーションに関しては, 従来から使われてきた FDTD 法 (1) が依然として広く用いられている. FDTD 法は, 2 次中心と Leap Rog 法
によって, 時間方向に staggered 配置さ
差分$(\mathrm{E}\mathrm{D}2)$$(\mathrm{L}\mathrm{F})$れた速度と圧力を安定に時間積分していく方法となっている. 簡単に中庸な精度の解が得られるものの, 高精度な解は得られない.
その–方で計算空力音響 (CAA) の分野では最近, Euler 方程式などをもとに,に代表される高解像度のスキームや, 高精度 Compact 差
分スキーム (3) と, 圧縮性流体むけに開発されてきた総変動減少スキームを組み合わせて, 4 次 Rllnge-Klltta 法
で時間積分するような計算法 (4) がしばしば採用されている.
そこで FDTD をもとに高精度な計算をする目的のもとに, 線形多段階積分法 (LMM) を適用して, 音響計算用の時間積分法の最適化をおこなった. 最適化の方法は Tamn] と Web による手法 (2) を用いた. ここでは精度の保持に必要以上な段階数を導入することによって, そこから出てくる余分の係数の値を変数として, 時間積分法のもつ実効角振動数に含まれる誤差を, 対象とする振動数領域内で最小とするように, その値を定めた. 時
間積分に LMM を選択する理由は, 時間方向に staggered 配置されている変数をそのまま Runge-Kutta 法で取り扱えないことによる. ここで扱う
種類の LMM は, 調べた範囲では過去の文献に取り扱われていない. ODEによく用いられる Adams-Molllton 法=Adams-Bushforth 法の予測子=修$\mathrm{D}\mathrm{R}\mathrm{P}^{(2)}$$(\mathrm{R}\mathrm{K}4)$

2
正子対など (5) と異なり, 時間方向に半メッシュずれた変数配置になっているため, 過去の文献に取り上げられる機会がなかったものと思われる.
計算誤差は時間, 空間両方の離散化から入ってくるので, 空間方向にも近似度を上げ, 不等間隔 4 次精度 compact 差分スキーム (CDuns4) を用いた. ただし本研究では, 複雑形状の実用計算に容易に応用可能なように, ステンシルの広がりが最小であるようなスキームに限った. 直交不等間隔格子を用いて複雑形状を表わす目的で, 格子間隔比とこのスキームの波数に含まれる誤差の関係を調べた. 差分スキームの実効波数は格子間隔が不均–なときに虚部成分をもつため, 代表的な格子生成関数に対して格子間隔比と位相誤差の関係を調べて, これを目標値以下におさえるために必要な条件を求めた.
ベンチマークテストの結果では最適化された LMM は, 従来法よりも 1 桁以上小さい誤差を示した. さらに, 最適化された LMM が DRP や$\mathrm{R}\mathrm{K}4$

に対してどの程度の競争力をもつのかを調べた.

2

計算方法

2.1

不等間隔格子に対する staggered comPact スキーム

1) 内点に対するスキーム
内点に対する staggered compact スキームとして最小ステンシル幅をもつ 4 次精度スキーム (CDllns4) のみを扱うことにする.

.

(1)$\alpha f_{1-1}’.+f’+\beta f_{\mathfrak{i}+1}’=a(f_{i+\frac{1}{2}}-f_{i-\frac{1}{2}})+\epsilon$

上式の係数は, 格子間隔を用いて表わされる. 精度は, 滑らかな格子 $h_{2}-$$h_{1}=O(h^{2})$ に対して 4 次,
, 滑らかでない格子に対して 3 次,を示す. また, 等間隔格子の場合 Pad\’e タイプの公式と–致する.
2) 境界点に対するスキーム
FDTD の計算の仕組みを利用すると, 境界点に対して 2 種類のスキームがと名づける. 境界上で必要となる. これらを便宜的に境界スキーム
速度が定義されているものとすると, 圧力, 速度に対してそれぞれ境界ス
はそれぞれ,
キームが適用される. 3 次精度境界スキーム$\epsilonarrow O(h^{4})$$\epsilonarrow O(h^{3})$$\mathrm{A},$$\mathrm{A},$$\mathrm{B}$$\mathrm{B}$$\mathrm{A},$$\mathrm{B}$$f_{i}’+\beta f_{i+1}’=a(f_{i+\frac{1}{2}}-f_{i-\frac{1}{2}})+b(f_{i+_{\overline{2}}};-f_{i-\frac{1}{2}})+\epsilon$

,

(2)

,

$(_{\iota}’\})$$f_{i}’+\overline{\beta}f_{i+1}’=\overline{a}(f_{i+\frac{}{2}}., -f_{i+2})+\overline{b}(f_{i+\frac{6}{2}}-f_{i+_{2}^{1}})+\overline{\epsilon}$はのように表わされる (図 1 参照). 等間隔格子に対しては, スキーム
それぞれ $\beta=-1,$ $a=2,$ $b=-1,\overline{\beta}=23,\overline{a}=26,\overline{b}=-1$ となる. 誤差の$\mathrm{A},$$\mathrm{B}$3
主要項の係数はまったく同じ (1/24) だが, 左辺非対角項の係数の大きさが異なる. スキーム の方が係数行列を解くときの条件が悪くなる. 以下の数値実験において, それは実際に誤差を大きくすることに反映される.
$\mathrm{B}$$i^{\frac{ihi+1}{-\frac{1}{2}i+\frac{1}{2}ig_{1}g_{2}g_{3}}}+ \frac{3}{2}ihi+1i+2g_{1}g_{2}g_{2}i+\frac{1}{2}i+\frac{\ovalbox{\tt\small REJECT}_{3}}{2}i+\frac{6}{2}$Case A
Case
Fig. 1: Boundary at grid point$\mathrm{B}$$i$3) 格子数と娯差いずれかを
内点 4 次精度スキーム CDuns4 に境界 3 次精度スキーム
組みあわせたスキームの誤差を, $y=x^{8},$ $y=\sin\pi x,$ $0\leq x\leq 1$ を試験関数として, 等間隔, 不等間隔 (hyperbolic) 格子の格子数を変えて求めた.
スキーム
どちらも平均で 4 次精度を示すが, 誤差の大きさはスキームの方が大きい また sine 波では 1 波長に 8 格子点あれば, 1 階導関数の程度にすることができる.
誤差を$\mathrm{A},\mathrm{B}$$\mathrm{A},\mathrm{B}$$\mathrm{B}$$10^{-3}$

Fig. 2: Error vs , Left :$h$(even grid, hyperbolic grid), right :
(
)
grid)
grid,

$\mathrm{C}\mathrm{D}1\ln_{\iota}\backslash ^{\backslash }4\cdot 3\mathrm{A}$$\mathrm{C}\mathrm{D}\mathrm{l}\mathrm{l}\mathrm{n}\mathrm{k}\mathrm{s}\sim 4’..3\dagger 3$($\mathrm{e}\mathrm{r}\mathrm{t})$$‘..\mathrm{v}\mathrm{e}^{1},\iota \mathrm{l}$$.7\mathrm{l}\mathrm{i}\mathrm{c}$$\mathrm{h}\mathrm{y}\mathrm{I}$

つぎに同じ順番で空間誤差分布を図 3 に示す. 理論が示すとおり, 境界ス
キームの誤差が内部に深く浸透しないことがわかる.

Fig. 3: Errorof(right), Even grid$\mathrm{e}\mathrm{l}\mathrm{i}\mathrm{s}|_{J}1.\cdot \mathrm{i}\mathrm{b}\mathrm{t}\mathrm{l}\mathrm{t}_{t}\mathrm{i}\mathrm{o}\mathrm{n}$$\mathrm{C}\mathrm{I}\mathrm{J}\mathrm{t}1\mathit{1}1\mathrm{b}^{\mathrm{t}}4+$$\mathrm{a}\mathrm{J}\mathrm{J}\mathrm{d}$scheme A (left) andbyperbolic grid})($|\ln(\mathrm{l}\mathrm{a}\mathrm{r}\mathrm{y}$$\mathrm{B}$44) 不等間隔格子上での CDuns4 の実効波数$-x_{i}=r(x_{i}-x_{i-1}),$
等比数列格子を例にとって格子間隔比がスキームの実効波数 w/ に及ぼす影響を調べた. それによるならば波数の範囲で誤差を $|Re(w’)-w|<10^{-2}$ ,と,$|I7n(w’)|<10^{-}2$ におさえることができる. 図 4 左の 2 図は等比数列格
子 $(r=1,1.02, \cdots, 1.1)$ , 右 2 図は等間隔格子 $(x\leq x_{i})$ を等比数列格子$(x..\geq x_{i})$ に接続させたつなぎ目の点における実効波数を示す.
$x_{i+l}$$7^{\cdot}\leq 1.02$

Fig. 4:$x_{i+1}-x_{i+\frac{1}{2}}=r(x_{i+\frac{1}{2}}-x_{i})$$\alpha\Delta x\leq 0.518$

of, righf,:$\mathrm{C}\mathrm{D}\mathrm{t}\mathrm{l}\mathrm{r}1_{\iota}\backslash .’4\mathrm{f}\mathrm{o}\iota.\cdot$$\mathrm{F}_{\lrcorner}\mathrm{f}\mathrm{f}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{i}.\mathrm{v}‘\supset,\mathrm{f}\mathrm{r}\epsilon^{1}.,\mathrm{q}\iota 1\mathrm{c}^{\mathrm{n}}.,\mathrm{r}\downarrow \mathrm{c}\mathrm{y}$no$\iota.n\cdot\iota \mathrm{i}\mathrm{f}\mathrm{o}\iota\cdot \mathrm{r}\mathrm{l}\mathrm{l}\mathrm{r}\mathrm{r}1\epsilon^{1},\mathrm{s}$$\mathrm{I}_{I}\mathrm{e}\mathrm{h}:\zeta_{J}^{1}\zeta].\mathrm{t}\mathrm{l}\mathrm{i}- Y\mathrm{e}\gamma_{\mathrm{I}}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{r}\mathrm{n}^{1}‘.,.\backslash ’ \mathrm{h}$$‘!,\mathrm{q}\iota.1\mathrm{i}-\Gamma’t\iota$

2.2

tio

)

$\mathit{1}\triangleright\iota.\iota \mathrm{n}\mathrm{i}\{\dot{(}$}$\mathrm{r}\mathfrak{r}\iota 1\mathrm{e}\mathrm{c}.$]$\mathrm{t}\mathrm{l}\mathrm{i}- Y_{\mathrm{f}}^{l}1_{t}\mathrm{t}\mathrm{i}\mathrm{o}\iota \mathrm{m}‘.d\backslash \backslash ’’.\mathrm{h}$,;ted to$\iota(.\prime c\leq.\prime r_{\mathrm{i}}.,)(^{\mathrm{u}}.,\mathrm{o}\mathrm{n}\mathrm{r}\iota\Leftrightarrow(^{\backslash }$$\mathrm{r}\mathrm{r}\mathrm{l}\mathrm{C}_{d\mathrm{L}}^{-\backslash \urcorner}‘,\mathrm{h}(x\geq.\prime L_{?}^{},)$

線形多段階積分法 (LMM) の最適化

最適化の手法には幾とおりかの方法が考案されているが, 基本的な考え方は, 精度の保持に必要以上な数の変数を導入し, 余分な自由度をパラメータとして変化させることによって, 対象とする波数の範囲に対して
数値積分スキームのもつ dissipation error と dispersion enor を最小にしつつ, 安定性の上限を最大化する点で–致する. ここでは, 速度と圧力の定義点が時間方向に半格子ずらされている配置の特徴を生かしたまま LF
を多段階に拡張し, Tam $Webb の方法 (2) で 2,3 次精度の最適化 LMM
を導いた.

2.2.1

移流方程式に対する線形多段階積分法 (LMM) の解析

1) 従来型 FDTD(Leap-Frog) 法 $(M_{1})$$\frac{f^{n+1}-f^{n}}{\Delta t}=F^{\mathrm{n}+_{2}^{1}}-\frac{\triangle t^{2}}{24}f_{ttt}$,$F=f_{t}$(4)5
簡単のために $f=0(t<0)$ と仮定して Laplace 変換をおこなうと,実効角振動数 は$\mathrm{L}\mathrm{F}$の$\overline{\omega}$$- \dot{l}\frac{i(e^{-i\omega\Delta t}-1)}{\Delta te^{-iv\Delta t/2}\{}\tilde{f}$$\frac{\tilde{d}f}{dt}\equiv-\dot{\iota}$$-\cdot$$- \tilde{f}arrow-.\Delta t_{J}=2si_{7l}(\frac{\omega\Delta t}{2})$$Re( \overline{\omega}\Delta t)=2sin(\frac{\omega\triangle t}{2})$,,$(_{\mathrm{t}}^{r_{)}})$(6)$I7n(\overline{\omega}\Delta t)=0$と求まり, 実数となることがわかる. このことが, 従来型 FDTD が広く用いられてきた理由と推測される.
2) 2 段階法 $(M_{2})$との数を 2 段階に拡張し, そのとき時間ステップ
つぎに
値流束の重み係数 を変数として変化させて, 時間積分法の性質を改善できないかどうか調べることにする. 時間積分法およびその実効角振動数は以下のようになる.
$n+ \frac{1}{2}$$M_{1}$$n- \frac{1}{2}$$b_{0}$(7)$\frac{f^{n+1}-f^{n}}{\Delta t}=b_{0}F^{n+\frac{1}{2}}+(1-b_{0})F^{n-\frac{1}{2}}+\Delta t(b_{0}-1)f_{tt}$$\overline{\omega}\Delta t=\frac{2sin\frac{\theta}{2}(cos\theta+b_{0}(1-cos\theta)+i(b_{0}-1)\mathrm{s}in\theta)}{1+2b_{0}(b_{0}-1)(1-cos\theta)}.$,(8)$\theta=\omega\Delta t$

ここでは Tam&Webb による最適化法を採用した. 対象とする角振動数において, 時間積分法の実効角振動数の実部, 虚部に含まれる 2を最小にするように の値を定める. 実部と乗誤差の加重平均の積分
虚部の重要さを表わす係数 の値は $0<\sigma<1$ である.の範囲$\eta$$b_{0}$$E_{1}$$\sigma$$E_{1}= \int_{-\eta}^{\eta}(\sigma Re(\overline{\omega}\Delta t_{J}-\omega\Delta t)^{2}+(1-\sigma)Im(\overline{\omega}\Delta t-\omega\Delta t)^{2})d(\omega\Delta t)$$\theta\sim 0$,$2sin \frac{\theta}{2}\sim 2(\frac{\theta}{2}-\cdots)$,$sin\theta\sim\theta-\cdots$,$Im(\overline{\theta})\sim(b_{0}-1)\theta^{2}$,$cos\theta\sim 1-\cdots$$\overline{\theta}=\overline{\omega}\Delta t$$Re(\overline{\theta})\sim\theta$

,

,

(9)

$E_{1}= \int_{-\epsilon}^{\epsilon}=\frac{2}{15}(1-\sigma)((1-b_{0})^{2}\epsilon^{2}+1)\epsilon^{3}$$\frac{dE_{1}}{db_{0}}=0arrow 2(b_{0}-1)\epsilon^{2}=0$$\theta=0$の近傍$|\theta|\leq\epsilon$$(M_{1})$$b_{0}=1$(1.0)

を極小化する. すなわち, 角振動が最良であることがわかる.

で $b_{0}=1$ が数が小さいときには従来法
3) 3 段階法,$E_{1}$$(M_{3})$$\frac{f^{n+1}-f^{n}}{\Delta t}=b_{0}F^{n+\frac{1}{2}}+(2-2b_{0})F^{n-_{2}^{1}}+(b_{0}-1)F^{n-\frac{\theta}{2}}+\frac{\Delta t^{2}}{24}(24b_{0}-25)f_{ttt}(11..)$

6
簡単のため $(f=0, t<\Delta t)$ と仮定して両辺を Laplace 変換し,角振動数を求めると以下のようになる.

$\mathrm{A}f_{3}$

の実効

(12)

$-i \frac{i(e^{-i\omega\Delta t}-1)}{\Delta t\sum_{j=0}^{2}b_{j}e^{-\dot{u}d(\frac{1}{2}-j)\Delta t}}\tilde{f}\approx\frac{\tilde{d}f}{dt}\equiv-i\overline{\omega}\tilde{f}arrow\overline{\theta}=\frac{2sin(\frac{\theta}{2})}{\Sigma_{j=0}^{2}b_{j}e^{ij\theta}}$

ただし$\omega\Delta t=\theta,\overline{\omega}\Delta t=\overline{\theta}$

で, 係数$\sigma$の値は文献 [2] に推薦されている値を用いた b。の最適値を 最小化によって決めたところ, つぎのような値, $b_{0}=1.00553341(\mathit{7}|=\pi/2, \sigma=0.36)$ が求まった. 積分誤差ののある値において実際にが最依存性を図 5(の左) に示す.

$(\sigma^{t}=0.36)$$E_{1}$$E_{1}$$b_{0}$$b_{0}$$E_{1}$

小となっていることが観察される. この値は対象とする角振動数の範囲にの実効角振動
よって, 表 1 のように $b_{0}=1$ のまわりで若干変化する.
数が $0.9\leq b_{0}\leq 1.1$ の範囲でどのように変化するのかを, 実部, 虚部にわけて図 6 に示す. 最適化によって, ごくわずかな位相誤差を許容するかわりに, 全体的な偏差が少なくなっていることがわかる.
$M_{3}$

Fig. 5:$\mathrm{I}^{\neg}.\{\mathrm{i}\mathrm{g}$.

$b_{0}(1\mathrm{e}^{1},\mathrm{p}(^{\backslash },\cdot\cdot \mathrm{n}\mathrm{t}\iota \mathrm{e}\mathrm{n}\mathrm{c}\prime \mathrm{y}$of integraJ deviation in$6:b_{0^{\mathrm{d}\mathrm{e}}1\prime \mathrm{n}\mathrm{c}\mathrm{y}}.\supset \mathrm{e}\mathfrak{n}(1‘.\backslash$of real part$\overline{\theta}$(left:$\dot{r}1\mathfrak{l}\mathrm{I}\mathrm{l}\mathrm{d}\mathrm{i}\mathrm{n}’\mathrm{l}\mathrm{a}|\mathrm{g}\mathrm{i}\mathfrak{n}\mathrm{a}\mathrm{r}\mathrm{y}$$M\prime 3\tau \mathrm{r}\mathrm{i}8\mathrm{h}\mathrm{l}J:\Lambda f_{4}$part of$\frac{\mathrm{T}\mathrm{a}\mathrm{b}\mathrm{l}\mathrm{e}1:\mathrm{I}\mathrm{n}\mathrm{f}\mathrm{l}\tau 1\mathrm{e}\mathrm{n}\mathrm{c}\mathrm{e}\mathrm{o}\mathrm{f}\eta \mathrm{o}\mathrm{m}b_{0}}{M_{3}\eta b_{0}M_{4}\eta b_{0}}$$\pi$$\pi/2$1.0
0.5

0.950878859
1.00553341
101911891
103340232$\pi$$\pi/2$

1.0
0.5

c.a.1.010380775

1.07340406
10843831$\overline{\theta}(\Lambda I.’;)$)

74) 4 段階法$(\Lambda f_{4})$ftttt,(13)$b_{3}= \frac{25}{24}-b_{0}$(14)$\frac{f^{n+1}-f^{n}}{\Delta t}=\sum_{j=0}^{3}b_{j}F^{n+\frac{1}{2}-j}+\frac{\Delta t^{4}}{12}(12b_{0}-13)$$b_{1}=-3b_{0}+ \frac{73}{24},$$b_{2}=3b_{0}- \frac{37}{12},$

3 段階積分法と類似の手順によって,最小化によって決の最適値をめると $b_{0}=$ 1.0380775 $(\eta=\pi/2, \sigma=0.36)$ のような値が得られた.
の実効角振動数の依存性は図 5(の右), 最適値の 依存性は表 1,$b_{0}$$E_{1}$$E_{1}$$b_{0}$$b_{0}$$\eta$$M_{4}$

は図 7 に示されている.

$\mathrm{F}^{\urcorner}\mathrm{i}\mathrm{g}$.

$7:b_{0}\mathrm{e}:1\mathrm{e}\mathrm{p}\mathrm{t}_{d}^{-\backslash }1$)($1.\mathrm{e}^{\mathrm{Y}}\mathrm{n}\mathrm{c}\mathrm{y}$of real$\mathrm{p}\mathrm{a}\mathrm{I}^{\cdot}\mathrm{t}_{1}\mathrm{a}\iota 1(1\mathrm{i}\mathrm{m}_{C}\backslash \mathrm{g}\mathrm{i}\mathrm{n}r11\}^{r}$part of$\overline{\theta}(\Lambda\cdot\cdot\ell_{4})$5) 多段階法 $(M_{3}, M_{4})$ の安定性限界
Tam と Webb による最適化手法には, 安定性に関する考慮が含まれていな
い. そのために, 時間積分法の増幅率 を計算して安定性の上限を求めた.
$r$$r= \frac{\tilde{f}^{n+1}}{\tilde{f}^{n}},$$r_{exad}=e^{-i\omega\Delta t},$$r_{M\circ}=1-i \omega\Delta t\sum_{j=0}^{\epsilon-1}b_{j}e^{-i(\#}-j)\omega\Delta t$$\frac{r_{M_{*}}}{r_{exac\ell}}=e^{1\omega\Delta t}-i\omega\Delta t\sum_{j=0}^{s}b_{j}e^{i(\frac{1}{2}+j)\omega\Delta t}=r_{(\tau)}e^{-i\delta}$

(15)
(16)

は増幅率比,\mbox{\boldmath $\delta$} は位相誤差を表わす (6). 安定性の上限はで定義される. 減衰, 分散誤差上限の基準はそれぞれ $|r_{(r)}-1|\leq 10^{-3}$ ,は $b_{0}<1.08$ であれば安定であると定めた. 図 8 によると,がにおいて不安定になるという結果を得た. これはあきら
が,かに経験的事実に反する. したがって, 安定性を議論するためには, 移流
ここで$|\delta|\leq 10^{-3}$$M_{3}$$r_{(\mathrm{r})}\leq 1$$r_{(r)}$$M_{4}$$b_{0}\simeq 1$方程式を 2 回組み合わせた波動方程式を対象として考えないと正しい結果に到達できない.

8Fig. 8: Accuracy limit of ampl.ification (red), stability limit (blue), right ::),accuracy limit of phase error (gi ) vs (denotes$r\mathfrak{U}1\mathrm{d}$$f_{J}0$$.\cdot\epsilon^{\mathrm{Y}},\mathrm{e}\mathrm{n}$$\mathrm{I}.\lrcorner \mathrm{e}\mathrm{f}\mathrm{f}_{1}$$M_{4}$$\lambda/\mathit{1}_{3}$$m$$\omega\Delta t$2.3

波動方程式に対する線形多段階積分法の安定性解析-Frog) 法
1) 従来型 FDTD (説明のために 1 次元の場合を考える. 圧力 に対する波動方程式に中間変数 を導入して, 以下のように連立方程式を解く. したがって, 時間空間
とも対称スキームで積分していることにな る.
$(M_{1}\cdot M_{1})$$\mathrm{L}\mathrm{e}\mathrm{a}\mathrm{p}$

$P$

$u$

$u_{t}=- \frac{1}{\rho_{0}}p_{x}$(17)$p_{tt}=c_{0}^{2}p_{xx}arrow$$p_{t}=-p_{0}c_{0}^{2}u_{x}$2) 線形多段階積分法 $(M_{q}^{2}\equiv M_{r+s-1}^{2}=M_{s}\cdot M_{f})$発展方程式ん $=H(f)$ を積分するとき, 中間変数 を介してゐ $=F(g)$ ,$g_{t}=G(f)$ のようにして解く. 各位の積分に多段階積分法
を用いる. その結果は, ん $=H$ を 2 階の方程式に対する多段階積分法
で積分したことと同じになる.
$g$$M_{s},$$M_{r}$$M_{\epsilon+r-1}^{2}$$f^{n+1}=f^{n}+ \Delta t\sum_{j=0}^{s-1}b_{j}F^{n+\frac{1}{2}-j}+\epsilon$

,

(19)

$g^{n+\frac{1}{2}}=g^{n-f}+ \Delta t\sum_{k=0}^{r-1}c_{k}G^{n-k}+\epsilon’1$,$(2())$.

(21)$f^{n+1}=2f^{n}-f^{n-1}+ \Delta t^{2}\sum_{l=0}^{q-1}d_{l}H^{n-l}+\epsilon’’$

ただし $H(f),$ $F(g),$ $G(f)$ は $H\equiv FG=GF$ をみたし, $q=S+r-1$ ,である.

$d_{l}=\Sigma_{j=0}^{l}b_{l-j^{C}j}$$M_{3},$$M_{4}$の組みあわせに対する$M_{q}^{2}$の段数と誤差の主要項を表 2 に示す.

9
Table 2: 2-parameter family of$M_{q}^{2}=M_{s}\cdot\Lambda f_{r}$.

$\frac{s\gamma\cdot q\epsilon’’}{335\triangle t^{2}(b_{0}+c_{0}-\frac{25}{12})f^{Iv}}$

4 4 7$\Delta t^{3}(b_{0}+c_{0}-\frac{13}{6})f^{V}$3) 波動方程式に対する線形多段階積分法 (LMM) の解析に対する実効角振動数はの組み合わせ以下のようになる.
$(M_{q}^{2}\equiv\Lambda I_{s+r-1}^{2}=M_{s}\cdot M_{r})$$M_{\mathit{8}},M_{r}$(22)$\frac{r_{M_{q-1}^{2}}}{r_{exact}}=2e^{\iota’\omega\Delta t}-e^{i2\omega\Delta t}-(\omega\triangle t)^{2}\sum_{l=0}^{q-1}d\iota e^{i(l+1)\omega\Delta t}=r_{(r)}e^{-i\delta}$$q=s+r-1,$$d_{l}= \sum_{j=0}^{\iota}b_{l-j^{C_{J}}j},$(23)$b_{j}=c_{j}$$|r_{(r)}-1|\leq 10^{-3},$安定性, 減衰, 分散誤差上限の基準をのように設定して, それらの 依存性を調べると図 9 のようになる. 安の範囲は,定領域で減衰, 分散誤差上限を大きくするようなに対してに対しておおよそ 1.005
1.043,$1.02\leq b_{0}\leq$ 1.676,の最となる. したがって, 上で導出した
適値は安定領域に入っていることが示された.
$|\delta|\leq 10^{-3}$$r_{(r\rangle}\leq 1,$$b_{0}$$M_{3}\cdot M_{3}$$b_{0}$$\omega\Delta t\leq 0.4,$$\leq b_{0}\leq$$M_{4}\cdot M_{4}$$\omega\Delta t\leq 0.4$$\mathrm{b}^{\backslash }\mathrm{i}\mathrm{g}.\backslash ‘):\Lambda_{\mathrm{C}\mathrm{C},11\mathrm{r}t\backslash (^{\backslash }}\cdot.,\mathrm{y}$$b_{0}$

limit of ”
linlit of)$\mathrm{m}\prime \mathrm{p}\mathrm{l}\mathrm{i}\mathrm{f}\grave{\iota}\mathrm{c}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{t}$$t$$\mathrm{p}1\mathrm{l}’\mathrm{e}1\mathrm{B}\mathrm{e}$$\mathrm{a}(^{\mathrm{Y}},.\mathrm{C}\mathrm{l}\mathrm{l}\mathrm{r}_{C}\backslash \mathrm{c}^{\backslash }\mathrm{y}$3

数値計算結果3.1

1 次元ベンチマーク問題(red), stability l.imit (blne) andcrro], (green) vs .

$\mathrm{n}$$b_{0}$

移流方程式の初期値問題$f_{t}+f_{x}=0$$f(0)= \frac{1}{2}\exp(-\ln(2)(\frac{x}{3})^{2}),$,$-20\leq x\leq 450,$ $\Delta x=1.0$(24)$(2_{\dot{\mathrm{t}}}^{\mathrm{r}_{)}})$

10の解 (4) を各種スキームを用いて計算して,において厳密解と図(2)10 に比較した. DRP は文献 の係数を用いた.
は Carpenter とKennedy による 5 段階 4 次精度最適化 Rllnge-Kutta 法 (7) で, 最小ステン
シルをもつ 4 次精度 regular compact スキームと組みあわせて使用した.
従来法 (FDTD) は Courant 数 1.0 で最良の結果を示し, 波形をほとんど変えず, 群速度がわずか 0.3% 小さい. DRP は Courant 数 0.25 以下で安定
は Courant 数 1.98 まで安定に積分できた. 従来法であった. 最適化の E2 を CDs4 に置換するだけで改善がみられた. 2 種類の最適化されたより良かった.
多段階積分法を比べると, 輪の方がで計算された解は Courant 数 0.2 以下で厳密解とほとんど–致する.
を Courant 数 0.1 において最小と数値解と厳密解の偏差の 2 乗和
$b_{\text{。}}=1.0101515$
した b。の値は$b_{0}=1.0405254$ であった.
最適化によって得られた値 M 勧 t, $b_{0}=1.00553341$ ,$b_{0}=1.0380775$
より若干小さめであった (図 11).
を Cotlrant 数に対してプロットしたのが図 12 である.の偏差がほぼ–定で広い安定領域が目立っている. 従来法 (FDTD) は Courant
数 1 をごくわずか超える約 1008 まで安定である. LMM の安定限界はCotlrant 数 0.83 程度である. しかし LMM は FDTD よりも 1-3 桁誤差がが小さい. 注目されることは, Courant 数 0.2 以下で, DRPよりも正確な点である.
$t_{\ovalbox{\tt\small REJECT}}=400$$\mathrm{R}\mathrm{K}4\mathrm{C}\mathrm{K}$$\mathrm{R}\mathrm{K}4$$M_{4}$$M_{3}$$dev_{2}$$M_{3}\cdot M_{3},$$M_{4}\cdot M_{4},$$M_{4\varphi t},$$dev_{2}$$\mathrm{R}\mathrm{K}4$$M_{3}\cdot M_{3}$3.2$\mathrm{R}\mathrm{K}4$球による 3 次元波面反射

音源 ( 点) から 1(m) の点に中心 ( 点) がある半径 $a=0.2(\mathrm{m})$ の球からの波面反射を従来法と LMM で計算して得られた数値解を厳密解と比較した (8). 腰元は以下のとおりである.
平均格子間隔: $x,$ $y,$ 方向 $\Delta x=0.02(\mathrm{m})$等間隔格子の場合: $\Delta x=0.02(\mathrm{m})$不等間隔格子の場合: $\Delta x_{m\iota’n}=0.0083$
格子点数:
時間刻み:
(msec) (固定)
音の振動数: 500, 1000, 2000厳密解は, 無限級数を加算する順序と数値誤差の入り方, 項数に十分注意$\mathrm{O}$$\mathrm{S}$$z$$(\mathrm{m})$$101^{3}$$\Delta t=0.\mathrm{O}1$$(\mathrm{H}\mathrm{z})$$(\mathrm{H}\mathrm{z})$$(\mathrm{H}\mathrm{z})$11

FDTD(LF+E2)$\mathrm{D}\mathrm{R}\mathrm{k}_{9:}^{)}.,$$f.\mathrm{t}\mathrm{K}4\mathrm{C}\mathrm{K}$$\Lambda\prime f_{1}\cdot M_{1}+\mathrm{C}’ \mathrm{D}\mathrm{s}43$$\mathrm{o}\mathrm{p}\mathrm{t}\mathrm{A}\mathit{1}_{4}\cdot M_{4}+\mathrm{C}\mathrm{D}\mathrm{s}43$$(\mathrm{P}\mathrm{t}\Lambda\parallel:s\cdot\Lambda f:;+\mathrm{C}\mathrm{I}\mathrm{J}\mathrm{s}^{\backslash }4_{\mathrm{t}}^{:}.i$Fig. 10:$\mathrm{C}\mathfrak{c}>\mathrm{n}$) )$e\mathrm{q},1^{\cdot}\mathrm{i}\mathrm{f};0\iota 1$$\mathrm{r}.$’of FDTD, )$\lfloor$$\mathrm{R}\mathrm{I}^{\supset},$$l.\mathrm{t}\mathrm{K}4\mathrm{C}\mathrm{K}$, opt$M_{3}\cdot\Lambda f:$and$(.)\mathrm{p}\mathrm{t}M_{4}\cdot\Lambda l_{4}$を払いながら以下の表式を数値計算した.
$p=ik \sum_{n=0}^{\infty}(.2n+1)P_{n}(\cos\theta)h_{n}^{(1)}(kr_{>})[i_{n}(kr_{<})-a_{n}’h_{n}^{(1)}(kr_{<})]$,: angle between OS andcoefficient $a_{n}’=j_{n}’(ka)/h_{n}^{(1)’}(ka)$ ,andsmaller ofand .
$\mathrm{O}\mathrm{P}$$\theta$$r_{>}:$$P_{n}$

larger of: Legendre 関数,$r_{OS}$$h_{n}^{(1)}$$r_{OP},$$r_{OS}$$r_{<}:$: 球 Hankel 関数,$j_{n}$$r_{OP}$: 球 Bessel 関数,

$P$

(26)
(27)
(28)

(29)

: 観測点.

代表的な結果を図 13 に周波数領域において表示した. 図 13(右上) より従来型 FDTD と今回の方法 LMM では最大程度の差があること,
下左右の図を比較することによって, LMM では FDTD と比べると音源の$1\cdot 5\mathrm{d}\mathrm{B}$12

Fig. 11:

Fig. 12:$d\epsilon^{j},v_{2}=2_{\lrcorner}^{\tau}|f_{n}-f_{\epsilon:x(xc,t}...|^{2}$)$\mathrm{C}\mathrm{o}\mathrm{t}\mathrm{l}:\iota\cdot \mathrm{a}\mathrm{n}\mathrm{f}|t1\mathrm{t}\mathrm{l}\mathrm{U}11$$\mathrm{e}^{1},\mathrm{r}\mathrm{v}\mathrm{s}$.

vs$()_{\{)}$$r\mathit{4}c_{i}v_{2}$近傍, 球の裏側などで厳密解との差が 1-2dB 小さくなっていることが指摘できる. このことから現実に近い 3 次元計算でも, 従来法よりも LMM と

comPact 差分を組み合わせた方法では改善がみられることがわかった.

4

おわりに

結果を要約すると,
は, 従1) 4 次精度 staggered compact スキームと最適化 3 段階法
来型 FDTD 法より 1 次元のベンチマーク問題に対して 1-2 桁少ない誤差を示した.
$(M_{3\varphi t})$2) 移流方程式に対する最適化多段階法は, 波動方程式に対しても安定領域および最適と考えられる領域にある.
3) 本計算法は Courant 数 0.2 以下で DRP および最小ステンシル幅の 4 次
より正確である.
精度 compact スキームと組みあわせた
4) 本計算法は既存の FDTD code の改良に適している.
ところで, 移流方程式を LMM によって, 続けて 2 回時間を進めると, 波
動方程式を 1 回時間積分したことと同じことになる $(M_{p+q-1}^{2}=M_{p}\cdot M_{q})$ .
そこで, 移流方程式の時間積分に対して最適化された LMM は, 波動方$\mathrm{R}\mathrm{K}4\mathrm{C}\mathrm{K}$

程式の時間積分に対しても最適解を与えるのかという疑問が生じる.

また, ある程度 p, q の値が大きくなると, 波動方程式に LMM を用いる方が13$\mathrm{a}.\iota_{\ovalbox{\tt\small REJECT}_{\theta_{-}^{f}}}.$$1\infty\cdot\underline{\ovalbox{\tt\small REJECT}_{\Re}}.\cdot$.

$.\cdot\overline{.\ovalbox{\tt\small REJECT}^{:}}$$n.0_{\iota_{\phi^{\mathrm{Y}}}^{\theta\S}\ _{\{}},$:$;_{p}.|^{\prec:}\sim:’\searrow\ovalbox{\tt\small REJECT}^{\sim_{-}}.-\ovalbox{\tt\small REJECT} m::$.

$\mathrm{U}.0^{\cdot}\ovalbox{\tt\small REJECT}^{1}\dot{:}\Delta^{=}\prime r_{\aleph}$Soun4$\mathrm{W}$

I

Sound

)

$\mathrm{P}\mathrm{r}\mathrm{o}\epsilon \mathrm{s}\mathrm{u}\mathrm{r}\epsilon 1\mathrm{d}\mathrm{B}$$1000\mathrm{B}\mathrm{z}$

Step
Uu
B. 20
nesh

$\mathrm{n}\uparrow \mathrm{r}\mathrm{u}\cdot$$–$$\mathrm{P}\mathrm{r}\cdot \mathrm{s}\mathrm{s}\mathrm{u}r\epsilon \mathrm{D}\mathfrak{l}2\uparrow*\mathrm{r}\cdot \mathrm{n}\mathrm{c}\epsilon$$1000\mathrm{M}\mathrm{z}$$1\mathrm{t}\mathrm{i}$.

$\mathrm{t}\cdot \mathrm{w}\mathrm{l}$

Wu 1

St

$\mathrm{t}\mathrm{i}$$\mathrm{U}\}\mathfrak{n}i\cdot \mathrm{u}\mathrm{r}\mathrm{n}\epsilon \mathrm{s}\mathrm{h}\in$

1dB)$\epsilon \mathrm{p}-\mathrm{F}\mathrm{D}\mathrm{T}\mathrm{D}(\mathrm{c}\mathrm{e}\mathfrak{n}\mathrm{t}\mathrm{r}*1)$

B. 20$(\mathrm{n}\cdot)-$.
$\mathfrak{g}.0\underline{.\mathrm{B}_{:_{3}}\mathrm{a}K}$$\mathrm{s}.0_{\mathrm{k},-}‘..4_{::}\prime\prime$.
$4.0^{\backslash }l.l.:_{::^{j}}^{\grave{\mathfrak{m}}_{1}}.\S\ovalbox{\tt\small REJECT}_{i_{r}}^{\backslash }\vee|\ovalbox{\tt\small REJECT}^{r}*\backslash ‘$:$0.t^{\mathrm{b}}\triangleleft.".\overline{...\ovalbox{\tt\small REJECT}^{*}s_{i}.\mathrm{R}arrow.\cdot.\theta_{:}h^{\backslash :}::::\|\backslash _{i}^{:::’}}$.

$\mathrm{V}\mathrm{i}n\mathrm{i}..\cdot \mathrm{u}\cdot\cdot\iota \mathrm{a}\mathrm{h}\mathrm{l}000\mathrm{H}2l\mathrm{u}\mathrm{I}\mathrm{t}\mathrm{i}$$\backslash t_{\mathrm{I}\mathfrak{n}.|\cdot 11-\cdot\cdot\epsilon \mathrm{h}\cdot 201\mathrm{r}\cdot)}^{000\aleph \mathrm{a}\mathrm{F}0\tau 0- \mathbb{E}\mathrm{X}l\mathrm{C}\mathrm{I}}$$\Gamma \mathrm{i}\mathrm{g}$.

}: Sound)$\mathrm{p}\mathrm{r}\mathrm{t}_{\text{ノ}^{}\mathrm{Y}}\dot{\mathrm{b}}\mathrm{S}’ 11\mathrm{r}‘..$$(\mathrm{H}7_{\mathrm{J}}).\mathrm{a}\mathrm{I}1(1(1\mathrm{i}\mathrm{f}l\mathfrak{k}^{1},1^{\backslash }\mathrm{C}^{>}\mathrm{n}\mathrm{c}^{\backslash }\mathrm{e})j\mathrm{i}\iota 1p\})\mathrm{e}\mathrm{t}\mathrm{W}’\mathfrak{k}^{1}.,\mathrm{e}\mathrm{n}$witll$b_{0}\simeq^{\vee}$

in

in
FDTD

$\prime \mathrm{d}\mathrm{i}\mathrm{s}’1\iota\cdot \mathrm{i}\mathrm{b}\iota 1\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}$$1_{\dot{\mathrm{i}}}$] $..04$ )(,(($11\mathrm{P}\mathrm{P}^{\mathfrak{k}^{\iota}\mathrm{r}),\Delta p}.$$\rceil_{(^{>},\mathrm{v}\mathrm{i}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}}$$(\mathrm{d}\mathrm{t}\}.)\mathrm{f}’()\mathrm{r}\mathrm{I}^{\neg}\{\mathrm{D}\mathrm{T}\mathrm{D}+\mathrm{C}\mathrm{D}$\iota 1\tau ).\iota ‘’ ノ\downar ow .

.

$\mathrm{t}\epsilon \mathfrak{p}- \mathbb{E}\mathrm{X}\mathrm{A}\mathrm{G}78.20(\cdot\cdot)$‘.

.:$\mathrm{t}\mathrm{h}\mathrm{c},$$\mathrm{f}\mathrm{r}.\mathrm{e}\mathrm{q}\mathrm{u}\mathfrak{c}_{d}^{\mathrm{j}}1\mathrm{J}\mathrm{C}\mathrm{y}$donlain at$1(\}\circ()$$\mathrm{a}\mathrm{r}\mathrm{l}(1\mathrm{p}\mathrm{r}\mathfrak{k}^{1}.,\iota 5’\mathrm{C}^{\Delta}11\mathrm{t},(\Lambda^{J}I_{3}’.\cdot\Lambda I_{3}.+\mathrm{C}\mathrm{D}\mathrm{t}11,.\backslash ^{\backslash }4$from$i_{-}1\mathrm{I}\mathrm{t}\mathrm{d}\iota$))$\mathrm{t}\cdot,\}1\epsilon_{d}$$1^{\alpha}\mathrm{C}^{\iota}.\dot{\mathrm{b}}’‘!.\iota\iota\uparrow,$$‘ 1,\mathrm{x}^{r}\mathrm{a}\mathrm{c}^{\backslash }\mathrm{t}\mathrm{s}^{\tau}\mathrm{o}\mathrm{l}\mathrm{l}\mathrm{l}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}$.

) hown$..\mathrm{s}$$\mathrm{n}\iota\in_{J}1\mathrm{t}\mathrm{l}\mathrm{l}\mathrm{O}(1(1\mathrm{c})\mathrm{W}\mathfrak{k}_{d}^{1}1)$

計算時の容量が少なくてすむ. このように中間変数 V を使用する方法と,速度を使わずに p(または ) についての波動方程式を直接解く方法のどちらが有利かについては不明である. さらに安定性限界と空間スキームの$\rho$

実効波数を最適化する際考慮に入れる必要がでてくる.
そこで今後の課題として,

1) 波動方程式に対する線形多段階積分法 (LMM) の最適化
2) staggered compact スキームの最適化 (anistrophy, eigen values, 境界条
件)3) 空間スキームの実効波数を考慮にいれた最適化法を検討する予定としている.
検討の結果として, 時間発展の実効角振動数を多項式で近似する Rungeほどの広い絶対安定領域は望めないが, DRP 的な波数の再Kutta 法$(\mathrm{R}\mathrm{K})$14
現性をもつスキームと組み合わせたときに, ある程度の Courant 数の範囲においては既存の 4 次精度 Runge-Kutta 法よりは正確な時間積分法をみつけることが可能であると予想される. このような最適化された LMM
は, Runge-Kutta 法よりは nlemory を使うが, flux の計算に時間のかかるよりも速く計算が可能である.
ような場合には$\mathrm{R}\mathrm{K}$参考文献(1) Yee, K.S. “ Nllmerical solution of initial bollndary value problems involving Maxwell’s equations in isotropic media,” IEEE Ransactions on
and propagation, vol.AP-14 no.3, pp.802-807, 1966.
(2) Tam, C.K.W. and Webb, J.C. “ Dispersion-relation-preserving schemes
for complltational acoustics,” J. Comput. Phys. vol.107, pp.262-281,$\mathrm{A}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{n}\mathrm{n}\mathrm{a}_{\iota}\mathrm{s}$1993.
(3) Lele, S.K., “ Compact finite difference schemes with spectral-like resolution,” J. Compllt. Phys., vol.103, pp.16-42, 1992.
Workshop on benchmark(4) Hardin, J.C. et al. (Eds.) “
problems in complltational aeroacollstics,” NASA CP-3300, 1995.
(5) Shampine, L.F. “Numerical solution of ordinary differential equations,” Chapman&Hall (New York), pp.187-198, 1994.
(6) Fu, F.Q., Hllssaini, M.Y. and Manthey, J.L. “ Low-dissipation and lowdispersion Rllnge-Klltta schemes for comsutational acollstics,” J. Compllt. Phys., vol.124, pp.177-191, 1996.
(7) Carpenter,M.K. and Kennedy, C.A., “ Fourth-order -storage RungeKutta schemes,” NASA TM 109112, 1994.
(8) 1. Iwatsll, R., Tsllrll, H. and Hirosawa, K., “Application of optimized
compact finite difference schemes on uneven grid to the complltation of
acoustic wave propagation,” proc. The Thirteenth International Congress
pages) Vienna, Allstria, Jllly
on Sound and Vibration, CD-ROM(
2-6, 2006.
$\mathrm{I}\mathrm{C}\mathrm{A}\mathrm{S}\mathrm{E}/\mathrm{L}\mathrm{a}\mathrm{R}\mathrm{C}$$2\mathrm{N}$$\mathrm{S}\mathrm{S}38arrow 3,8$

「複雑流体の数環とシミュレーション」 RJMS 研究集会 (京大数理解析研宍所) 11 月

$16\cdot 1\mathit{7}\mathrm{B},$$2\mathfrak{m}6$

年