44 神奈川工科大学研究報告 46 るので Eulr-Mcluri の総和公式を使えば 数値積分を容易に計算できる 逆に 台形公式の計算部分を級数と見なし 積分 部分が解析的にできるか何等かの方法で簡単に計算できる場合積分の値を微分を含む部分で補正することによって 収束の遅い無限級数が計算でき る

Transcription

1 Eulr-Mcluri の総和公式を利用した数値積分 平山 加藤 43 [ 研究論文 ] Eulr-Mcluri の総和公式を利用した 数値積分 平山弘 加藤俊二 自動車システム開発工学科 Nuricl grio y Eulr-Mcluri Forul Hirosi Hiry d Sui Ko Asrc Eulr-Mcluri suio orul c rgrdd s rror vluio wi vlu o igrio d vlu o uricl igrio y rpzoidl rul. Sic diril coicis o ucio wr coid i rror vluio orul, uricl igrio ws o prord uil ow usig is. Tylor sris od wic is id o uoic diriio is usd i ordr o solv is prol, diril coicis r clcull wi suici ccurcy. is od is usd, i c pcd civ uricl-igrio orul c giv. is ppr, i is sow Eulr-Mcluri suio orul wi Tylor sris od wic giv ccurcy diril coicis cos civ uricl igrio od s s s or ldig uricl igrio d quivl grds. Sic Tylor sris od c clcul vlu o ucio i sigulr poi o ucio wi sigulriy o pprc wi suici ccurcy, i s ur wic c clcul uricl igrio o ucio wi ppr sigulr wi suici ccurcy. Kywords: Eulr-Mcluri suio orul, Tylor sris, C++ progr はじめに Eulr-Mcluri の総和公式は 周期関数の一周期に渡る積分が高精度で計算できることや 積分の積分区間の両端で関数値や微分係数が になる場合 台形公式を使うと高精度で計算できることを説明するためによく使われる公式である この公式には 被積分関数の高階微分係数を含むため 実際の数値積分に使われることは今までほ とんどなかった 高階微分係数を差分法で計算すると 桁落ちが生じ 精度の良く計算ができない このため これらの公式は積分計算に使われることがほとんどなかった 関数の微分係数を計算するために差分を使わない方法として自動微分法 が知られている この方法と原理的には同じ方法である Tylor 展開法もある これを使えば 高精度で微分係数が計算でき Tis docu is providd y JAXA.

2 44 神奈川工科大学研究報告 46 るので Eulr-Mcluri の総和公式を使えば 数値積分を容易に計算できる 逆に 台形公式の計算部分を級数と見なし 積分 部分が解析的にできるか何等かの方法で簡単に計算できる場合積分の値を微分を含む部分で補正することによって 収束の遅い無限級数が計算でき る このような研究は長田等によって行われている この場合 微分係数の計算に 数式処理システムが使われている 微分係数の計算に数式処理システムを使う場合 途中に人間が介入することになるので 少し大きな問題になると かなり大きな作業になる また途中で人的な誤りが入る可能性が生じる 計算の状況に応じて 微分係数等を計算するような場合は 実際上不可能になる Eulr-Mcluri の総和公式 数値積分や級数和の計算でよく使われる公式に Eulr-Mcluri の総和公式がある この公式は 43-4 いろいろな文献等で紹介されているが 詳しい説明がない ここでは主に長田 を参考にし その導出を行った 関数 は区間 [, ] で連続で微分可能とする 積分区間 [, ] を 等分して 台形公式を適用した計算値と厳密な積分値との差を微分係数を含む 項で近似すると 次の式が成り立つ! d ここで は roulli 数 roulli Nur であり d! である. roulli の多項式と roulli 数 roulli の多項式 を次のように定義する! の両辺に の Tylor 展開式を掛け 展開 すると! ここで!!!!!! である の式の の係数を等しいと置くと次の式が得られる 3 3 に を代入すると 入すると次のようになる この式から り返すと 次の式が得られる を代 となる この操作を繰 , roulli の多項式の定数部分を と記述し roulli 数と呼ぶ すなわち であ る を 3 に代入すると 次の roulli 数 の漸化式が得られる この式を使って 浮動小数点数で roulli 数を 計算すると 桁落ちが生じ 精度良い計算ができない悪条件の計算であることが知られている に を代入すると roulli 数 を次の ように定義することもできる! 上の式に を加え に を代入すると 4 となり 上の式は偶関数であることがわかる したがって の奇数次の係数はゼロとなり Tis docu is providd y JAXA.

