重み付高精度DE-Sinc近似アルゴリズムとその応用

全文

(1)Vol.2009-HPC-120 No.5 2009/6/12. 情報処理学会研究報告 IPSJ SIG Technical Report. の帯状領域. Dd := {z ∈ C| |Imz| < d}. 重み付高精度 DE-Sinc 近似アルゴリズムとその応用. (1). (d は正数) において解析的な関数を対象とする．数値計算において，関数近似公式は理論・実用両面でも最も基本的である．主な関数近似. アヒニヤズ. ヌルメメット†1. 公式としては多項式補間やスプライン補間，Fourier 多項式による近似公式があるが，近年. Sinc 関数を用いた Sinc 関数近似が注目されている1),2),4)–6) ．Sinc 関数は次のような基底本稿では，従来の DE-Sinc 近似に重みをつけることで，境界条件がゼロではない関数においても高精度の近似が得られ，積分方程式等の近似解にもそのまま適用できることを報告する．この近似法は 2 変数関数の近似においても有効である．. 関数である．. sin[π(t − jh)/h] , [π(t − jh)/h] ここで，S(j, h)(t) は，sinc 関数 ( S(j, h)(t) =. sinc(t) =. Highly Accurate DE-Sinc Algorithm with Weight Function and It’s Applications. sin(πt) πt. j = 0, ±1, ±2, · · ·. (2). x 6= 0. 1 x=0 を用いて S(j, h)(t) = sinc(t/h − j) とも書ける．R 上の関数 f (t) を，R 上に等間隔にとった値によって近似する無限和. Ahniyaz Nurmuhammad†1. f (t) ≈. ∞ X. f (jh)S(j, h)(t). (3). j=−∞. 有限で打ち切った次の式が Sinc 関数近似である．. In this paper, we will report about applications of the DE-Sinc approximation with weight function. And this method also give highly accuarate numerical approximate for the function of two variables.. fN (t) ≈. N X. f (jh)S(j, h)(t).. (4). j=−N. Sinc 関数に基づく数値計算法は Sinc 法と呼ばれる．C(f, h) が通信理論・信号処理への工学的応用目的で利用され，Shannon のサンプリング定理等がある．それから 30 年程の研究. 1. はじめに. によって，(4) に基づく Sinc 法は，関数近似，導関数の近似，定積分および不定積分の近. 1.1 Sinc 近似法. 似，常微分方程式および偏微分方程式の解の近似に適用されるようになった1)–6) ．Sinc 法. 近年，計算機の性能の上向により，理工学をはじめとする様々な分野で数値計算による模. が注目を集める理由の一つは，他の補間を用いた手法に比べてはるかに高精度であることに. 擬実験が多く行われるようになり，数値計算の精度に対する要求も高まっている．その数値. ある．(4) に現れる関数評価回数 N 対し，関数 f の帯状領域での解析性の下で，無限遠点. 計算の対象として，初等的な微積分学の範囲での数学的計算は多くの割合を占めている．微. での減衰の仕方によって，Sinc 法は大体 O(exp(−c0 N/ log(c00 N ))) の精度を実現すること. 積分学は連続変数の関数を対象とするが，このような関数を対象として数値計算を行う場. が知られている．この精度を実現するのに関数が無限遠点で次のような減衰することが条件. 合には，有限個の数値によって対象の関数を近似し，これらの数値を用いた計算結果から，. になる．. 既成の適当な理論に基づいて微積分などの演算結果に関する情報を得ることになる．このた. ∀t ∈ R : |f (t)| ≤ β exp(−α exp(γ|t|)). めには，対象の関数が適当な領域で解析的あることが必要になる．本報告では，複素平面上. (5). この (5) が二重指数 (DE) 的減衰と呼ばれる．. 1.2 DE-Sinc 近似法近似式 (4) は R 上の関数に対するものであり，関数が帯状領域で解析的で，かつ実 †1 東京電機大学 (Tokyo Denki University). 1. c 2009 Information Processing Society of Japan °.

