Euler-Maclaurin の総和公式を利用した数値積分

[研究論文]

Euler-Maclaurin の総和公式を利用した

数値積分

平山弘・加藤俊二

自動車システム開発工学科

Numerical Integration by Euler-Maclaurin Formula

Hiroshi Hirayama and Shunji Katoh

Abstract

Euler-Maclaurin summation formula can be regarded as error evaluation with the value of integration and the value of the numerical integration by the trapezoidal rule. Since the differential coefficients of the function were contained in the error evaluation formula, numerical integration was not performed until now using this.

If the Taylor series method which is a kind of automatic differentiation is used in order to solve this problem, the differential coefficients are calculable with sufficient accuracy. If this method is used, it can be expected that an effective numerical-integration formula can be given.

In this paper, it is shown that Euler-Maclaurin summation formula with the Taylor series method which give the accuracy differential coefficients becomes an effective numerical integration method as same as other leading numerical integration and equivalent grades.

Since the Taylor series method can calculate the value of a function in the singular point of a function with the singularity on appearance with sufficient accuracy, it has the feature which can calculate the numerical integration to a function with the apparent singular with sufficient accuracy.

Keywords: Euler-Maclaurin summation formula, Taylor series, C++ program

1 はじめに Euler-Maclaurin の総和公式は、周期関数の一周 期に渡る積分が高精度で計算できることや、積分 の積分区間の両端で関数値や微分係数が 0 になる 場合、台形公式を使うと高精度で計算できること を説明するためによく使われる公式である。 この公式には、被積分関数の高階微分係数を含む ため、実際の数値積分に使われることは今までほ とんどなかった。高階微分係数を差分法で計算す ると、桁落ちが生じ、精度の良く計算ができない。 このため、これらの公式は積分計算に使われるこ とがほとんどなかった。 関数の微分係数を計算するために差分を使わない 方法として自動微分法 12)_{が知られている。この方} 法と原理的には同じ方法である Taylor 展開法もあ る。これを使えば、高精度で微分係数が計算でき

るので、Euler-Maclaurin の総和公式を使えば、 数値積分を容易に計算できる。 逆に、台形公式の計算部分を級数と見なし、積分 部分が解析的にできるか何等かの方法で簡単に計 算できる場合積分の値を微分を含む部分で補正す ることによって、収束の遅い無限級数が計算でき る。このような研究は長田 11)_{等によって行われて} いる。この場合、微分係数の計算に、数式処理シ ステムが使われている。微分係数の計算に数式処 理システムを使う場合、途中に人間が介入するこ とになるので、少し大きな問題になると、かなり 大きな作業になる。また途中で人的な誤りが入る 可能性が生じる。計算の状況に応じて、微分係数 等を計算するような場合は、実際上不可能になる。 2 Euler-Maclaurin の総和公式 数値積分や級数和の計算でよく使われる公式に Euler-Maclaurin の総和公式がある。この公式は いろいろな文献 2)4)13-14)_{等で紹介されているが、詳} しい説明がない。ここでは主に長田 10)_{を参考にし、} その導出を行った。 関数

f

(

x

)

は区間

[

a

,

b

]

で連続で微分可能とする。 積分区間

[

a

,

b

]

を

n

等分して、台形公式を適用した 計算値と厳密な積分値との差を微分係数を含む

m

項で近似すると、次の式が成り立つ。 m k m k k k k b a n k

R

a

f

b

f

h

k

B

dx

x

f

b

f

kh

a

f

a

f

h

1 ) 1 2 ( ) 1 2 ( 2 2 1 1

)

(

)

(

)!

2

(

)

(

)

(

2

1

)

(

)

(

2

1

ここで、

B

_nは Bernoulli 数(Bernoulli Number)、

n

a

b

h

であり、 b a m m m m

_m

f

t

dt

B

h

R

(

)

)!

2

2

