数理解析研究所講究録
第 1568 巻 2007 年 127-134

127

Green 関数・DE-Sinc-Nystr\"om 法を用いた
Sturm-Liouville 型固有値問題の数値解法
久保隆貴
TAKAKI KUBO*

筑波大学数理物質科学研究科
GRADUATE SCHOOL OF PURE AND APPLIED SCIENCES, UNIVERSITY OF TSUKUBA

1

はじめに二重指数変換 (double exponential transformation, 略して DE 変換) は, 解析関数の定積分を高精度に能率良く計算するために, 1974 年に高橋秀俊・森正武によって提案された変換である [9].
この変換によって導かれる二重指数公式 (double exponential formula, 略して DE 公式) は, 最も少ない分点数で最も高い精度を達成するという意味で最適な公式 [8] であることが知られており,土木工学や金融工学に至る広い分野で利用されるようになっている. 一方, 近年の二重指数変換に関する著しい研究成果の一つに, 1980 年代に F. Stenger らによって構築された Sinc 関数近似の理論 [5] への応用があり, Sinc 関数近似の場合にも二重指数関数型の減衰がある意味で最適性を実現することが 1996 年に杉原正顕によって証明された [7]. この Sinc 関数近似に対する研究は,
Sinc-Galerkin 法および Sinc-collocation 法といった微分方程式の数値解法, 数値不定積分や不定畳み込みなどの近似公式を, より高性能なものへと進化させることになった. 本稿では, 二重指数変換に基づく数値 Green 関数法 [1] と $Sinc- Nystr\ddot{o}m$ 法 [2] とを組み合わせた Sturm-Liouville型固有値問題に対する数値解法を提案する. なお, 二重指数変換とその応用は計算代数の分野では馴染みがうすいので, 下記ではやや長い解説をすることを御容赦願いたい.

2

数値不定積分2.1

不定積分に対する二重指数公式

DE 変換は定積分だけでなく不定積分の数値計算にも適用ができる. その基本は区間における関数近似であり, 基底関数として使用する Sinc 関数は次の形で定義される.
$S(j,h)(t)= sinc(\frac{t}{h}-j)=\frac{\sin\frac{\pi}{h}(t-jh)}{\frac{\pi}{h}(t-jh)}$Sinc 関数 (2.1) のグラフを図 1 に示す.
luboimath.tsukuba. ac.jp$(-\infty, \infty)$

(2.1)

128

図 1: Sinc 関数 sinc $(t/h-j)=S(j, h)(t)$数値不定積分の代表的な例として, 1 変数解析関数 $f(x)$ の不定積分$I(s)= \int_{a}f(x)dx$,$a<s<b$(2.2)

に対して DE 公式を示す. 被積分関数 $f(x)$ は $a<x<b$ に特異性を持たないと仮定するが, 端点$x=a,$$b$に特異性を持っていてもよい. 不定積分 (2.2) に対して DE 変換$x= \phi(t)=\frac{a-b}{2}$tanh$( \frac{\pi}{2}$sinh $t)+ \frac{a+b}{2}$(2.3)

を適用すると, 次のようになる.
$I(s)= \int_{a}^{\iota}f(x)dx=\int_{-\infty}^{\tau}u(t)dt$$u(t)=f(\phi(t))\phi’(t))$,$\tau=\phi^{-1}(s)$

(2.4)

.

無限区聞の不定積分 (2.4) の被積分関数 $u(t)$ が, 複素平面の帯状領域 $|{\rm Im} t|<d$ で正則であり,かっとをある正の定数として$\alpha’$$c$1を満たすとする. 刻み幅$h$($u(t)|\leq\alpha’$ $\exp$(-cexp 1 $t|$ )),$|t|\prec\infty$(2.5)と打ち切り項数 $N$ を次式を満たすように決める.
$h= \frac{1}{N}$log$( \frac{\pi dN}{\alpha})$(2.6)( は任意の小さい正数) である. (2.4) の被積分関数を Sinc 基本展開すると,不定積分に対する次の DE 公式が得られる [4].
ただし,$\alpha=\alpha’-\epsilon$$\epsilon$$I_{Dg}(s)=h \sum_{j=-N}^{N}f(\phi(jh))\phi’(jh)(\frac{1}{2}+\frac{1}{\pi}$

Si $( \pi\frac{\phi^{-1}(s)}{h}-\pi j))+o(\exp(-\frac{\pi dN}{\log(\pi dN/\alpha)}))(2.7)$