(2) Vol.2009-HPC-120 No.5 2009/6/12. 情報処理学会研究報告 IPSJ SIG Technical Report. 軸上では (5) のような条件を満たすとき，前述のような高精度が実現出来る．一般には，. 重み関数 Φ(c − t) として DE 変換の微分 cosh(c − t) ¡ ¢ |Φ(c − t)| = 2 π cosh µ 2µsinh(c − t) ¶¶ π(1 − ²) exp |c − t| as t → ±∞ = O exp − 2 を適用し，c = t とすると，関数 f (x) に対する次の Sinc 近似式が得られる．. 関数 f : (a, b) −→ R に対し (4) に基づいた高精度近似を実現するには，適当な変換. ψ : R −→ (a, b) を用いて f を変換したものがそれらの条件を満たせばよい．この方針に基 8) づいて，ここで用いるのは次の二重指数変換 (DE Ã r変換) である．. ³. ´. ³. ´2. !. 1 1 x−a x−a log + log +1 , (6) π b−x π2 b−x ³ ´ (b + a) (b − a) π x = ψ(t) = φ−1 (t) = tanh sinh t + . (7) 2 2 2 有限区間 (a, b) で定義された関数 f (x) が帯状領域 Dd で解析的で，かつ端点では t = φ(x) = log. j=dφ(x)/he+N. N X. fN (x) ≈. f (ψ(jh))S (j, h) ◦ φ(x).. f (ψ(jh))S(j, h) ◦ φ(x). j=bφ(x)/hc−N. cosh(φ(x) − jh) ¡ ¢ cosh2 π2 sinh(φ(x) − jh). (11). この近似式において，刻み幅を. log(πdN/α) , (12) N のようにとると（α は関数の減衰によって定まる），Sinc 近似式 (11) に基づく数値計算で h=. f (a) = f (b) = 0 であるならば，この DE 変換を用いた Sinc 近似公式は次のようになる．. X. f (x)N ≈. (10). は O(exp(−c0 N/ log(c00 N ))) の精度を実現することが出来る．. (8). このように関数の減衰を加速させる手法によって得られる Sinc 近似式としては，重み関. j=−N. 数に Gauss 関数を用いた Sinc-Gauss 近似もある7) ．. この公式に基づく数値計算法は DE-Sinc 法と呼ばれ，公式を有効にするための条件 (5) は. ·. j=dφ(x)/he+N. f (ψ(t)) に課せられることになる．. X. f (x) =. f (ψ(jh))S(j, h) ◦ φ(x) exp −. j=bφ(x)/hc−N. 2. 1 変数関数における重み付 DE-Sinc 近似. (φ(x) − jh)2 . 2r2 h2. ¸ (13). この近似式を用いた関数近似では O(exp(−c0 N )) 程度の精度が実現できるとされているが，. 関数 f (x) に対する条件 (5) は，Sinc 近似式 (refeqn:DE-sinc) の打ち切り誤差を小さく. 必要精度を実現するのに刻み幅を小さく取る必要があり，基底関数 (サンプリング点) の数. |f (ψ(jh))| よって評価すれば. が dφ(x)/he + N によって定められるため，サンプリング点の数が大幅に増えてしまうこ. 良いのである．一般には，関数 f (x) が区間 (a, b) の端点で f (a) = fa 6= 0, f (b) = fb 6= 0. とになる．Sinc-Gauss 近似を用いて積分方程式等の数値解を求めるときには，最適である. するためのもであり，この場合には打ち切り誤差を. P. |j|>N. のときに，変換後の関数 f (ψ(t)) が減衰条件 (5) を満たさないため，前述の制度が実現で. DE 積分公式を適用すると，積分公式の刻み幅と Sinc-Gauss 近似の刻み幅の統一性が出来. きなくなる．このような場合には，Sinc 近似式 (8) に適当な補助関数 wa (x), wb (x). なくなる．それに，Sinc-Gauss 近似を微分方程式等の数値解法に用いるときには，近似式での基底関数の数が各点に依存して dφ(x)/he + N によって一定に定まらないため，予め近. wa (a) = 1, wa (b) = 0 , wb (a) = 0, wb (b) = 1. 似式での基底関数の数を決めておき，. を加えて，次の近似式が用いられた．. N X. fN (x) ≈ fa wa (x) +. X. j=N. f (x) =. f (ψ(jh))S (j, h) ◦ φ(x) + fb wb (x).. j=−N. j=−N. ·. (φ(x) − jh)2 f (ψ(jh))S(j, h) ◦ φ(x) exp − . 2r2 h2. ¸ (14). のような近似式を用いると，前述の O(exp(−c0 N )) 程度の精度が得られなくなる．一方，. しかし，その適当な補助関数 w(·) の自由度があり，さらに関数の端点での減衰が Sinc 近似の刻み幅に影響するため，最適な補助関数を定めるのが実際の応用で難しいのが現実であ. 我々が提案した重み付 DE-Sinc 近似では，近似式での基底関数の数を予め決めておいても，. る．そこで，変換後の関数にある重み関数 Φ(c − t) (c はパラメ−タ, Φ(0) = 1) をつけて，. 上述のような問題は起こらない．したがって，実際の応用では，次の式を使用することに. 関数 f (ψ(t))Φ(c, t) が減衰条件 (5) を満たすようにして Sinc を近似する．. する．. f (ψ(t))N Φ(c − t) ≈. X. f (ψ(jh))S(j, h)(t)Φ(c − jh). f (x)N ≈. (9). X. j=N. j=dt/he+N. j=−N. j=bt/hc−N. 2. f (ψ(jh))S(j, h) ◦ φ(x). cosh(φ(x) − jh) ¡ ¢ cosh2 π2 sinh(φ(x) − jh). (15). c 2009 Information Processing Society of Japan °.