(

) 2 2 ( 2 2 2 2 である。

2.1 Bernoulli の多項式と Bernoulli 数

Bernoulli の多項式 Bn(t)を次のように定義する。 0

!

)

(

1

_k k k x tx

k

x

t

B

e

xe

₍₁₎ (1)の両辺に

e

x

1

の Taylor 展開式を掛け、展開 すると 1 1 0 1 1 0 1

!

)

(

!

)

(

!

)!

1

(

k k j k j k k k k k k k k

k

x

t

B

j

k

k

x

t

B

k

x

k

x

t

(2) ここで、

)!

(

!

!

j

k

j

k

j

k

である。(2)の式の

x

k の係数を等しいと置くと次の式が得られる。 1 0 1

)

(

k j j k

_B

_t

j

k

kt

(3) (3)に

k

1

を代入すると、

(

)

1

0

t

B

、

k

2

を代 入すると次のようになる。

)

(

2

)

(

2

1 0

t

B

t

B

t

この式から、

2

1

)

(

1

t

t

B

となる。この操作を繰 り返すと、次の式が得られる。

6

1

)

(

2 2

t

t

t

B

、

B

t

t

t

t

2

1

2

3

)

(

3 2 3 , … Bernoulli の多項式の定数部分を

B

_nと記述し、 Bernoulli 数と呼ぶ。すなわち、

B

_n

B

_n

(

0

)

であ る。

t

0

を(3)に代入すると、次の Bernoulli 数 の漸化式が得られる

1

0

B

、

0

1 0 k j j

B

j

k

この式を使って、浮動小数点数で Bernoulli 数を 計算すると、桁落ちが生じ、精度良い計算 ができない悪条件の計算であることが知られてい る。 (1) に

t

0

を代入すると Bernoulli 数

B

_n を次の ように定義することもできる。 0

!

1

_k k k x

B

x

_k

e

x

₍₄₎ 上の式に

2

x

_を加え、

_x

_に

_x

_{を代入すると}

2

1

)

1

(

2

2

)

(

1

)

(

) (

x

e

x

e

x

xe

x

e

x

x x x x となり、上の式は偶関数であることがわかる。し たがって、

x

の奇数次の係数はゼロとなり、

0 2 2

)!

2

(

2

1

_k k k x

_k

x

B

x

e

x

となる。このことから、

n

が 3 以上の奇数のとき、

0

n

B

である。 (1)の式に

t

1

を代入すると 0

!

)

1

(

1

_k k k x x

k

x

B

e

xe

₍₅₎ この式は、(4)から、次のように変形できる。 0 ) (

!

)

1

(

1

)

(

1

_k k k k x x x

k

x

B

e

x

e

xe

₍₆₎ (5) と(6)の

x

の係数を比較すると、次の関係式が 得られる。

B

_k

(

1

)

(

1

)

k

B

_k (7) k

B

は

k

が奇数なら、

k

1

以外ではゼロであり、 偶数なら

B

_k

(

1

)

と等しくなる。すなわち

)

,

1

,

0

(

)

0

(

)

1

(

₂ 2

B

B

k

B

k k k (8)

)

1

(

0

)

0

(

)

1

(

_{2 1} 1 2

B

k

B

k k となる。 (1)の式の両辺を

t

で微分すると 0 2

!

)

(

1

_k k k x tx

k

x

t

B

e

e

x

₍₈₎ (1)と(8)の係数を比較することによって、次の関 係式が得られる。

(

)

(

)

1

t

nB

t

B

n n (10) 微 分 の 公 式 (10) を 使 う と 、 Bernoulli の 多 項 式

)

(

2

x

B

_n を区間

[

0

,

1

]

で、次の様に Fourier 級数に展 開できる。Fourier 級数の

cos(

2

k

x

)

の係数

c

_kは, 次の式で表される。 1 0 2

2

cos

)

(

2

B

x

k

xdx

c

k n これから、次の式が得られる。 1 2 2 1 2

2

cos

)