129
ただし, Si は次式に示される積分正弦関数である.

Si

2.2$(t)= \int_{0}^{t}\frac{\sin\tau}{\tau}d\tau$

Si $(-t)=-SI(t)$

,

(2.8)

積分正弦関数の数値計算法

積分正弦関数 (2.8) について, $t=0,$$\infty$のときは Dirichlet の積分定理から

Si (0) $=0$ , Siを得る. また, sinSi$\tau$(2.9)$( \infty)=\frac{\pi}{2}$を Taylor 展開して (2.8) に代入すると(2.10)$(t)= \int_{0}^{t}\frac{1}{\tau}(\tau-\frac{\tau^{3}}{3!}+\frac{\tau^{5}}{5!}-\frac{\tau^{7}}{7!}+\cdots)d\tau=\sum_{k=1}^{\infty}\frac{(-1)^{k-1}t^{2k-1}}{(2k-1)(2k-1)!}$を得る. ただし, 実際の計算では無限和の上限を $k=n$ で打ち切る必要がある. (2.10) の右辺のが小さいときは収束が速いので直接計算によって高精度に計算できるが,が大きいときは次の性質を利用して計算する必要がある.
積分正弦関数 Si は (2.9) から

級数は,

$t$

$t$

Si

$(t)= \int_{0}^{t}\frac{\sin\tau}{\tau}d\tau=\frac{\pi}{2}-\int^{\infty}\frac{\sin\tau}{\tau}d\tau$(2.11)と表現できる. この右辺の第二項を漸近級数展開し, 無限和を積分表示で表すと

Si $(t)= \frac{\pi}{2}-f(t)$ cos $t-g(t)$ sin ,(2.12),(2.13)$t$$f(t)= \int_{0}^{\infty}\frac{te^{-x}}{x^{2}+t^{2}}dt$$g(t)= \int_{0}^{\infty}\frac{xe^{-x}}{x^{2}+t^{2}}dt$を得る. (2.13) は定積分であるから, 半無限区間の DE 公式(2.14)$t=\phi(u)=\exp(u-\exp(-u))$を採用し, 等間隔刻み幅$h$の台形公式を適用すれば, 高精度に計算することができる.

3

Sturm-Liouville 型固有値問題と Green 関数

3.1

Sturm-Liouville 型固有値問題$L[u]$

を Sturm-Liouville 作用素 $\frac{d}{dx}(p(x)\frac{du(x)}{dx})+q(x)u(x)$ とし, 微分方程式$\frac{d}{dx}(p(x)\frac{du}{dx})+q(x)u+\lambda\rho(x)u=0$を考える. ただし,$\lambda$,$a<x<b$(3.1)

は固有値で, 境界条件として次式を与える.

$u(a)=0$ ,$u(b)=0$(3.2)

130境界条件を満足する微分方程式の解 $u(x)$ と固有値とを求める問題は Sturm-Liouville 型固有値問題である. $\rho(x)$ は荷重関数と呼ばれ, 連続実関数で $\rho(x)>0$ である. 荷重関数は直交条件に際して重要となるが, $\rho(x)=1$ の場合が多いため, 本稿ではこれを採用する. 境界条件 (3.2) はDirichlet 境界条件と呼ばれ, 微分方程式の解 $u(x)$ が両端 $x=a,$ で固定されていることを意味$\lambda$

$b$

している.

3.2

Sturm-Liouville 作用素 $L[u]$ には Green 関数の存在が知られている.

Green 関数