(3) Vol.2009-HPC-120 No.5 2009/6/12. 情報処理学会研究報告 IPSJ SIG Technical Report. 表 2 例 1. における (11) と (13) による最大誤差，刻み幅，基底関数の数 DE-Sinc with weight Gauss-Sinc. 3. 数値例. N 2 4 8 16 32 64 128. この節では，次の関数を前節で述べたいくつの Sinc 近似式を用いて近似する． √ 例 1. f (x) = 1 − x2 , x ∈ (0, 1) この関数は Dπ/2 で解析的ではあるが，f (0) = 1 6= 0 であるため，二重指数的減衰条件を満たしていない．しかし，関数 f (ψ(t))Φ(c − t) は二重指数的に減衰する．Φ(c − t) は (10) で定義したものである．比較のために，この関数を近似式 (8)，(11)，(13)，(15) によって. maxerror 2.73e-02 6.03e-03 1.90e-04 1.52e-07 2.67e-13 2.13e-22 3.75e-27. h 5.72e-01 4.60e-01 3.16e-01 2.02e-01 1.22e-01 7.20e-02 4.14e-02. φ(x)/h 2 3 4 7 12 21 36. maxerror 2.53e-02 8.09e-04 1.09e-06 2.72e-12 2.36e-23 1.19e-30 2.47e-30. h 4.95e-03 2.48e-03 1.24e-03 6.19e-04 3.09e-04 1.55e-04 7.73e-05. φ(x)/h 308 617 1235 2471 4942 9884 19769. 近似し，(0, 1) 区間を 1000 当分した各点 (端点を除く) si = i/1000, i = 1, 2, · · · , 999 でのこの表 2 からも見て取れるように，近似式 (11) では φ(x)/h は最大でも 36 で，表 1 で. 最大誤差. max |f (si ) − fN (si )|. 1≤i≤999. の結果からはそれは無視できることが見える．それにに対して，同じ程度の精度を実現す. p. N/π, h = (1 − ²)/N, ² =. るためには近似式 (13) ではそれが 9884 となっている上，刻み幅も前者よりはるかに小さ. 1/100 とするが，他の近似式では (12) に従い，そこで α = 1 とした．全ての計算は 4 倍. くなっている．後者の場合には，実際の応用ではそれを無視する事が出来ないことがうかが. 長精度で行われた．. える．. を調べる．近似式 (13) において，Tanaka. 7). にしたがって r =. 重みつき DE-Sinc 近似の精度と効率を確かめるために，従来の DE-Sinc 近似式 (8) と. 4. 2 変数関数における重みつき DE-Sinc 近似. φ(x)/h を無視した重みつき DE-Sinc 近似式 (15) による計算結果を表 1 に示す．φ(x)/h. 2 変数関数 F (x, y) が有界領域 D = (a, b)×(c, d) で必要な解析性を満たすならば，F (x, y). を無視した (14) による近似は収束しないの表 1 には示していない．. がその領域の境界でゼロではなくても重み関数をつけた DE-Sinc 近似によって 1 変数関数 N 2 4 8 16 32 64 128. 表 1 例 1. における DE-Sinc with weight maxerror 7.22e-01 4.51e-02 1.32e-03 1.44e-05 1.60e-09 2.35e-17 3.75e-27. (8), (15) による最大誤差，刻み幅 Standard DE-Sinc h maxerror 5.72e-01 7.38e-02 4.60e-01 3.92e-02 3.16e-01 1.98e-02 2.02e-01 9.94e-03 1.22e-01 4.97e-03 7.20e-02 2.49e-03 4.14e-02 1.24e-03. 同様の精度の次の近似式が得られる．. j=dφ(x)/he+Nx k=dφ(y)/he+Ny. h 5.72e-01 4.59e-01 3.16e-01 2.02e-01 1.22e-01 7.20e-02 4.14e-02. X. FNx ,Ny (x, y) ≈. X. F (ψ(jh), ψ(kh))S(j, h) ◦ φ(x)S(k, h) ◦ φ(y). j=bφ(x)/hc−Nx k=bφ(y)/hc−Ny. ×. cosh(φ(x) − jh) cosh(φ(y) − kh) ¡π ¢ ¡ ¢ 2 π cosh 2 sinh(φ(x) − jh) cosh 2 sinh(φ(y) − kh) 2. (16) ここでも ψ(·), φ(·) は，式 (6), (7) で定義されたものを利用する．この近似式を用いて次の関数を近似する．例 2. . この表 1 での結果から，関数 f (x) が端点ゼロではない場合には，従来の Sinc 近似では. F (x, y) = exp(x) exp(y),. (x, y) ∈ (0, 1) × (0, 1). 関数 F (x, y) が定義された領域の境界ではゼロではないので，重みなしの従来の Sinc 近. 高精度の実現は出来ないことが明らかである．前節でも述べたように，特に関数近似以外の. 似を用いるのは無意味である．実際の計算では Nx = Ny として，刻み幅は (12) に従い，. 応用では，関数評価回数 N を予め決める必要がある．そのために，重み付 DE-Sinc 近似. α = 1 とした．この関数を格子点 (境界点を除く) si , tl , i, l = 0, 1, · · · , 99 において，. (11) と Gauss-Sinc 近似式 (13) での φ(x)/h の数を調べた結果を表 2 に示す．ここでの関. max |F (si , tl ) − FNx ,Ny (si , tl )|. 数評価回数は (2N + 1) に |φ(x)/h| の整数部分が加わることになる．. 1≤i, l≤99. を重み付 DE-Sinc 近似式 (16) を用いて行った計算結果を次の表 3 に示す．. 3. c 2009 Information Processing Society of Japan °.