2

(

)!

2

(

2

)

1

(

)

(

k n n n n

x

n

_k

k

x

B

(11)

)

(

1 2

x

B

_n のように奇数次数の Bernoulli の多項式 は、このようには Fourier 級数に展開できない。 この Fourier 級数に

x

0

を代入すると 1 2 2 1 2

1

)

2

(

)!

2

(

2

)

1

(

k n n n n

n

_k

B

が得られる。

n

1

のとき、級数部分は 2 より小さ いから、次の関係式が得られる。

B

_n

n

_n 2 2

)

2

(

)!

2

(

4

(12) (11)の Fourier 級数から、区間

[

0

,

1

]

では、

x

0

の時が

B

_n

(

x

)

の絶対値が最大値になるこ とがわかる。すなわち

B

_n

x

B

_n 2

)

(

(

0

x

1

)

(13)

2.2 Euler-Maclaurin

の総和公式

次の積分

I

_jk , を定義する。 _jk h _k

f

k

x

_j

t

dt

h

t

B

k

I

0 ) 2 ( 2 ,

(

)

)!

2

(

1

こ こ で 、

x

_j

hj

、

n

a

b

h

で あ る。この式で

1

k

の場合を考える。この式は部分積分を使って 次のように変換できる。 1

)

(

1

)]

(

)

(

[

2

1

)

(

)

(

!

2

)

(

!

2

1

)

(

!

2

1

2 1 1 2 0 2 2 2 0 2 1 , j j x x j j j j h j h j j

dt

t

f

h

x

f

x

f

h

x

f

x

f

B

dt

t

x

f

B

h

t

h

t

dt

t

x

f

h

t

B

I

(14) 上の式(14)を

j

0

,

,

n

1

について加えると、 次の式が得られる。 b a n k n j j

dt

t

f

h

b

f

kh

a

f

a

f

h

a

f

b

f

B

I

)

(

1

)

(

2

1

)

(

)

(

2

1

1

)]

(

)