Sturm-Liouville 型固有値問題では固有値 が特定の値をもつときのみ解が存在し, 解は固有関数と呼ばれる. 解は Green 関数を用いて表現できることが知られており, 今の場合の Green 関数$\lambda$$G(x,\xi)$は作用素 $L[u]$ に対する同次壇界値問題$\{\begin{array}{ll}L[u]=0 y(a)=0, y(b)=0\end{array}$(3.3)

に関して, 次の 3 条件を満たす関数である.
$\bullet$$\bullet$$\bullet$は境界値問題 (3.3) を満たす連続関数である.
$x\neq\xi$ では $G(x,\xi)$ のに関する 1 階, 2 階導関数が連続である.
$G(x,\xi)$ は $x=\xi$ で次を満たす.
$G(x,\xi)$$x$$( \frac{dG}{dx})_{x\prec+0}--(\frac{dG}{dx})_{x=\xi-0}=-\frac{1}{p(\xi)}$

上記の条件を満たす Green 関数は, 境界値問題 (3.3) の 2 つの基本解 $y_{1}(x),y_{2}(x)$ を用いて$G(x,\xi)=\{\begin{array}{ll}y_{1}(\xi)y_{2}(x) (\xi\leq x)y_{1}(x)y_{2}(\xi) (x\leq\xi)\end{array}$(3.4)と表現できる [11]. $y_{1}(x),$ $y_{2}(x)$ はそれぞれ $x=a,$ での境界条件を満たす. 固有値問題 (3.1) の解は $G(x,\xi)$ を用いて次の積分方程式の形で表せる (Green の定理) .
$b$$u(x)= \lambda\int_{a}^{b}G(x,\xi)u(\xi)d\xi$

実際の計算では,(3.5)

Green 関数の $x=\xi$ における特具性を避けて解析性を保っために$u(x)=\lambda\{\int_{a}^{\xi}G(x,\xi)u(\xi)d\xi+\int_{\xi}^{b}G(x,\xi)u(\xi)\mathfrak{X}\}$(36)として積分を分けて計算する必要がある. Sturm-Liouville 作用素から得られた Green 関数 $G(x,\xi)$は完全な積分核であることが知られており, 量子的に飛び飛びな無限個の固有値と固有関数をもつことがわかっている.

131ta4$DE-Sinc$-Nystr\"om4.1

二重指数変換に基づ \langle Green 関数法

Sturm-Liouville 型固有値問題の解を表す積分方程式 (3.6) に二重指数変換に基づいた数値不定積分公式を適用する. 右辺の$\frac{1}{\lambda}u(x)=$$v_{1}(t)$$\lambda$を左辺に移項した式に対して適用すれば,$y_{2}(x)h \sum_{j=-N}^{N}v_{1}(jh)\sigma_{+}(x;j)+y_{1}(x)h\sum_{j\approx-N}^{N}v_{2}(jh)\sigma_{-}(x;j)+E_{N}$$y_{i}(\phi(t))u(\phi(t))\phi’(t)$$=$を得る [1]. ただし, (4.1) の$\phi(t)$,(4.1),

は不定積分に対する二重指数変換$\xi=\phi(t)=\frac{b-a}{2}$tanh$( \frac{\pi}{2}$(4.2)

sinh $t)+ \frac{b+a}{2}$であり. \mbox{\boldmath $\sigma$}\pm (x; のは積分正弦関数 Si $(t)$ を用いて次式で表される.
$\sigma\pm(x;j)=(\frac{1}{2}\pm\frac{1}{\pi}$

刻み幅

$h$

Siと片側求解点数 $N$ を (2.6) を満たすように決めるとき, 誤差項$E_{N}=O( \exp(-\frac{\pi dN}{\log(2dN/\alpha)}))$

また,

$x$

(4.3)

$( \frac{\pi\phi^{-1}(x)}{h}-\pi j))$,$\alpha$と$d$$E_{N}$

は正の定数.

(4.4)

が Sinc point $x_{k}=\phi(kh)$ であれば $\sigma\pm(x_{k};j)$ は次式のように簡単になる.

(4.5)$\sigma\pm(x_{k};j)=Si(\pi(k-j))$4.2

は次式となる.

Nystr\"om の方法への遁用

積分方程式 (3.5) において固有値$G$$\lambda$を改めて $1/\nu$ とおけば, 積分作用素; $u(x) rightarrow(Gu)(x)=\int_{a}^{b}G(x,\xi)u(\xi)\mathfrak{X}$$a\leq x\leq b$(46)に対応する固有値問題は次式となる.
$(Gu)(x)= \int_{a}^{b}G(x,\xi)u(\xi)d\xi=\nu u(x)$,$a\leq x\leq b$(4.7)

が積分作要素 $G$ の固有値と固有ベクトルであるならば, 関数解析の適当な枠組みを用いて,$(1/\nu,u)$ は随伴する微分形の問題の固有値と固有ベクトルであることが容易に確かめられる. 積分方程式 (4.7) を数値積分公式によって近似する解法は, Nystr\"om の方法 [10] と呼ばれ, 本研究では数値積分公式として二重指数公式を採用する.と $N$ が (2.6) を満たすとして二重指数変換 (4.2) を適用する. 変換 (4.2) の変(4.7) に対し,数 の変域 $(-\infty, \infty)$ で等間隔に並ぶ Sinc point$(\nu,u)$

$h$

$t$

$t_{k}=kh$,$k=-N,$ $-N+1,$$\cdots$

, $N$

(4.8)

132をとれば,$x$の変域 $(a, b)$ における標本点は, 変数$x_{k}= \phi(t_{k})=\frac{b-a}{2}$となる.

tanh

$x$

での対応する Sinc point$( \frac{\pi}{2}\sinh kh)+\frac{b+a}{2}$

,

(4.9)

$kh=\phi^{-1}(x_{k})$この Sinc point を二重指数公式 (4.1) の標本点として選べば,右辺の有限和とちょうど同じ分点数を持つことになり,$\{u(x_{k})\}$

二重指数公式の標本点がに関して$\nu u(x_{k}\ovalbox{\tt\small REJECT} y_{2}(x_{k})h\sum_{j=-N}^{N}v_{1}(jh)\sigma_{+}(x_{k};j)+y_{1}(x_{k})h\sum_{j=-N}^{N}v_{2}(jh)\sigma_{-}(x_{k};j)+E_{N}$(4.10)を得る. この等式から行列 $M$ , ベクトル $u$ をそれぞれ$M_{kj}=h\{y_{2}(x_{k})y_{1}(x_{j})\sigma_{+}(x_{ki}j)+y_{1}(x_{k})y_{2}(x_{j})\sigma_{-}(x_{ki}j)\}\phi’(jh)$$u$,$=(u(x_{-N}),u(x_{-N+1}),$ $\cdots,u(x_{N}))^{t}$$k,j=-N,$ $-N+1,$ $\cdots N$(4.11)と定義すれば, (4.10) は次式となる.
$Mu=\nu_{m}u$(4.12) を満たす $2N+1$ 個の固有値,(4.12)$m=1,$ $\cdots 2N+1$と固有ベクトル $u$ は, 既存の線形計算の算法によって計算することができる. 一つの固有値について, $x\neq x_{k}$ を含む任意のに対する固有ベクトル $u_{m}(x)$の値, すなわち Sturm-Liouville 型固有値問題 (3.1) の解 $u(x)$ は, $\nu_{m}\neq 0$ という条件の下では、
Nystr\"om の自然補間公式を用いて, 算式$\nu_{m}$$x$$u_{m}(x)= \frac{h}{\nu_{m}}\{y_{2}(x)\sum_{j=-N}^{N}v_{1}(jh)\sigma_{+}(x;j)+y_{1}(x)\sum_{j=-N}^{N}v_{2}(jh)\sigma_{-}(x;j)\}+E_{N}$で計算することができる. これは微分形の固有値と積分作用素の固有値
Groen 関数法に用いる二重指数公式 (4.1) に他ならない.
$\lambda$

5$\nu$(4.13)の関係を考えれば,数値例
簡単な熱伝導方程式を変数分離法によって解く場合, 常微分方程式(5.1)$\frac{d^{2}u}{dx^{2}}=-\lambda u$の形の Sturm-LioudUe 型固有値問題がよく現れる. ここで, 境界条件は次式で与えられる.

(5.2)$u(O)=u(1)=0$

微分方程式 (5.1) に対し, 前節の DE-Sinc-Nystr\"om 法を適用して数値実験を行った. ただし,と $N$ の関係を (2.6) とし, $\alpha=1,$ $d=\pi/2$ と決めた. (4.12) の行列は非対称行列となることから, 固有値の計算にはべき乗法を用いた. また, $L[u]=u”$ と選べば, $L[u]$ に対応する Green関数を構成する基本解 $y_{1}(x)$ と $y_{2}(x)$ は次式となる.
$h$$M_{kj}$$y_{1}(x)=x$,$y_{2}(x)=1-x$

(5.3)

133
固有値問題 (5.1) を積分作用素 (4.7) の形とした場合の固有値$\nu_{m}=\frac{1}{m^{2}\pi^{2}}$である.

したがって,

この$\nu$,$m=1,2,$$\nu$の真の解は(5.4)$\cdots$に随伴する微分形の問題 (5.1) の固有値$\lambda_{m}=m^{2}\pi^{2}$,$m=1,2,$$\lambda$の真の解は(5.5)$\cdots$の誤差の減衰を図 2 に示す. 次に, DB
Sinc-Nystr\"om 法における Nystr\"om の自然補間公式 (4.13) を用いて, 3 つの固有値 $\lambda_{m}(m=1,2,3)$に関する (5.1) の数値解を計算した. 正規化した真の解はとなる. 大きい方から 5 つ $(m=1,2,3,4,5)$ の固有値

$u_{m}(x)=\sqrt{2}$sin $m\pi x$ ,$\lambda_{m}$$m=1,2,$である. 区間 $[0,1]$ を $1W$ 等分した点, すなわち $x=001,$$\cdots$$\cdots$

099 で数値解を計算した.

(5.6)
この99 点での真の解との絶対娯差の最大値を図 3 に示す.

N
図 2: 固有値の誤差の減衰$N$

図 3: 数値解の最大誤差の減衰

参考文献[1] T. Echigo, M. Mori, Numerical Green’s Function Method Based on The DE Transformation,
J. Comput. Appl. Math. 23 (2006) 193-205.

[2] B. Bialecki F. Stenger, Sinc-Nystr\"om method for numerical solution of one-dimensional
Cauchy singular integral equations given on a smooth arc in the complex plane, Math.
Comp. 51 (1988) $13\succ 165$ .
[3] M. Muhammad, A. Nurmuhammad, M. Mori and M. Sugighara, Numerical solution ofintegral equation by means of the Sinc collocation method based on the double exponential
transformation, J. Comput. Appl. Math. 177 (2005) 269-286.

[4] M. Muhammad, M. Mori, Double exponential formulas for numerical indefinite integration,
J. Comput. Appl. Math. 161 (2003) 431-448.

134[5] F. Stenger, Numerical Methods Based on Sinc and Analytic lPttrnctions, Springer-Verlag,
Berlin, New York, 1993[6] J. Lund, K.L. Bowers, Sinc Method for Quadrature and Differential Equations, SIAM,
Philadelphia, PA, 1992
[7] 杉原正顕, 二重指数関数型変換を用いた Sinc 関数近似, 京都大学数理解析研究所共同研究集
会「科学技術における数値計算の理論と応用 」 予稿集 (1996) 26-27$\mathbb{I}$[8] M. Sugihara, Optimality of the double exponential formula–functional analysis approach-, Numer. Math. 75 (1997) 379-395.
[9] H. Takahasi, M. Mori, Double exponential formula for numerical integration, Publ. RIMS
Kyoto Univ. 9 (1974) 721-741.
[10] E. J. Nystr\"om, Uber die praktische Aufl\"osung von Integralgleichungen mit Anwendungenauf Randwertaukaben, Acta Math. 54 (1930) 185-204.
[11] R. Courant, D. Hilbert, Method of Mathematical Physics Vol.1, Interscience Publishers INC.
New York, 1953