積分区間分割による数値積分の効率的な計算法
平山弘・加藤俊二
自動車システム開発工学科
Efficient Computation Method of Numerical Integration
by Integral Interval Division
Hiroshi HIRAYAMA, Shunji KATOH Abstract
When an integral interval [a; b] of the integral with the bad condition is divided into two or more intervals and higher order numerical integral formula is applied every area and the whole one, good results are sometimes obtained with few sample points.
In this paper, we show that you can get the condition to define convergence speed of numerical
integration and partition method from numerical integration method of contour integral on complex plane. Keywords: Numerical Integration, Division of Integral Interval
はじめに 次のような有限区間[𝑎𝑎, 𝑏𝑏]に渡る関数𝑓𝑓(𝑥𝑥)の積分を考え る。 𝐈𝐈 = ∫ 𝒇𝒇(𝒙𝒙)𝒅𝒅𝒙𝒙𝒃𝒃 𝒂𝒂 ここでは、関数𝑓𝑓(𝑥𝑥)は積分区間の近くに特異点を持つよ うな条件の悪い積分とする。 式 で表される積 分は、&DXFK\の積分表示を使うと、周期関数の１周期に 渡る積分に変換できる。周期関数の１周期に渡る積分 は、台形公式を使って効率的に計算できることが知られ ているので、この変換によって効率的に数値積分の計算 できると期待できる。積分区間内で特異性を持つ場合、 被積分関数は、𝑓𝑓(𝑥𝑥) = 𝑤𝑤(𝑥𝑥)𝑔𝑔(𝑥𝑥)のように表現することが できる。 ここで、𝑤𝑤(𝑥𝑥)は特異性を持つ部分、𝑔𝑔(𝑥𝑥)は積分区間に 内で正則な部分である。𝑤𝑤(𝑥𝑥)と𝑔𝑔(𝑥𝑥)の選び方は一意性は ないが、以下の議論ではどのようにとっても同じ結果と なる。 &DXFK\の積分表示を利用すると、次のように周回積分 に変換することができる。 𝑰𝑰 = ∫ 𝒇𝒇(𝒙𝒙)𝒅𝒅𝒙𝒙 = ∫ 𝒘𝒘(𝒙𝒙)𝒃𝒃 𝟐𝟐𝝅𝝅𝝅𝝅𝟏𝟏 𝒂𝒂 ∫ 𝒈𝒈(𝒈𝒈) 𝒈𝒈 𝒛 𝒙𝒙 𝒄𝒄 𝒅𝒅𝒈𝒈𝒅𝒅𝒙𝒙 𝒃𝒃 𝒂𝒂 = 𝟐𝟐𝝅𝝅𝝅𝝅 ∫ 𝒈𝒈𝟏𝟏 (𝒈𝒈) [∫𝒃𝒃𝒘𝒘(𝒙𝒙)𝒈𝒈 𝒛 𝒙𝒙 𝒅𝒅𝒙𝒙 𝒂𝒂 ] 𝒅𝒅𝒈𝒈 𝒄𝒄 = 𝟐𝟐𝝅𝝅𝝅𝝅 ∫ 𝒈𝒈𝟏𝟏 (𝒈𝒈)𝑯𝑯(𝒈𝒈)𝒅𝒅𝒈𝒈 𝒄𝒄 積分路cは、図 のように積分区間を囲み、積分区間以 外に特異点を含まな いような閉積分路で ある。
図積分路 [研究論文]括弧>@の部分𝐻𝐻(𝑧𝑧)は、関数𝑤𝑤(𝑥𝑥)の +LOEHUW変換と呼ば れる。たとえば、𝑤𝑤(𝑥𝑥) = 1の場合 ∫𝑥𝑥𝑥𝑥𝑤𝑤(𝑤𝑤)𝑧𝑧𝑧𝑤𝑤𝑑𝑑𝑥𝑥 = log𝑤𝑤𝑧𝑥𝑥𝑤𝑤𝑧𝑥𝑥であるから、次のようになる。 ∫ 𝑓𝑓(𝑥𝑥)𝑥𝑥 𝑥𝑥 𝑑𝑑𝑥𝑥 = 1 2𝜋𝜋𝜋𝜋 ∫ 𝑓𝑓(𝑥𝑥)log 𝑧𝑧 − 𝑎𝑎 𝑧𝑧 − 𝑏𝑏 𝑑𝑑𝑧𝑧 𝑐𝑐 𝑤𝑤(𝑥𝑥) = 1のとき以外の有限区間の +LOEHUW 変換例を表１ に示す 議論を単純化するために、次のような変数変換を行う。 𝑧𝑧 = 𝑤𝑤 +𝑥𝑥+𝑥𝑥2 と置くと、積分は𝑤𝑤の積分になり積分区間は [−𝑒𝑒, 𝑒𝑒] 𝑒𝑒 𝑥𝑥𝑧𝑥𝑥2 となる。さらに、ジューコフスキー -RXNRZVNL 変換 𝑤𝑤 =12 (𝑢𝑢 +𝑒𝑒𝑢𝑢 )2 を行うと、元の積分区間は、原点を中心とする半径eの 円になる。このような変換を行えば、特異点がどこにあ っても、常に積分路を原点を中心とする円にとることが できる 。 本論文では、このように変形された積分は、高精度で計 算でき、どの程度の速さで収束するか計算できる。収束 の速さを高める方法として、積分区間を分割する方法を 提案する。この手法を使えば、これまで多くの自動積分 法で２～ 桁程度しか計算出来なかった .DQHU の 番目 の問題も問題なく ２桁以上の精度で計算可能である。 積分の収束率 式 のように変形され、平行移動や -RXNRZVNL 変換が なされた被積分関数をF(𝑧𝑧)とする。この関数は、次のよ うにローラン /DXUHQW 展開できるものとする F(𝑧𝑧) = ∑ 𝑎𝑎𝑘𝑘𝑧𝑧𝑘𝑘 ∞ 𝑘𝑘=0 + ∑ 𝑏𝑏𝑘𝑘𝑧𝑧𝑧𝑘𝑘 ∞ 𝑘𝑘=1 次のような積分を計算することを考える。 I =2𝜋𝜋𝜋𝜋 ∫ 𝐹𝐹1 (𝑧𝑧)𝑑𝑑𝑧𝑧 𝑐𝑐 ＝ 1 2𝜋𝜋𝜋𝜋 ∫ ∑ 𝑎𝑎𝑘𝑘𝑧𝑧𝑘𝑘 ∞ 𝑘𝑘=0 𝑑𝑑𝑧𝑧 𝑐𝑐 + 1 2𝜋𝜋𝜋𝜋 ∫ ∑ 𝑏𝑏𝑘𝑘𝑧𝑧𝑧𝑘𝑘 ∞ 𝑘𝑘=1 𝑑𝑑𝑧𝑧 𝑐𝑐 の第 項の積分𝐽𝐽は、原点を中心とした半径𝑟𝑟の円を 積分路𝑐𝑐にとり、𝑁𝑁等分に分割し台形公式を使って計算 するとその値を𝐽𝐽𝑁𝑁とすると 𝐽𝐽𝑁𝑁=𝑁𝑁 ∑ 𝑎𝑎1 𝑘𝑘𝑟𝑟𝑘𝑘+1 ∞ 𝑘𝑘=0 (∑ 𝑒𝑒𝑖𝑖(𝑘𝑘+1)ℎ𝑗𝑗 𝑁𝑁𝑧1 𝑗𝑗=0 ) ここで、h =2𝜋𝜋 𝑁𝑁である。この式の括弧 の部分Sは、等 比数列なので、容易に計算出来る。𝑘𝑘 + 1が𝑁𝑁の倍数にな ると、すべての項が となり和は𝑁𝑁となる。𝑘𝑘 + 1が𝑁𝑁の 倍数でない場合、等比数列として和を計算すると分子が 零になることから となる。すなわち、次のようにな る。 𝑆𝑆 =1 − 𝑒𝑒1 − 𝑒𝑒𝑖𝑖(𝑘𝑘+1)ℎ𝑁𝑁𝑖𝑖(𝑘𝑘+1)ℎ = {𝑁𝑁 mod0 mod(𝑘𝑘 + 1, 𝑁𝑁) ≠ 0(𝑘𝑘 + 1, 𝑁𝑁) = 0 したがって、次の式が得られる。 𝐽𝐽𝑁𝑁= ∑ 𝑎𝑎𝑘𝑘𝑁𝑁𝑧1𝑒𝑒𝑘𝑘𝑁𝑁 ∞ 𝑘𝑘=1 ここで、 𝑎𝑎𝑚𝑚=2𝜋𝜋𝜋𝜋 ∫1 𝐹𝐹(𝑧𝑧)𝑧𝑧𝑚𝑚 𝑐𝑐1 𝑑𝑑𝑧𝑧 であるから、 𝐽𝐽𝑁𝑁=2𝜋𝜋𝜋𝜋 ∫ 𝐹𝐹(𝑧𝑧) ∑ (1 𝑟𝑟𝑧𝑧) 𝑘𝑘𝑁𝑁 ∞ 𝑘𝑘=1 𝑑𝑑𝑧𝑧 𝑐𝑐1 =2𝜋𝜋𝜋𝜋 ∫1 𝐹𝐹(𝑧𝑧) (𝑟𝑟𝑧𝑧) 𝑁𝑁 1 − (𝑟𝑟𝑧𝑧)𝑁𝑁 𝑑𝑑𝑧𝑧 𝑐𝑐1 この積分値は正則関数𝐹𝐹(𝑧𝑧)の積分だから、コーシー &DXFK\ の積分定理により零である。したがって、 の𝐽𝐽𝑁𝑁が の 項から生じる誤差となる。 の積分を 評価するために、積分路𝑐𝑐1として原点を中心とする半径 𝑅𝑅1の円周上の積分を考える。半径𝑅𝑅1と半径𝑟𝑟は自由に選 べるから、r < 𝑅𝑅1とし、半径𝑟𝑟の円周上での𝐹𝐹(𝑧𝑧)の絶対値 の最大値を𝑀𝑀1とすると次のように評価できる。 |𝐽𝐽𝑁𝑁| ≤ ( 𝑟𝑟𝑅𝑅 1) 𝑁𝑁 1 − ( 𝑟𝑟𝑅𝑅 1) 𝑁𝑁𝑀𝑀1< 𝑀𝑀1(𝑅𝑅𝑟𝑟 1) 𝑁𝑁 𝑅𝑅1が大きければ、この誤差は、小さく出来る。𝑅𝑅1は、 原点に最も近い特異点までの距離にとると最大のにな る。 同様に式 の第 項の積分から𝑘𝑘 = 1の項を除いた積分 𝐾𝐾𝑁𝑁＝ 12𝜋𝜋𝜋𝜋 ∫ ∑ 𝑏𝑏𝑘𝑘𝑧𝑧𝑧𝑘𝑘 ∞ 𝑘𝑘=1 𝑑𝑑𝑧𝑧 𝑐𝑐 も評価できる。𝑘𝑘 = 1の項の積分値は、𝑏𝑏1となるので、 その項を除いて評価すると誤差評価になる。半径𝑅𝑅2上で の関数𝐹𝐹(𝑧𝑧)の最大値を𝑀𝑀2とすると |𝐾𝐾𝑁𝑁| < 𝑀𝑀2(𝑅𝑅𝑟𝑟 )2 𝑁𝑁 𝑅𝑅2が小さくすれば、この誤差も、小さく出来る。𝑅𝑅2 は、原点からの積分区間までの距離𝑒𝑒になった時最小のに なる。 したがって、誤差𝐸𝐸𝑟𝑟は、次のようになる。 𝐸𝐸𝑟𝑟 = |𝐽𝐽𝑁𝑁| + |𝐾𝐾𝑁𝑁| = 𝑀𝑀1(𝑅𝑅𝑟𝑟 1) 𝑁𝑁 + 𝑀𝑀2(𝑅𝑅𝑟𝑟 )2 𝑁𝑁
7DEOH )LQLWH+LOEHUW7UDQVIRUPDWLRQ
:KHUHΓ(𝑥𝑥)LV*DPPDIXQFWLRQ F(𝑎𝑎, 𝑏𝑏; 𝑐𝑐; 𝑥𝑥)LV +\SHUJHRPHWULFIXQFWLRQ誤差𝐸𝐸𝐸𝐸が最小になるのは、相加相乗平均の式を適用す ると、第 項と第 項の式が一致した場合である。これ から、 𝐸𝐸 = √𝑀𝑀𝑀𝑀2 1 2𝑁𝑁 √𝑅𝑅1𝑅𝑅2 のとき、誤差が最小になる。𝑁𝑁が十分大きいとき r = √𝑅𝑅1𝑅𝑅2 となる。このとき、誤差𝐸𝐸𝐸𝐸は次のようになる。 Er = (𝑀𝑀1+ 𝑀𝑀2) (𝑅𝑅𝑅𝑅2 2) 𝑁𝑁 2 通常は -RXNRZVNL 変換で半径 の円に変換するので、 𝑅𝑅2= 1 𝑒𝑒 = 1 となる。𝑅𝑅1はこのように変換したときも原 点から最も近い特異点までの距離となる。したがって、 誤差𝐸𝐸𝐸𝐸は、 Er = (𝑀𝑀1+ 𝑀𝑀2)𝑒𝑒−(12log𝑅𝑅1)𝑁𝑁= 𝑂𝑂(𝑒𝑒−𝑞𝑞𝑁𝑁) となる。この𝑞𝑞 =1 2log𝑅𝑅1を収束率と定義する。収束率𝑞𝑞 が大きいほど、速く収束する。𝑞𝑞は特異点が積分区間に 近いほど小さくなる。積分区間が広くなるほど、相対的 に特異点は原点に近づくため、小さくなり、収束が悪く なる。 式 のように与えられた積分の収束率𝑞𝑞は、次のような 手順で計算出来る。与えられた被積分関数𝑓𝑓(𝑥𝑥)の特異点 を複素数𝑧𝑧とし、積分区間を[𝑎𝑎, 𝑏𝑏]とし、&言語風の関数 で書くと以下のようになる。この関数を𝑞𝑞(𝑧𝑧, 𝑎𝑎, 𝑏𝑏)とすと GRXEOHT FRPSOH[]GRXEOHDGRXEOHE ^ FRPSOH[X ] DE ED FRPSOH[Y VTUW X X GRXEOHW ORJ PD[ DEV XY DEV XY UHWXUQW ` 𝑞𝑞(𝑧𝑧, 𝑎𝑎, 𝑏𝑏)が求める収束率である。 積分区間の分割 積分区間を分割する前の収束率を𝑞𝑞、標本点数を𝑁𝑁とす る。積分区間を２分割したとき、それぞれの区間積分の 収束率を𝑞𝑞1および𝑞𝑞2としそれぞれの区間の標本点数を 𝑁𝑁1、𝑁𝑁2とする。 min (𝑞𝑞1𝑁𝑁1、𝑞𝑞2𝑁𝑁2) > 𝑞𝑞𝑁𝑁 であり、かつ 𝑁𝑁1+ 𝑁𝑁2< 𝑁𝑁 が成り立つならば、分割数が少なくて、精度の良い結果 が得られると期待できる。３以上の分割も同様に、求め られる。 上の式からわかるように、積分公式の誤差が𝑂𝑂(𝑒𝑒−𝑞𝑞𝑁𝑁) である公式が対象となる。すなわち、ガウスの数値積分 やここで扱った周回積分法などが対象となる。この収束 率が分割によって非常に効率的になる場合は、多くの数 値積分公式でも有効である場合多い。 実際この分割を、決定するには、𝑞𝑞1= 𝑞𝑞2と仮定して、 区間を求め、 と の条件を満たすかどうかを判定 する。この計算は、一般に連立非線形方程式となる。
数値例 数値例 次の積分 𝐼𝐼1= ∫ 𝑑𝑑𝑥𝑥 𝑥𝑥2+ 1 100 1 −1 = 20 tan −110 = 29.422 ⋯ を計算する。この被積分関数を図２に示す。𝑥𝑥 = 0で鋭 いピークを持つ関数である。 よく知られているように、このような積分は、>@と >@と分割すると効率的に計算出来る。実際収束率を計 算すると区間>@の積分の収束率𝑞𝑞(0.1𝑖𝑖, 0,1)は、 となり、区間>@の収束率𝑞𝑞(0.1𝑖𝑖, −1,1)は、 と なった。 収束率が約 倍で 倍以上であるから、積分区間を分 割すれば、効率的に計算可能と思われる。これを計算し た結果を表 に示した。 ガウスの数値積分公式を使って計算したものが非分割 計算で、標本点数の半分の標本点数のガウスの数値積分 公式を使って、 分割した部分の積分値を計算し加算し たものが 分割計算である。たとえば、全標本点数 の 場合、区間>@の積分を 点のガウスの数値積分で 計算したものが非分割計算で、 点のガウスの数値積分 法を区間>@と>@に適用して計算したのが 分割 計算である。誤差は相対誤差である。 7DEOH,QWHJUDWLRQRI𝐼𝐼1E\GLYLVLRQ 全本点数 1 ２分割誤差 非分割誤差 H H H H H H H H ガウス型数値積分公式の分点と重みは 'DYLV にしたが って、多倍長数 を利用して計算した。 分割する場合、区間>@を>F@、>FF@、 >F@と 分割し、 つの区間で収束率を同じになるよ うに、F と を選ぶ。この条件で F を選ぶと F となる。対称性から F とな る。 つの区間で集束率は で、区間>@の収 束率の約 倍となり 倍以上であるから、分割して計 算するのは意味があると思われる。実際計算した結果を 表 に示した。 図２ 関数 [A のグラフ7DEOH,QWHJUDWLRQRI𝐼𝐼1E\GLYLVLRQ 全本点数 1 分割誤差 非分割誤差 H H H H H H H H この問題を二重指数型数値積分法 で要求精度 H で計算すると、分割なしの時標本点数 、誤差 H となり、原点で 分割すると、標本点数の合計が 、誤差が H であった。この問題に対しては、 ２分割は二重指数型数値積分法に有効であることがわか る。 上の例題と同じように、積分区間を 分割し、要求精度 H で計算すると、標本点数の合計が で、誤差 が H となった。この問題に対して、３分割は２分 割より効果がないことがわかる。 数値例 次の積分 𝐼𝐼2= ∫ 𝑑𝑑𝑑𝑑 𝑑𝑑 + 1100 1 0 = log 101 = 4.61512051 ⋯ を計算する。このような積分は、積分区間を分割して計 算することは通常ないが、以下のように分割することに よって、効率的に計算できることがわかる。この積分を 分割して積分する。この場合、分割点 F となる。これを計算すると表 のようにな る。 分割された区間の積分の収束率 、全区間 の収束率は となり約 倍で 倍以上な ので、効率的に計算できると推定できる。 7DEOH,QWHJUDWLRQRI𝐼𝐼2E\GLYLVLRQ 全本点数 1 ２分割誤差 非分割誤差 H H H H H H H H この積分を 分割して積分する。この場合、分割点 F 、F となる。これを計 算すると表 のようになる。 分割された区間の積分の収束率 、全区間 の収束率は となり約 倍で 倍以上な ので、効率的に計算できると推定できる。 7DEOH,QWHJUDWLRQRI𝐼𝐼2E\GLYLVLRQ 全本点数 1 分割誤差 非分割誤差 H H H H H H H H この積分は積分区間の端点近くに特異点がある。多くの 数値積分公式は、積分区間の端点における特異性があっ ても効率的に計算できるので、このような分割は不要で ある。 ここで使ったガウス型数値積分公式は、積分区間の端点 における特異性は考慮していない公式なので分割は非常 に有効である
。
この問題を二重指数型数値積分法で計算すると、分割な しの時標本点数 、誤差 H となり、𝑑𝑑 = 0.0904987475で 分割すると、標本点数の合計が 、 誤差が H であった。この問題に対しては、２分割 は二重指数型数値積分法には有効でないことがわかる。 𝑑𝑑 = 0.0𝑥65700877と𝑑𝑑 = 0.206877𝑥911で積分区間を 分割し、要求精度 H で計算すると、標本点数の合 計が で、誤差が H となった。この問題に対し ては、２分割でも 分割しても効果がないことがわか る。 数値例 .DKDQHU の 番目の数値積分問題を考える。この問題 は多くの自動数値積分法で扱われている問題で、これま での自動数値積分法では 桁以上の精度で計算出来なか った問題である。この積分は、 𝐼𝐼3= ∫ 𝑓𝑓0(𝑑𝑑)𝑑𝑑𝑑𝑑 1 0 ここで、𝑓𝑓0(𝑑𝑑) = 𝑓𝑓1(𝑑𝑑) + 𝑓𝑓2(𝑑𝑑) + 𝑓𝑓3(𝑑𝑑) 𝑓𝑓1(𝑑𝑑) =(cosh 10(𝑑𝑑 − 0.2))1 2 𝑓𝑓2(𝑑𝑑) =(cosh 100(𝑑𝑑 − 0.4))1 4 𝑓𝑓3(𝑑𝑑) =(cosh 1000(𝑑𝑑 − 0.6))1 6 である。この積分は解析的に計算出来て 𝐼𝐼3= 0.2108027𝑥5500549277𝑥7564𝑥 ⋯ が得られる。これを関数のグラフを図 に示す。𝑑𝑑 = 0.4と𝑑𝑑 = 0.6に鋭いピークを持つことから、この近辺に特 異点を持つと推定できる。 𝑑𝑑 = 0.4付近では、z = 0.4 ±100𝜋𝜋 𝑖𝑖、𝑑𝑑 = 0.6付近では z = 0.6 ±1000𝜋𝜋 𝑖𝑖の特異点を持っている。特に𝑑𝑑 = 0.6近くの特 異点の位置から収束率を計算すると が得ら れる。数値例１では収束率が約 でも分割が有効だっ たので、この収束率では当然分割すべきである。したが って積分は、区間>@と>@の少なくとも二つに 分割すべきである。 これまでのように収束率が分割された区間で同じになる ように分割すると分割点は となるが、このよう 図 関数𝑓𝑓0(𝑑𝑑)のグラフに特異点が積分路に近い場合は、特異点に最も近い積分 路の点で分割することは、非常に良い近似になる。 積分区間>@の間にも積分区間に近い特異点があ る。これを使って収束率を計算すると となる。 ここでも分割することにすると、積分は区間>@、 >@、>@と分割する。ここでも収束率を使っ て計算すると となり十分精度の良い近似になっ ている。 このように分割してから二重指数型積分公式を使って計 算すると、分割なしの場合、標本点数 で、誤差 H で約 桁の精度が得られる。上のように積分区 間を 分割すると、総標本点数 で、誤差は H で約 桁の精度で得られる。 同じように分割された積分に 点のガウスの数値積 分公式を適用すると、分点数は で、誤差が H で約 桁の精度が得られた。 この計算結果は、日比野等 の計算（$4('）結果 と比較すると、この問題に関しては非常に高精度の結果 であることがわかる。日比野等では、要求精度 H の 時、標本点数 、誤差 H である。約２～３桁程 度の結果である。 日比野等で比較対象になっている $411 や '4;* のプ ログラムでも同様である。 まとめ 積分区間の近くに被積分関数の特異点があるような条 件の悪い積分を計算する場合、積分区間を分割する方法 が非常に有効であることを示した。経験的には、特異点 の近くで分割すれば、多くの場合、事足りる場合が多 い。 特異点と積分区間が決まれば積分値の計算効率を示す 収束率を定義し、収束率が大きくなるように積分区間を 分割すれば、これまでの数値積分法では実際上不可能な 数値積分も容易に計算が可能であることを示した。ここ では、そのような例として .DKQHU の 番の問題を挙げ た。 本計算では、積分路付近の特異点の位置が必要である。 このためには、被積分関数の微分係数の零点を求める必 要がある。この部分を拡充することがこれからの目標と なる。
参考文献