(

[

!

2

2 1 1 2 1 0 1 , (15) (8) と(10) の関係式を利用すると h j k k k j

_k

B

_h

t

f

x

t

dt

I

0 ) 2 ( 2 ,

(

)

)!

2

(

1

1 , 2 ) 1 2 ( 1 ) 1 2 ( 2 ) 1 2 ( 0 2 1 0 ) 1 2 ( 2

1

)

(

)

(

)!

2

(

)

(

)!

1

2

(

1

)

(

)!

2

(

1

k j j k j k k j k h k h j k k

I

h

x

f

x

f

k

B

dt

t

x

f

h

t

B

h

k

t

x

f

h

t

B

k

したがって、次のようになる。

)

(

)

(

)!

2

(

) 1 2 ( 1 ) 1 2 ( 2 2 , 2 1 , j k j k k k j k j

x

f

x

f

k

h

B

I

h

I

(16) (15) の 1 , j

I

に(16)の関係式を代入すると m k m k k k k b a n k

R

a

f

b

f

h

k

B

dx

x

f

b

f

kh

a

f

a

f

h

1 ) 1 2 ( ) 1 2 ( 2 2 1 1

)

(

)

(

)!

2

(

)

(

)

(

2

1

)

(

)

(

2

1

を得る。ここで、 1 0 ) 2 2 ( 0 2 2 2 2 1 0 1 , 2 2

)

(

)!

2

2

(

n j j m h m m n j jm m m

dt

t

x

f

h

t

B

m

h

I

h

R

である。(13)を使い上の式の

R

_mを評価すると以下 のようになる。 b a m m m m

_m

f

t

dt

B

h

R

(

)

)!

2

2

(

) 2 2 ( 2 2 2 2 (17) さらに(12) 式を適用すると、 b a m m m

h

f

t

dt

R

(

)

2

4

(2 2) 2 2 (18) ここで、等式が成り立つのは、右辺の積分が0 の ときである。多くのEuler-Maclaurin の総和公式 の 誤 差 評 価 は 、 (17) の 式 で 与 え て い る が 、 Bernulli 数の大きさがわからないため直感的にそ の大きさがわかりにくい。(18) の式は少し過大評 価であるが、直感的に分かりやすくなっている。 この式からたとえ周期関数の一周期に渡る積分で も、微分係数が非常に大きくなる関数は、効率的 に求められない事になる。たとえば 2 0

)

sin

50

cos(

x

dx

I

(19) の場合、1 回微分する毎にだいたい50 倍位の数値 になる。このような場合、(18) の誤差評価からあ まり急速に収束するようにはならないと推定でき る。分割数90でほぼ15 桁の精度が得られる。もし、 積分が 2 0

)

cos(sin

x

dx

I

(20) だったら、分割数16 で15 桁の精度が得られる。 台形公式を使った場合、分割数を2, 4, 8,16 と増 やして計算する場合が多いので分割数が偶数にな ることが多い。このためあまり気づかれていない が、分割数を奇数にすると非常に高精度の結果が 得られる。(20) の問題の場合、分割数が9でも16 の場合よりも高精度の結果が得られる。分割数が7 の場合でも16 の場合とほぼ同精度の結果が得られ る。(19) の問題の場合でも同様な結果が得られる。 分割数が90 の結果より45 の場合が良い結果が得 られる。 3 Euler-Maclaurin の総和公式による数値積分 Euler-Maclaurin の総和公式を使って数値積分 を行うには、(1)の式からわかるように、最初に、 積分区間の端点において、Taylor 展開を計算する。 このTaylor 展開を計算するプログラムは、通常の 関数値の計算と宣言部分を除いてほぼ同じになる。 このTaylor 展開を行うためのテンプレート・プロ グラムも公開されている6)_{のその計算は容易である。} テンプレート機能3)_{を使えばTaylor 展開と関数値} を計算するプログラムは多くの場合１個のプログ ラムで記述できる。実際この論文で示した計算例 は、C++言語のテンプレート機能を使って１個の関 数 プ ロ グ ラ ム の 形 で 記 述 し た 。 た と え ば 、

x

x cos

cosh

92

.

0

_{をプログラムで記述するには、} 次の様に書く。 template<typename T> T func( const T& x )

{ T s ; s =0.92*cosh(x)-cos(x) ; return s ; } このように1 個関数プログラムを書けば、関数計 算やTaylor展開の計算に使える。次に分割数

n

を 決める。分割数がある程度推定できる場合、その 値にする。推定が出来ない場合、最小値の2 にす る。この分割数を利用して、台形公式を利用して、 積分の近似値を計算する。近似値を次の項を利用 して、補正する。この級数は漸近級数なので漸近 級数の手法を使って計算する。

)

(

)

(

)!

2

(

)

(

2

h

2

f

(2 1)

b

f

(2 1)

a

k

B

k

c

k k k k この補正項

c

(

k

)

の絶対値が要求精度以下になるま で、補正項を加算して行く。補正項の絶対値が要 求精度より小さくなったら、その項を使って補正 し計算を止める。この計算値が求める積分値にな る。補正項の絶対値が要求精度より小さくならな い場合や逆に大きくなった場合、補正するのを止 め、分割数

n

を 2 倍にして、台形公式による積分 の計算に戻る。分割数を 2 倍にすると、その以前 の台形公式計算値を利用し、計算量を少なく出来 るためである。この手順を繰り返す事によって積 分値を計算する。

3.1

簡単な計算例

簡単な例として、次の積分を計算する。これを 要求精度 9

10

で計算する。 1 0

₁

1

dx

x

I

最初に、積分区間の両端点で Taylor 展開をする。 展開は 20 次まで行った。その計算結果を以下に示 す。x = 0 における Taylor 展開を 15 次まで表示 すると次の様になる。 15 14 13 12 11 10 9 8 7 6 5 4 3 2

1

x

x

x

x

x

x

x

x

x

x

x

x

x

x

x

1

x

における Taylor 展開を 5 次まで表示すると 次の様になる。 5 4 3 2

1)