3 Eulr-Mcluri の総和公式を利用した数値積分 平山 加藤 45! となる このことから が 3 以上の奇数のとき である の式に を代入すると! 5 この式は 4 から 次のように変形できる! 6 5 と 6 の の係数を比較すると 次の関係式が 得られる 7 は が奇数なら 以外ではゼロであり 偶数なら と等しくなる すなわち,, 8 となる の式の両辺を で微分すると! 8 と 8 の係数を比較することによって 次の関係式が得られる 微分の公式 を使うと roulli の多項式 を区間 [,] で 次の様に Fourir 級数に展開できる Fourir 級数の cos の係数 c は, 次の式で表される c d cos これから 次の式が得られる! cos のように奇数次数の roulli の多項式 は このようには Fourir 級数に展開できない この Fourir 級数に を代入すると! が得られる のとき 級数部分は より小さ いから 次の関係式が得られる 4! の Fourir 級数から 区間,] の時が の絶対値が最大値になるこ とがわかる すなわち [ では 3. Eulr-Mcluri の総和公式 次の積分 ここで, を定義する!, d である この式で の場合を考える この式は部分積分を使って次のように変換できる,! d! d 4! [ ] d 上の式 4 を,, について加えると 次の式が得られる, [! d ] 8 と の関係式を利用すると!, d 5 Tis docu is providd y JAXA.

4 46 神奈川工科大学研究報告 46!!! したがって 次のようになる,!, d 5 の, に 6 の関係式を代入すると! を得る ここで!,, 6 d d である 3 を使い上の式の を評価すると以下 のようになる! さらに 式を適用すると 4 d d 7 8 ここで 等式が成り立つのは 右辺の積分が のときである 多くのEulr-Mcluri の総和公式の誤差評価は 7 の式で与えているが rulli 数の大きさがわからないため直感的にその大きさがわかりにくい 8 の式は少し過大評価であるが 直感的に分かりやすくなっている この式からたとえ周期関数の一周期に渡る積分でも 微分係数が非常に大きくなる関数は 効率的に求められない事になる たとえば cos 5si d の場合 回微分する毎にだいたい5 倍位の数値になる このような場合 8 の誤差評価からあまり急速に収束するようにはならないと推定できる 分割数 でほぼ5 桁の精度が得られる もし 積分が cossi d だったら 分割数 6 で5 桁の精度が得られる 台形公式を使った場合 分割数を, 4, 8,6 と増やして計算する場合が多いので分割数が偶数になることが多い このためあまり気づかれていないが 分割数を奇数にすると非常に高精度の結果が得られる の問題の場合 分割数がでも6 の場合よりも高精度の結果が得られる 分割数が7 の場合でも6 の場合とほぼ同精度の結果が得られる の問題の場合でも同様な結果が得られる 分割数が の結果より45 の場合が良い結果が得られる 3 Eulr-Mcluri の総和公式による数値積分 Eulr-Mcluri の総和公式を使って数値積分を行うには の式からわかるように 最初に 積分区間の端点において Tylor 展開を計算する このTylor 展開を計算するプログラムは 通常の関数値の計算と宣言部分を除いてほぼ同じになる このTylor 展開を行うためのテンプレート プログラムも公開されている 6 のその計算は容易である テンプレート機能 3 を使えばTylor 展開と関数値を計算するプログラムは多くの場合 個のプログラムで記述できる 実際この論文で示した計算例は C++ 言語のテンプレート機能を使って 個の関数プログラムの形で記述した たとえば. cos cos をプログラムで記述するには 次の様に書く pl<yp T> T uc cos T& Tis docu is providd y JAXA.

5 Eulr-Mcluri の総和公式を利用した数値積分 平山 加藤 47 { } T s ; s =.*cos-cos ; rur s ; このように 個関数プログラムを書けば 関数計算やTylor 展開の計算に使える 次に分割数 を決める 分割数がある程度推定できる場合 その値にする 推定が出来ない場合 最小値の にする この分割数を利用して 台形公式を利用して 積分の近似値を計算する 近似値を次の項を利用して 補正する この級数は漸近級数なので漸近級数の手法を使って計算する c! この補正項 c の絶対値が要求精度以下になるま で 補正項を加算して行く 補正項の絶対値が要求精度より小さくなったら その項を使って補正し計算を止める この計算値が求める積分値になる 補正項の絶対値が要求精度より小さくならない場合や逆に大きくなった場合 補正するのを止め 分割数 を 倍にして 台形公式による積分の計算に戻る 分割数を 倍にすると その以前の台形公式計算値を利用し 計算量を少なく出来るためである この手順を繰り返す事によって積分値を計算する 3. 簡単な計算例 簡単な例として 次の積分を計算する これを 要求精度 で計算する d 最初に 積分区間の両端点で Tylor 展開をする 展開は 次まで行った その計算結果を以下に示す = における Tylor 展開を 5 次まで表示すると次の様になる における Tylor 展開を 5 次まで表示すると次の様になる として 台形公式で数値積分すると となる これを補正する この場合の補正項は収束が非常に遅く要求精度を満たすことができない 8 で逆に絶対値が増加し始める このため分割数を 倍の 4 として 台形公式を使って数値積分を行う このときの積分値は となる この結果を再度補正する ここでのは補正項も十分速く小さくなるとは言えないが 7 で補正項が 3. となり要求精度より小さくなった そこまでの計算値が積分値となる 積分値は となった 要求精度通りの結果が得られた このときの関数計算回数は 5 回でその内 回が Tylor 展開を計算するための関数計算であった 分割数 を増加させると は小さくなるので補正項 c はすばやく小さくなる このため 高 次の微分係数はあまり必要としなくなる Tylor 展開は低次数でも Eulr-Mcluri の総和公式を利用できることになる 逆に Tylor 展開が高次数であるならば 分割数 を小さく出来る 上の計算では 次の Tylor 展開を利用したが それが適切かどうか考慮する必要がある 3. 数値例 Eulr-Mcluri の総和公式を利用した計算法の性能を評価するために Kr の問題 8 の中から両端点で特異性を持つため Tylor 展開できない問題等を除いた 3 問題について計算を行った 計算した問題を表 に示す 番号は Kr の問題の番号である それを計算した結果を示す 実行結果を表 に示す N は分割数すなわち関数の計算回数 Error はその積分ルーチンが出力した誤差である 分割数の内 回は Tylor 展開の計算である EM は Eulr-Mcluri の総和公式を利用した計算 DE は二重指数型積分公 5 Doul Epoil orul DAQN は 適応型ニュー 8 トン コーツ法の結果はその参考文献による結 Tis docu is providd y JAXA.

6 48 神奈川工科大学研究報告 46 果である M 型 64 ビットの浮動小数点演算による結果である DE のプログラムとして ネットで公開されている大浦 の C 言語用 DE プログラムを C++ 言語に変換し使用した コンパイラーとして Microso Visul C++ を使用した 4 まとめ Eulr-Mcluri の総和公式を利用して 数値積分する方法を提案した この方法は 単純でわかり易い計算法であるにも関わらず 多くの場合 有力な数値積分法と同等程度またはそれ以上の性能を発揮する この計算を行うには 被積分関数の微分係数を精度良く計算する必要がある この方法として 自動微分法があるが 一般のユーザーにそれほど普及していない この方法が使われるためにも自動微分の普及することを期待したい 5 参考文献 []Arowiz M. d Sgu. A.: Hdoo o Micl Fucios,Dovr,7 [] Dvis P.J., iwiz P. 森正武訳 : 計算機による数値積分法, 日本コンピュータ協会, 8 [3] Dvid Ars, Alsy Gurovoy: C++ Tpl Mprogrig, Addiso Wsly, 5 [4] Gisl Egl-M ullgs, Fr Ulig: Nuricl Algoris wi Forr, Sprigr, 6 [5] Hrici P.: Applid d Copuiol Copl Alysis, Vol., Jo Wily & Sos, Nw Yor, 74 [6] 平山弘 館野裕文 浅野直之 川口隆史 : Tylor 級数演算ライブラリの使用法, 東北大学情報シナジーセンター大規模科学計算システム広報 SENAC, [7] 平山弘 小宮聖司 佐藤創太郞 :Tylor 級数法による常微分方程式の解法, 日本応用数理学会論文誌, pp. 8 [8] 二宮市三 : 適応型ニュートン コーツ積分法の改良, 情報処理学会誌,8,54 53 [] 大浦拓哉 : Oour s Micl Sowr Pcgs, p:// [] 長田直樹 : 対数収束級数の漸近展開と加速法, 情報処理学会論文誌 88, 56 6 [] 長田直樹 : お話 : 数値解析第 3 回オイラー マクローリンの公式, p:// 8-7.pd [] ll,l..: Auoic Diriio Tciqu d Applicios, Lcur Nos i Copur Scic, Vol., Sprigr Vrlg, rli-hidlrg-nw Yor, 8 [3] 篠原能材 : 数値解析の基礎, 7 章, 日新出版, 78 [4] 杉原正顕, 室田一雄 : 数値計算法の数理, 章, 岩波書店, 7 [5] Tsi, H. d Mori, M.: Doul poil orul or uricl igrio, Pul. MS, Kyoo Uiv., 74, 4 Tis docu is providd y JAXA.

7 Eulr-Mcluri の総和公式を利用した数値積分 平山 加藤 4 Tl Ts Prol No. grl d cos cos d d d d si d d d si34.5 d d si5 5 d coscos 3si cos 3si 3cos3 d d Tis docu is providd y JAXA.

8 5 神奈川工科大学研究報告 46 Tl Copriso o Prorc o Qudrur ouistolrc = EM DE DAQN No. N Error N Error N Error Tis docu is providd y JAXA.