(4) Vol.2009-HPC-120 No.5 2009/6/12. 情報処理学会研究報告 IPSJ SIG Technical Report 表 3 例 2. における (16) による最大誤差，刻み幅 DE-Sinc with weight. N 2 4 8 16 32 64. maxerror 7.44e-00 1.23e-00 6.19e-02 3.83e-04 6.86e-08 3.22e-15. において近似解を統一して次のように置くことができ，高精度の数値解が得られる． j=N X cosh(φ(x) − jh) ¡ ¢ (22) uN (x)N = uj S(j, h) ◦ φ(x) cosh2 π2 sinh(φ(x) − jh). h 3.47e-01 3.47e-01 2.60e-01 1.73e-01 1.08e-01 6.50e-02. j=−N. 課題はまだあるものの，このような手法を用いて偏微分方程式の等の高精度の数値解が得られることが期待できる．積分方程式を含め，これらの結果を次の機会で発表する．. 参. 与えることが見て取れる．. 5. 今後の展望このように重みつきの DE-Sinc 法は関数近似において超高精度の数値計算結果を与えるので，これを用いて，積分方程式，常微分方程式，偏微分方程式等の数値解法に適用し，高精度の数値解を得るべく研究をすでに初めている．従来の Sinc 法では，次のような. Z. x. k(x, ξ)u(ξ)dξ = f (x), x, ξ ∈ [a, b]. (17). a. 積分方程式の数値解を得るためには，それの境界条件. Case1 : u(a) = 0, u(b) = 0. Case2 : u(a) 6= 0, u(b) = 0. Case3 : u(a) = 0, u(b) 6= 0. Case4 : u(a) 6= 0, u(b) 6= 0. に応じて，近似解を. Case1 : uN (x) =. N X. uj S(j, h) ◦ φ(x). 文. 献. 1) A.Nurmuhammad, M.Muhammad, M.Mori and M.Sugihara, Double exponential transformation in the Sinc-collocation method for a boundary value problem with fourth-order ordinary differential equation, to appear in J. Comput. Appl. Math. 2) J.Lund, K. L.Bovers, Sinc Methods for Quadrature and Differential Equations, Society for Industrial and Applied Mathematics, Philadelphia 1992. 3) M．Muhammad，M．Mori，Double exponential formulas for numerical indefinite integration，J. Comput．Appl．Math．161 (2003) 431-448. 4) M．Sugihara，Near optimality of the sinc approximation，Math．Comput．72 (2002) 767-786． 5) M.Sugihara, Double exponential transformation in the Sinc-collocation method for two-point boundary value problems, J. Comput. Appl. Math. 149 (2002) 239–250. 6) F. Stenger, Numerical Methods Based on Sinc and Analytic Functions, SpringerVerlag, Berlin, New York, 1993. 7) K. Tanaka, M. Sugihara and K. Murota, Complex-Analytic Approach to the SincGauss Sampling Fourmula, Japan Journal of Industrial and Applied Mathematics, Vol. 25, No. 2, 209–231. 8) H．Takahasi，M．Mori，Double exponential formulas for numerical integration， Publ，RIMS，Kyoto Univ. 9 (1974) 721-741．. この表 3 からも重みつき DE-Sinc 近似式 (16) は 2 変数関数においても高精度の近似を. u(x) − λ. 考. (18). j=−N +1. Case2 : uN (x) = ua wa (x) +. N X. uj S(j, h) ◦ φ(x). (19). j=−N +1. X. N −1. Case3 : uN (x) =. uj S(j, h) ◦ φ(x) + ub wb (x). (20). j=−N. X. N −1. Case4 : uN (x) = ua wa (x) +. uj S(j, h) ◦ φ(x) + ub wb (x). (21). j=−N +1. と置く必要があった．それにたいして，重みつき DE-Sinc 近似法を用いると，全ての Case. 4. c 2009 Information Processing Society of Japan °.

(5)