-0.015625(

-1)

-0.03125(

1)

-0.0625(

-1)

-0.125(

1)

-0.25(

-0.5

x

x

x

x

x

2

n

と し て 、 台 形 公 式 で 数 値 積 分 す る と 、 0.708333333333333 となる。これを補正する。 この場合の補正項は収束が非常に遅く要求精度 を満たすことができない。

k

8

で逆に絶対値が増 加し始める。このため分割数を 2 倍の

n

4

として、 台形公式を使って数値積分を行う。このときの積 分値は、0.697023809523809 となる。この結果を 再度補正する。ここでのは補正項も十分速く小さ く な る と は 言 え な い が 、

k

7

で 補 正 項 が 9

10

1

.

3

となり要求精度より小さくなった。そこ ま で の 計 算 値 が 積 分 値 と な る 。 積 分 値 は 0.6931471804863029718 となった。要求精度通りの 結果が得られた。このときの関数計算回数は、5 回でその内 2 回が Taylor 展開を計算するための 関数計算であった。 分割数

n

を増加させると、

h

は小さくなるので 補正項

c

(

k

)

はすばやく小さくなる。このため、高 次の微分係数はあまり必要としなくなる。Taylor 展開は低次数でも Euler-Maclaurin の総和公式を 利用できることになる。逆に Taylor 展開が高次数 であるならば、分割数

n

を小さく出来る。上の計 算では 20 次の Taylor 展開を利用したが、それが 適切かどうか考慮する必要がある。

3.2

数値例

Euler-Maclaurin の総和公式を利用した計算法 の性能を評価するために、Kahaner の問題8)_の中か ら両端点で特異性を持つため Taylor 展開できない 問題等を除いた 13 問題について計算を行った。計 算した問題を表 1 に示す。番号は Kahaner の問題 の番号である。それを計算した結果を示す。 実行結果を表 2 に示す。N は分割数すなわち関 数の計算回数、Error はその積分ルーチンが出力 した誤差である。分割数の内 2 回は Taylor 展開の 計算である。EM は Euler-Maclaurin の総和公式を 利用した計算、DE は二重指数型積分公 15）_(Double Exponential formula)、DAQN9 は、適応型ニュー トン・コーツ法 8) _{の結果はその参考文献による結}

果である。IBM 型 64 ビットの浮動小数点演算によ る結果である。DE のプログラムとして、ネットで 公開されている大浦9)_{の C 言語用 DE プログラムを} C++言語に変換し使用した。 コンパイラーとして Microsoft Visual C++ 2012 を使用した。 4 まとめ Euler-Maclaurin の総和公式を利用して、数値 積分する方法を提案した。この方法は、単純でわ かり易い計算法であるにも関わらず、多くの場合、 有力な数値積分法と同等程度またはそれ以上の性 能を発揮する。 この計算を行うには、被積分関数の微分係数を 精度良く計算する必要がある。この方法として、 自動微分法があるが、一般のユーザーにそれほど 普及していない。この方法が使われるためにも自 動微分の普及することを期待したい。 5 参考文献

[1]Abramowitz M. and Stegun I. A.: Handbook of Mathmatical Functions,Dover,(1972)

[2] Davis P.J., Rabinwitz P.(森 正武訳): 計算 機による数値積分法, 日本コンピュータ協会, (1981)

[3] David Abrahams, Aleksey Gurtovoy: C++ Template Metaprogramming, Addison Wesley, (2005)

[4] Gisela Engeln-M ¨ ullges, Frank Uhlig: Numerical Algorithms with Fortran, Springer, 1996

[5] Henrici P.: Applied and Computational Complex Analysis, Vol. 1, John Wiley & Sons, New York, (1974)

[6] 平山 弘、館野 裕文、浅野 直之、川口 隆史: Taylor 級数演算ライブラリの使用法, 東北大学 情報シナジー セ ン タ ー 大 規 模 科 学 計 算 シ ス テ ム 広 報 SENAC, 40(2007) 29 68 [7] 平山 弘、小宮 聖司、佐藤 創太郞:Taylor 級 数法による常微分方程式の解法, 日本応用数理 学会論文誌、12(2002), pp.1 8 [8] 二宮 市三: 適応型ニュートン・コーツ積分法 の改良, 情報処理学会誌,21(1980),504 513 [9] 大浦 拓哉: Ooura’s Mathematical Software

Packages, ttp://www.kurims.kyotou.ac.jp/ooura/index-j.html [10] 長田 直樹: 対数収束級数の漸近展開と加速 法, 情報処理学会論文誌 29(1988), 256 261 [11] 長田 直樹:お話：数値解析第 3 回オイラー・ マ ク ロ ー リ ン の 公 式 , http://www.cis.twcu.ac.jp/osada/rikei/rikei 2008-7.pdf

[12] Rall,L. B.: Automatic Differentiation Technique and Applications, Lecture Notes in Computer Science, Vol. 120, Springer Verlag, Berlin-Heidelberg-New York, (1981) [13] 篠原 能材: 数値解析の基礎, 7 章, 日新出

版, (1978)

[14] 杉原 正顕, 室田 一雄: 数値計算法の数理, 12 章, 岩波書店, (2007)

[15] Takahasi, H. and Mori, M.: Double exponential formula for numerical integration, Publ. RIMS, Kyoto Univ., 9(1974), 121–141

Table 1 Test Problem

No.

Integral

1 1

_1.71828182

₈₅

0

e

dx

x 4 1

₀

_.

₉₂

_cosh

_cos

_0.47942822

₆₆₉

1

x

xdx

5

37

1.58223296

9

.

0

1

1 1

x

4

x

2

dx

8

734

0.86697298

1

1

1 0

x

4

dx

9

90

1.15470066

)

4159

.

31

sin(

2

2

1 0

x

dx

10

056

0.69314718

1

1

1 0

x

dx

11

304

0.37988549

1

1

1 0

e

x

dx

12 1 0

411

0.77750463

1

dx

e

x

x 13 10 0

52566

0.00909864

14159

.

3

)

159

.

314

sin(

dx

x

x

16 10 0 2

287

0.49936380

)

1

2500

(

14159

.

3

50

dx

x

17 1 01 . 0 2

963

0.11213956

14159

.

3

50

)

50

sin(

50

dx

x

x

18 0

338

0.83867632

)

3

cos

3

2

sin

3

2

cos

2

sin

3

cos(cos

x

x

x

x

x

dx

20

41

1.56439644

005

.

1

1

1 1

x

2

dx

Table 2 Comparison of Performance of Quadrature Routines(Tolerance =

9

10

)

EM DE DAQN9

No. N Error N Error N Error

１ ３ 3.5e-11 33 3.4e-91 25 5.2e-18

４ ３ 2.5e-10 67 2.2e-9 25 3.5e-17

５ ９ 2.2e-10 65 6.3e-9 61 1.8e-11

8 9 1.7e-10 65 3.5e-9 91 8.9e-12

9 33 5.4e-10 519 3.7e-9 81 1.0e-5

10 5 3.1e-10 33 1.4e-9 21 3.9e-13

11 3 2.6e-10 33 7.6e-11 21 3.3e-18

12 3 7.3e-11 33 1.6e-9 21 2.2e-18

13 129 7.9e-10 461 1.2e-8 321 6.2e-7

16 5 3.2e-10 141 4.2e-9 91 5.2e-7

17 129 1.9e-10 555 5.7e-9 101 7.4e-4

18 33 1.1e-10 131 1.7e-8 51 5.4e-6