Byrne型Pseudo-Runge-Kutta法の特性について利用統計を見る

(1)

論文

Byrne型Pseudo−Runge−Kutta法の特性について

（平成元年8月31日受理）小野俊治田中正次山下茂

On Byrne's type of Pseudo-Runge-Kutta methods

ToshiharuONO MasatsuguTANAKA ShigeruYAMASHITA

Abstract Byrne and Gruttke presented theL optimal parameters for 3rd，4th and 5th−order Pseudo −Runge−Kutta methods． Concerning the same methods We show the regions of absolute stabi1・ ity and the contour maps of criteria on the size of truncation error． By the observation of their maps， we claerly see that the optimization of the above authors are proper， and what charac −terristics the their formulas have． 1．はじめに 1966年ByrneとLambertは、 R皿ge−Kutta法のstep当りの関数計算回数の削減を主眼としたPseudo−Runnge −Kutta法を提案した1）。これは一段解法（Runge−Kutta）及び線形多段階法のいずれの解法にも属さないが両者の特徴をもつ解法で、Costabile、田中、新谷、中島らに引き継がれ改良あるいは発展を遂げ、今日に至っている5）・6）・7）・8）。Byrneが提唱した解法は2段階法であるt）・2）。彼及びGruttkeは、微分方程式が単一の場合について打ち切り誤差の観点による2段数、3段数及び4段数法の最適公式を導いている2）ほ）。ここでは両者が報告した最適公式の正当性の視覚的な確認と、別の基準による連立及び単一の微分方程式に対する最適公式の導出（使用計算機はACOS 850）について報告する。

2． Byrne型2段階法の一般形1）

連立の常微分方程式の初期値問題を

dyi

＝Fi（y（x）） x∈［a，b］

dx

yi（a）＝ηi （i＝1，2，…q）とする。ここで、 y1（x） y2（x） y（x）＝：： y。（x）（2．1．1）（2．1．2）である。初期値問題の理論解Y（x）は一意でクラスcm＋2 に含まれるものとする。但し、mは解法の次数でここでは3∼5である。またY（Xn），x．C［a，b］はx＝Xnにおける理論解で、ynはその数値解をあらわし、刻み幅は一定でh＝Xi−Xi．iを満たすものとする。 ynからy 。＋iを求めるために必要な1ステップ間での関数F［Y （x。）］の評価回数を段数と呼ぶ。p段数の一般公式は（2．2．1）で表せる。 ko（yn）＝hF［yn］ k，（yn）＝hF［yn＋μ1k◇（yn）］ k2（Yn）＝hF［Yn＋（μ2一μ2、1）k◇（yn）＋μ211kt（Yn）］ kp＿1（yn）＝hF［Yn＋（μP＿t一μP＿t、t−… 一μP−t，P−2）ko（Yn）＋SL，，t，、k、（Yn）＋… ＋μP−11P−2k，−2（Yn）］ P−t P−1 k（Yn）＝Σαiki（Yn）＋Σ βik，（Yn−t） i＝o i＝o yn刊＝yユ＋k（yn）（2．2．1）ここで、 Ft（y（Xn）） F2（y（Xn）） F（y（Xn））＝、：： F，（y（xぬ））（2．2．2） αi，βi（i＝1，2，3），μj，μk，．（j，k，r＝1，2，3，k＞r）は定＊電子情報工学科，Department of Electrical Engineering and

(2)

数で、変数XがXnからX。＋tに増すときの理論解の増分関数D（Y（Xn））は D（Y（Xn））＝Y（Xn＋1）−Y（x丑）（2．3．1）また、このときの局所打ち切り誤差T［Y（x。）］は T［Y（xn）］＝k（Y（xn））−D（Y（xn））（2．3．2）で与えられる。増分関数D（Y（Xn））が、（2．2．1）に示される数値解の増分関数k（Y（x。））とhqの項まで一致するとき（2．2．1）で示されるPseudo−Runge−Kutta法は q次法であるといわれる。一31＋40（μ1＋μ3）−50μtμ3 β2＝ α2＝一β2， 120μ2（μ2一μ1）（μ2一μ3）， β3＝一一31＋40（＋）50 120μ3（μ3一μ1）（μ3一μ2）， βo＝一（βt＋β2＋β3）−1／2，』 μ2，1＝

一31＋40

240β・2μi（μ2一μ3）’ α3＝一β3， α◇＝1一βo， 3．解系 Byrne型2段階法はP段数でP＋1次の精度を得ることが知られている（P＝2，3，4）D・3）。5次法yn＋tの次数条件式を求めるために必要なD（Y（x。））とk（Y（Xn））の展開の係数を表1．に示す3）。記法はByrneが使用したものと同じであるD・2）。 Fy＝∂F（Xn， yn）／∂y F，j＝∂F／∂Y」， F，，j＝∂2F／∂Yi∂Yj また、繰り返された添え字は和をとることを表す。 9_{F，．jVj＝Σ （∂Fi／∂Yj）Vj} j＝1 P段数P＋1次法の解系は次のようになる。（1） 4段数5次法表1．よりμt，・PL・2を自由パラメータとした次の5次法の解系が得られる。

＿＿一

β1 ， α1＝一β1， 120μt（μビμ2）（μ1一μ3）表1． 5次法誘導のための係数

μ・it＝｛＿1（ll）

@μ1［−31＋40（μt＋μ2）−50μ1μ2］｝

×｛＝＋一『L｝

＿ −31＋60μt μ3・2−@ 360μ2β3（μ2一μt）・ μ3＝62／85 （3．1．1）但し、μ1≠0，μ1≠62／85，μi≠31／60，μ2≠0， μ2≠62／85，μ2≠（40μ仁31）！（50μ仁40），μ1≠μ2 とする。（2） 3段数4次法表1．においてα3，β3，μ3，μ311，μ312＝0とおくこび_{とによりμt，μ2を自由パラメータとした次の4次法の} 解系が得られる。 α1＝

u可・ β・＝一…

項係数 D（Y（x。）） k（Y（x，））

hF

1 α0＋αt＋α2＋α3＋β0＋β1＋β2＋β3 h2F， 1！2 α1μ1＋α2μ2＋α3μ3一βo＋β1（μ1−1）＋β2（μ2−1）＋β3（μ3−1）（h3／2）F2） 1！3 α1μ12＋α2μ22＋α3μ32＋β日＋β1（μ1−1）2＋β2（μ2−1）2＋β3（μ3−1）2 （h3／2）F，。F㌶ 0 一（α1＋β1）μ12＋（α2＋β2）（2μ2．1μ1一μ22）＋（α3＋β3）（2μ3．1μ1＋2μ3．2μ2一μ32）（h4！6）F《3， 1！4 α1μ13＋α・μ23＋α・μ33一β・＋β1（μ「1）3＋β，（μ2−1）3＋β，（μ3−1）・（h4／6）F，。F3（2） 0 一｛α1μ13＋β1μ12（μ「3）｝＋｛α・2（3μ2．！μ12一μ23）＋β2［3μ2，1μ1（μ「2）一μ22（μ2−3）］｝＋｛α3［3μ3，1μ12＋3μ3．2μ22一μ33］＋β3［3μ3．1μ1（μ「2）＋3μ3，2μ2（μ2−2）ごμ32（μ3−3）］｝（h4／2）F．stFs，Ft 0 一｛α1μ13＋β1μ12（μ1−1｝｝＋｛α2（2μ2．1μ「μ22）μ2＋β2（2μ2，1μ1・・μ22）（μ2−1）｝＋｛α3（2μ・，1μ1＋2μ・，2μ2一μ32）μ・＋β，（2μ3，1μ1＋2μ，．，μ2一μ32）（μ3−1）｝（h4！2）F，sFs，aF8， 0 一（α・＋β・）μ・，1μ12−（α3＋β・）（2μ・．1μ12−2μ，，2μ，．1μ1＋μ，．，μ22）（h5／24）F“） 1！5 α1μ14＋α2μ24＋q3μ34＋βo＋β1（μド1）4＋β2（μ2−1）4＋β3（μ3−1）4 （h5！24）F．。F。《3） 0 一｛α1μ14＋β1μ12（μ12−4μ1＋6）｝＋｛α2（4μ2．1μ13一μ24）＋β2［4μ2，1μ1（μ！2−3μ1＋3）一μ22（μ22−4μ2＋6）］｝＋｛α3［4μ3．1μ】3＋4μ3．2μ23一μ34］＋β3［4μ3．1μ1（μ12−3μ1＋3）＋4μ3．2μ2（μ22−3μ2＋3）一μ32（μ32−4μ3＋6）］｝（h5／6）F．3tFsl2）Ft 0 一｛α1μ14＋β1μ12（μ1−3）（μ1−1）｝＋｛α2（3μ2．1μ1・2一μ23）μ2＋β2［3μ2，1μ1（μ1−2）一μ22（μ2−3）（μ2−1）］｝＋｛α3［3μ3．1μ12＋3μ3．2μ22一μ33］μ3 ＋β3［3μ3．1μ1（μ・ド2）＋3μ1＋3μ3．2μ2（μ2−2）一μ32（μ3−3）］（μ3−1）｝（h5！8）F．。tF。’Ft’ 0 一｛α1μ14＋β1μ12（μ12−4μ1＋2）｝＋｛α2（4μ2，12μ12一μ24）＋β2［4μ2，1μ！（μ2．tμ「2μ2＋1）一μ22（μ22−4μ2＋2）］｝＋｛α3［4（μ3，1μ1＋μ3．2μ2）2一μ34］＋β3［4（μ3．1μ1＋μ3．2μ2）（μ3．1μ1＋μ3，2μ2−2μ3＋1）一μ32（μ32−4μ3＋2）］｝（h5！4）F．st．Fs’FtF、 0 一｛α1μ14＋β1μ12（μド1）2｝＋（α2（2μ2．1μドμ22）μ22＋β2（2μ2，1μ1一μ22）（μ2−1）2｝＋｛α3（2μ3．1μ1＋2μ3，2μ2一μ32）μ32＋β3（2μ3．1μ1＋2μ3．2μ2一μ32）（μ3−1）2｝（h5／2）】F．sF5，abFa，Fb 0 一｛α・μ・．1μ13＋β2μ2．1μ12（μ「1）｝一｛α3（μ3．1μ13＋μ，，，μ23−2μ3，2μ1μ，）＋β3［μ3，1μ12（μド1）＋μ3．2μ22（μ2−1）−2μ3．2μ3．1μ1（μ2−1）］｝（h5／6）F．sF。．。Fa《2） 0 一｛α・μ・．1μ13＋β2μ2，1μ12（μ1−3）｝一｛α・（μ3，1μ13＋μ，，，μ23−3μ，．2μ2．1μ12）＋β3［μ3パμ12（μ1−3）＋μ3．2μ22（μ2−3）−3μ3，2μ2．1μ1（μ「2）］｝（h5／2）F，ヨtFsFt．aFa，’ 0 一｛α2μ2．1μ12μ2＋β2μ2．1μ12（μ2−1）｝一｛α3（μ3，1μ12＋μ3，2μ22−2μ3．2μ2．1μ1）μ3 ＋β3［μ3．1μ12＋μ3，2μ22−2μ3．2μ2．1μ1］（μ3−1）｝

(3)

α2 ＝ 5μ1−4 12μ2（μ1一μ2）， β・＝一．1．LE．21，1la：．IE．，．L一 μ2・t＝@ μt（4−5μ1） β2＝一α2， αo＝1一βo （3．1．2）但し、μ1≠0，μ2≠0，PL・t≠μ2，μ1≠4／5とする。（3） 2段数3次法同様にして表1．よりμ1を自由パラメータとした次の 3次法の解系が得られる。 βo＝ 5一μ1 12μ1 iB ，ニーαt

…＝1一β・・α1＝

ﾅ1・

但し、lt．t≠0 （3．1．3） 4．打ち切り誤差先ず、記法を次のように定める。

D＝∂／∂x十F∂／∂y，F＝F（Xn， yn），

k段数p次Pseudo−Runge−Kutta法において、微分方程式が単一及び連立の各場合の打ち切り誤差の主項をそれぞれγ。hp＋t及びテ。hp＋1とすれば、 γ3＝a3tD3F＋a32F予D2F＋a33Fr2DF＋a34DFDF7 （4．1．1） γ4＝a41D4F＋a42FyD2F＋a43F72D2F＋a44Fr3DF ＋a45DF7D2F＋a46D2FyDF＋a47FyDFDFy ＋a48Fy7（DF）2 （4．1．2） γ5＝astDsF＋a52F7D4F＋a53DFアD3F＋a54F72D3F ＋assD2F7D2F＋as6FyDFyD2F＋as7F73D2F ＋a58（DFr）2DF＋a5g F 72DFrDF＋a5toF74DF ＋a51tF7yDFD2F＋a512DFyy（DF）2 ＋a5t3F7F77（DF）2＋ast4F，D2FyDF＋astsD3FyDF （4．1．3） γ3＝a3tF（3）＋a 32F， jkFj Fk ’＋a33F．jFjtt （4．2．1） 74＝a4tF（4）＋a42F，jmFjFm，，＋a4，F，jkmFjFkFm，＋a44F，jFj（3）＋a4sF，jFj．kmFx，Fm＋a46F，jk ×Fj，Fk，＋a47F．jmFj．kFk，Fm＋a4sF，jFj．kFk，（4．2．2）てア5＝a5tF（5）＋a52F，sFs〈4》＋as3F．etFs（3）F，＋含54F，。tF，，，Ft＋a55F．stnFs（2）FtFn ＋as6F．stuFs，Ft，Fu＋a57F，stuvFs，F，F竃Fv ＋a5sF．。F。．ab。Fa，FbF。＋a sgF，。F。， abFa，Fb，＋astoF．sF。．abFa，iFb＋asttF．。F。， aFa（3）＋ast2F．stF。Ft，abFa，Fb＋as13F，。tF。Ft．aFa” ＋白5t4F，。tF。 ’Ft ，aFa，＋a5t5F，。tuF。FtFn．aFa，＋a5t6F，。F。，］F1，abF。，Fb＋ast7F，。Fs．］Fi．a ×Fa，，＋a5t8F．，F。．1mFIFm．aFa，＋a5tgF，．tFe XF，，1FLaFa，（4．2．3）と表すことができる。ここで，aij，a．tは公式の係数のみの関数である。次に、微分方程式が単一及び連立の各場合について、各次数の公式の打ち切り誤差の大小を判定するための尺度として次の5種類の判定基準を使用する。r次法の単一及び連立の基準をそれぞれA，i，▲riで表し、これらを公式の打ち切り精度判定基準と呼ぶ。これらの基準が各公式の打ち切り誤差の大小の判定に有効なことについては多くの実績がある2）’4）15）。

A31＝81a311＋31a321＋31a331＋41a341，

4 4

A32＝Σla3il， A33＝Σ（a3i）2，

’；1 ’；1 ゑ32＝￨、la・・1・▲33＝蓬、（a・・）2（4．3．1）

A4t＝161a411＋81a421＋31a431＋31a441

＋81a45 1＋8 1 a46 1＋4 1 a47 1＋4 1 a48 1 ，

A42＝Σla4il， A43＝Σ（a4i）2，

i；1 i＝1 8 A・・＝O11a・・1・A・・＝三1（a4i）2（4．3．2）

A51＝161a511＋81a521＋81a531＋41a541

＋81a55 1＋4 1 a 56 1＋2 1 a 57 1＋4 1 a 58 1 ＋21a5g l＋la5io l＋41a51t l＋41a5i2 1 ＋21a513 1＋41ast4 1＋81a515 1，ユら

A52＝Σ la5il，

i；1 ゑ52＝Σ la5il， i；1 5．安定性 tら A53＝Σ 1；’ ゑ53＝Σ i；1 （a5、）2 （a5i）2 （4．3．3） Pseudo−Rnge−Kutta法を用いて数値解y。を求めるとき、ステップが進むにつれて誤差が減少し、数値解が理論解に次第に近づいて行くことを安定といい、その逆を不安定という。連立常微分方程式の解法の安定性は、単一の微分方程式の解法の安定性に帰着する。テスト方程式 y’＝λy （λは負の実部をもつ複素数）（5．1．1）に（2．2．1）を単一にした場合の式を適用することより Yユ＋t＝Hllp（hλ）yn十H2．p（hλ）yn＿t （5．2．1）が得られる。ここで、Ht，p（hλ），H2，p（hλ）はB，を P−1 Bt（α，P）≡1＋Σαi 但し、 P≧2（5．2．2） i；O P−t B2（α，P）≡Σμiαi P≧2（5．2．3） i＝1 P−1 P−2 B3（α，P）≡Σ Σμ」μi、jαi P≧3（5．2．4） i；2 」＝t i）j B，（α，P）≡ μtμ2，1μ3，〉αp＿t P≧4（5．2．5）とすると、 PHt，p（hλ）＝Σ ｛Bi（α，P）（hλ）i｝ ’Bt H2、p（h》し）＝Σ ｛Bi（β，P）（hλ）i｝−hλ（5．2．6） i＝t と表すことができる。（5．2．1）のY。＋t、Y。及びyn−t の代りに真の解Y（Xn＋1）、 Y（x。）及びY（Xn．t）を代入すると Y（Xn＋t）＝H1．p（hλ）Y（Xn）＋H2．p（hλ）Y（Xn−t）一一 Tn＋t （5．3．1）が得られる。ここでT。刊はYn＋1の局所打ち切り誤差である。（5．2．1）から（5．3．1）を辺々引算しen＝ YD−Y （Xn）とおき、 Tn＋tの代りにlT。＋tl≦Eを満足するEを代入すれば、集積打ち切り誤差e。の差分方程式 en＋t＝Ht．p（hλ）en＋H2．p（hλ）en＿t十E （5．3．2）が得られる。（5．3．2）の同次形に対しe。＝ξとおくことにより、P段数の特性方程式（5．4．1）を得る。 ξ2−Ht．p（h；し）ξ一H2．p（hλ）＝0 （ξ（≠0）は複素数）（5．4．1）さらに（3．1．1）∼（3．1．3）の解系を（5．2．2）∼（5．2．5）に代入することよりHi，p（hλ）はμj及びμk，。を含まない

(4)

ことがわかる。よって、絶対安定領域は自由パラメータに関係なく一意に定まることがわかる。特性方程式（5．4．1）の2根をξ，，ξ2とするとe。に関する差分方程式の一般解は e。＝Ctξ tn＋C2 9 2n＋｛（5．3．2）の特解｝と表される。従って、n→。。としたときe。が有界であるための必要十分条件は 1ξ11＜1，1ξ21＜1 である。（1） 4段数5次法の場合絶対安定領域S5は（3．1．1），（5．3．1）及び（5．4．1）より、 S5＝｛Zl720ξ2i −Ht．4（Z）ξi−H214（Z）＝0， Hl，4（Z）＝720＋1080Z＋300Z2＋120Z3＋31Z4， H2，4（Z）＝一（360Z＋300Z2＋120Z3＋31Z4）， 1 ξil＜1 （i＝1，2）， Z（≡C｝（5◆5◆1）で表せる。Z空間におけるS5を図示すると図1の閉曲線PRK5の内部になる。また有効絶対安定領域に関する数値を表4に示す。（2） 3段数4次法の場合絶対安定領域S4は、α3＝β3＝μ3＝μ3，t＝・PL・3，2 ＝0と（3．1．2），（5．3．1）及び（5．4．1）より S4＝｛Zl 12ξ2i−Hll3（Z）ξi−H2，3（Z）＝0， Hi，3（Z）＝12＋18Z＋5Z2＋2Z3， H2，3（Z）＝一〔6Z＋5Z2＋2Z3〕， 1ξil＜1 （i＝1，2），Z∈C｝（5．5◆2）で表せる。Z空間におけるS4を図示すると図1の閉曲線PRK4の内部になる。また有効絶対安定領域に関する数値を表3に示す。（3） 2段数3次法の場合絶対安定領域S3は、α2＝β2＝μ2ニμ3＝0と（3．1．3），（5．3．1）及び（5．4．1）より S3＝｛Zl12ξi2−Ht，2（Z）ξi−Hz2（Z）＝0， Ht，2（Z）＝12＋3Z＋5Z2，H2，2（Z）＝一〔6Z＋5Z2〕， 1ξil＜1 （i＝1，2），Z∈C｝（5．5．3）で表せる。Z空間におけるS3を図示すると図1の閉曲線PRK3の内部になる。また有効絶対安定領域に関する数値を表2に示す。図1 Byrne型2，3，4段数の絶対安定領域 i pRK5 一■・・一”t・鼈黶h ’ @ ’ @，’噤f f

絢〆

、i ’i 0．5

PRK5

’−

oRK4

∼←“ ■ ← ← ← _← 3．0 ．・Q．0 一1．0 実軸一．・A・・、 @ 、、／一〇．5

PRK3

ノi ’．、、＼ i P

ｺノ

“ “舶一一・←一・・…

PRK5

1虚軸 6．丸め誤差に関する性質の判定 2段数3段数及び4段数の各公式の丸め誤差に関する特性を判定する尺度として、次に定義する数量を使用する。 t 1 R3＝Σ ｜αil＋Σ 1β11＋21μ11， ’；° ；＝° （6．1．1）

R4＝Σ1αi1＋Σ1β，1十21μ11十21μ21

i＝o i；o 十41μ21il，（6．1．2） 3 3 3

R5＝Σ1αi1＋Σ1βil＋2Σ1μiI

i＝o iニo i＝t ＋41μ2、il＋41μ3、il＋41Pt 3，21 （6．1．3）これを2段数，3段数及び4段数の各公式に対する丸め誤差特性判定基準と呼ぶ。丸め誤差にとって、ki（y n）やki（Yn．t）にかけられる係数が小さい程好ましいことは自明であろう。通常のように常微分方程式の数値解が倍精度で計算され、また打ち切り誤差が丸め誤差と比べてかなり大きい場合には、丸め誤差特性判定基準が相当大きくてもその影響は軽微である。しかし打ち切り誤差が非常に小さく丸め誤差と同等の大きさを持つ場合や、判定基準が著しく大きい場合に一これは最適化によって顕著な打ち切り精度の改善がみられる場合にしばしば発生する一は警戒が必要である。 7． Byrne型Pseudo−Runge−Kutta法の特性 2段数3次法の打ち切り精度判定基準と丸め誤差特性判定基準の両者の等高線を図2に示す。この図では横軸に自由パラメータμ1を、縦軸に各判定基準の値をとっている。

図2

4．0 Byrne型2段数3次法の判定基準A31 と丸め誤差特性判定基準R3の等高線図 i丸め誤差特性判定基準R3 i 公式・2i 公式33 i公式31 打ち切り精度判定基準A3， 0．2 0．4 0．6 0．8 1・0 1・2 1．4 1令6 1．8 2．0 ロ 3段数4次法の打ち切り精度判定基準の等高線と丸め誤差特性判定基準の等高線の両者を図3に、また4 段数5次法のそれを図4示す。図3，4では自由パラメータ Pt・tとμ2をそれぞれ横軸、縦軸に選び打ち切り精度判定基準の等高線を実数値の付いた太い線で、丸め誤差特性判定基準の等高線を整数値の付いた細線で示している。

(5)

図3 Byrne型3段数4次法の判定基準A41と丸め誤差特性判定基準R4の等高線ミ6 8．既知公式現在文献等で紹介されているByrne型の公式は、著者の知る限りではByrne自身による2段数3次法、3段数4次法とGruttke による4段数5次法のみであるので、これを特性の解明と評価の対象とする2），3）。表 3∼5にByrne型の各公式の特性値を示す。 9．公式の導出 Byrne型の場合、安定性は自由パラメータに無関係に一意に定まるので、図2∼図4 を参考にして打ち切り精度判定基準について特に有意義な公式を選び、ついで丸め誤差の観点を考慮して公式を導いた。Byrne 型P次法の打ち切り精度判定基準Ap2，Ap3 ，▲p2，▲p3が最小となる自由パラメータから得られる公式を順に公式p1，公式p2，公式p3，公式p4と名付ける。 2段数3次法の場合、丸め誤差特性基準は表2に示すように低く、打ち切り精度判定基準値が最良となる公式31∼公式34を得る。但し、公式33はμ＝0．8でByrne3 （Byrneによる3次法）と一致する。 3段数4次法の場合、丸め誤差に関する性質が劣化せず、打ち切り精度判定基準A 42，A43，ム43が最良となる公式41，公式42 ，公式44を導いた。打ち切り精度判定基準．o．30 ，0．10 ．．．50 ．．・10 A42を最良とする公式，3は，丸め誤差の点 μ1 で問題になるので、打ち切り精度を若干犠図4 Byrne型4段数5次法の判定基準As 1と丸め誤差特性判定基準R5の等高線牲にし、丸め誤差特性判定基準を小さくした公式，3°PTに改良する。 N μ1 4段数5次法の場合、Gruttkeが使用した手法3）（荒い刻み幅で最良区間を絞り込み最良点を見つけ出す手法、以降区間縮小法と呼ぶ）により、判定基準A52，A53，As2 ，A53がそれぞれ最小となる4公式の公式5 1，公式52，公式53，公式54を導いた。 10．数値実験と考察誘導が正しく行われたことを確認するために、数値実験を行った。数値実験の結果を表5に示すが、紙面の都合上4段数5次法のみとする。 10．1 2段数3次法（1）各公式の性能はほぼ同じで優劣付けがたい結果になっている。 10．2 3段数4次法（1）Byrneは、 A4t（Lotkinの誤差限界の公式の係数のみに依存する係数）を判定基準にして打ち切り精度最良な公式By rne4としてpt．t＝0．541，μ2＝0．763を選んでいるが、図3によればこの判断は適切であることがわかる。しかも、丸め誤差特性基準も最良に近いので極めて好ましい公式といえるだろう。

（2）公式41は、表3から明らかなよう

に打ち切り精度判定基waA，，においても最小であるから、公式Byrne4の最適化でもあ

(6)

る。実際、数値実験を行うとわずかながら公式，1が優れている場合が多い。（3）公式，2と公式，4は、Byrne4および公式41より常に若干悪いことが観察され不適当である。（4）公式，3は、打ち切り精度をよくするために、 PL・tsμ2を共に0の近傍に取っているが、丸め誤差に関する性質が極端に悪化するので不適当である。』このことは数値実験結果からも（刻み幅を小さくすると他の公式と比較して誤差が大きくなっていることから）確かめられる。（5）公式，3を丸め誤差の観点から最適化した公式4 3°PTは、非線形連立微分方程式において他の公式よりも優れている例を見つけることができ、この点で有効と考えられる。 10．3 4段数5次法（1）Gruttkeの区聞縮小法では打ち切り精度判定基 waA，，を最小にする点を見逃す恐れがある。事実、表 4から明らかなように公式51の方が小さいし、よい結果が得られる場合が多くみられる。（2）数値例の選び方にも問題があるかも知れないが、微分方程式が単一であるかあるいは連立であるかを問わず、公式52，公式54の場合によい結果が得られている。これらの公式は桁落ちの危険の少ない安定な公式でもある。

11．まとめ

（1） Byrne型2段階法の2段数3次法、3段数4次

法及び4段数5次法は、安定性が自由パラメータに関係なく一意に定まる方法である。従って、公式の自由パラメータによる特性の制御は、打ち切り精度と丸め誤差の両面に関して行われることになる。

（2）Byrne型 2段階法の2段数3次法、3段数4

表2 Byrne型2段数3次法の特性凡例は各判定基準での最小値を示している。打ち切り精度判定基準安定性公式名自由 pラメータ @ μ A31 P一，1ρtkin A32 P一，Σll A33 P一，Σ02 A32

A立，Σ1｜ A33A立，Σoε 丸め誤差ﾁ性判定譓?l 有効絶対安定ﾌ域の面積区間の @ 長さ選出判定 @基準 Byme3 0．8 ₅₇₃ ki666ハ 112 i＝0．5） 1／12 i＝0．0833…） 1！3 i司．3333…） 1！18 i＝0．0555…） 3．6 A31 公式31 1．2 _5鴻 i1666）蹴2

_癌ﾚ

1／16 i司．0625） 7！12 i司．5835…） 19！144 i＝0．1319…） 4．70555 A32 公式32 1．12 _5灘 i1666） 13！30 i＝0．43333…） 11z180 i二〇〇6i1〕 8！15 i司．533訣・・） 17／150 i＝0．1133…｝ 4．24 3．90 2．26 A33 公式34 4！15 31！9 i＝3．444…） 516 i司．83333…） 71！324 i＝0．2191…）万9_i間泌） 1／54 i灘0185） 3．78333 A33 表3 Byrne型3段数4次法の特性凡例は各判定基準での最小値を示している。打ち切り精度判定基準安定性公式名自由 pラメータ纈i：μ1 ｺ段：μ2 A41P一，1ρtkin A42 P一ΣlI A43

P一，Σ02 A42A立，Σll ん3A立，Σ02 丸め誤差ﾁ性判定譓?l 有効絶対安定ﾌ域の面積安定区間ﾌ長さ選出判定 @基準 Byme4 ．541 D763 ．63193611d＋00 ．17511583d＋00 ．10781563d−Ol ．14848403d＋00 ．61427398d−02 7，073 A41 公式41 ．54229 D76219 631幽亘σo 1750髄1ユ翻oo ．10769081d−01 ．14909971d＋00 ．61737754d−02 7，073 A42 公式42 ．87061 D76488 ．88769062d杖）0 ．17846662d＋00 61949697dO2 ．28145192dH）0 ．17119794d−01 7，421 A43 公式43 一．00001 D0㎜2

．64582511d＋Ol ．10333211d＋Ol ．21132573d＋00

_43⊇OI

．18537330d・02 ．3333d＋10 _3．39 _1．34 A42

公式430PT ．13000

D78000

．13944444d＋Ol ．30294444d＋00 ．28153621d・・Ol ．99944444d−01 ．25246150d・02 13．43 A42，R4 公式44 ．01425 D25939 ．47944520d＋01 ．79587453d＋00 ．12438385d＋00 ．64003398d−01 _蛯248蹴oz 149．8 A43 表4 Byrne型4段数5次法の特性凡例は各判定基準での最小値を示している．打ち切り精度判定基準安定性公式名自由 pラメータ纈i：μ1 @ ： 2 As 1 P一 Lotkin A52 P一Σll A53 P一Σ02 A52 ｧΣll A53 ｧΣ（2 丸め誤差ﾁ性判定譓?l 有効絶対安定 @ の面安定区間ﾌさ選出判定 @基準 Gruttke5 ．250 D869 2．3069 0．3073 0．01654 0．1527 0．2676 16，569 A5t 公式51 ．225 D789 2096正 02894 0．01442 0．09929 0．001045 17，892 A52 公式52 ．410 D839 2．4416 0．3006 001348 0．2290 0．008505 16，429 _4．09 _2．87 A53 公式53 ．061 D406 3．3991 0．4396 0．02954 0．04537 0．0003069 16，856 A52 公式54 ．085 D444 3．3081 0．4283 0．02859 0．04653 00002819”∴’ 13，176 A53

(7)

次法及び4段数5次法について、自由パラメータが実用的に有意義な変域において変動するとき、打ち切り精度判定基準及び丸め誤差特性判定基準がどのように変動するかを図式的にとらえた。このグラフの利用により我々は、さらに性能の良い公式開発の可能性の有無の判定や、任意の公式が与えられたときその特性の直観的な評価等が容易に出来る。（3）打ち切り精度判定基準にはいろいろな尺度があ表5．数値実験の結果 1．初期値問題理論解 y，＝−y十sin （2x）， sin（2x）−2cos（2x）

y−

@ 5 り、全ての尺度に対して最良な公式は存在しない。（4）Byrne型の3段数4次法の公式は、他の2段数3 次法及び4段数5次法の公式に比べ安定区間が非常に短いので刻み幅を小さくとる必要がある。（5）既知公式Gruttke5（Gruttkeによる5次法）は、図4より打ち切り精度の面で改良可能であることがわかる。公式誘導に当たっては、公式周辺の丸め誤差特性判定基準の動向に十分注意する必要がある。 y（0）＝−O．4 （一

ﾌく・）

ステップ数100 誤差（刻み幅h＝0．1）誤差（刻み幅h＝0．05）誤差（刻み幅h＝0．01）公式名選出

ｻ基

1st ste Hax． 1ast 1st ste Hax． 1ast 1st ste Ha叉． 1ast

Gruttke5 A51 ．562691−09 ．301976−06 ：i：：ユ姶951：≒仰1 ：i：：マ聯；ユOi・．943561−08 ．476634−08 ．81田1揮・14．．248496−11 ．248496−11

公式51 A52 ．561468−08 ．270345−06 ．357403−07 ．143734−09 ．844042−08 ．480339−08 ．120375−13 ．208669・11 ．208669−11 4段数

@5次法公式52_公式53 A53_A52 ‘魏42681氾o．．381249−07 ．175706−06．176161−06 ．169460−06．103440−06 ．579567−09．353918−09 ．542005−08．554380−08 ’．6671訂畑9：．371698−08 ．361456−13．212504−13 ．152578−11．178303−11 ．178019−11．125917−11

公式54 A53 ．335452−07 ：i：‘ 泊6．’ ．150705−06 ．517219−09 ’1：64’・．389023−08 ．326378−13 ．47 ・1：1：： ’；層：？6 層そ2． 2．初期値問題理論解 y，＝y＋Ain（2x）， sin（2x）＋2cos（2x）

y−

@ 5 y（0）＝−0．4

（∂f

_{ﾝy〉・）} ステップ数100 選出 _{誤差（刻み幅h＝O．1）} 誤差（刻み幅h＝0．05）誤差（刻み幅h＝0．01）公式名 _基準

1st ste Hax． 1ast 1st ste Hax． 1ast 1st ste Hax． 1ast

Gruttke5 A51 ．172207−07 ．518161−02 ．518161−02 ：．：；2056’二〇9’ ．109479−05 ．109479−05 D98449−14．．667518−11 ．667518−11 公式S1 A52 ．195803−07 ．474489−02 ．474489−02 ．253318−09 ．100503−05 ．100503−05 ．134422−13 ．634556−11 ．634556−11 4段数 @5次法公式52 A53 ’i：：’ P鰯8じo？．159713−02 ．159713−02 ．296667−09 ．328486−06 ．328486−06 ．205168−13 ”湖1鋤・i2．．871310−12

公式53 A52 ．330662−07 ．529130・03 氾29130ro3；．539781−09 層纐2識91iそ｛輝．．ia23鋤i≒頒9 ．356351−13 ．241408−11 ．241407−11

公式54 A53 ．309146−07 ．113915−02 ．113915’02 ．496478−09 ．251865−06 ．251865−06 ．326378−13 ．114756−11 ：：：．マ68598〔1 3．初期値問題理論解 y，＝y＊cos （x）， y （O）＝1．0 ＿ sin《x） y − e （一

恊U動）

ステップ数100 選出誤差（刻み幅h＝0．1）誤差（刻み幅h＝0．05）誤差（刻み幅h＝0．01）公式名

敦定基 1st ste Hax． 1ast 1st ste Hax． 1ast 1st ste Hax． 1ast

Gruttke5 A51 B755085二田i：：；：：824238’；08i’”：：9_{Q囎22泰09：：} _{．551479−10}

．824238−08 ．249221−09 ．567948−14 ．253372−11 ．253372−11 公式51 A52 ．400565−08 ．236757−06 ．464920−07 _潟7期・1Q 、661413−08 ．1雛鋤一〇9 抽田2458剥6・．214156−11 ．214156−11 4段数 @5次法公式52_公式53 A53_A52 ．700571−08．314666−07 ．220706−06．304689−06 ．693297−07．111067−06 ．505363−09．147972−09 ．537962−08．618044−08 ．232539−09．274670−09 ．324185−13．112704−13 ．146516・11．176754−11 ．176754−11．362501−12 ハ式54 A53 ．280840−07 ．217443−06 ．867581−07 ．446757−09 ・：41 ・．151354−09 ．284321−13 ’：：湿06 ；11：・：㍉U ．・3’ x2y2 4．初期値問題 y’＝− 3 ， 9_理論解

y＝

@x3＋1 ・ y（2）＝1．0 ∠ （非線形）ステップ数 50 選出 @基誤差（刻み幅h＝0．1）誤差（刻み幅h＝0．05）誤差（刻み幅hニ0． 01）公式名

1st ste Hax． 1ast 1st ste Hax． 1ast 1st ste Hax． 1ast

Gruttke5 A51 ．630314−06 ．100201−05 ．822456−08 ．112627−07 ．329058−07 ．260307−09 ．802269−12 ．109900−10 ．641846−11

公式51 A52 ．543231−06 ．872240−06 ．742017−08 ．101129−07 ．298611−07 ．242393−09 ．737176−12 ．102151−10 ．602435−11 4段数

@5次法公式52_公式53 _A52A53 ；；1．盤1354言06；．．，423053疏．546309−06 ．933835−06 詩正捌鯉8；：，”1554箔1≒櫨．882470−08 ．927011−08 ：i・：ユ531？o≒征．．296420−07 二．1㎜1二〇9・．272884−09 G魏546正．12∫．632020−12 i・：5卿223一控’．963622−11 E）301536三11．607521−11

ハ式54 A53 ．577225−06 ．976320−06 ．904930−08 ．983450−08 ．310363−07 ．280014−09 ．672997−12 ．100985−10 ．630344−11 参考文献評価への応用について 1） Byrne， G◆D．＆ Lambert，R．J．： Pseudo−Runge− Kutta 匝ethods involving two points J． Assoc． Comput． Mach．，13，（1966），114−123） 2） Byrne， G．D．： Parameters for pseudo Runge− Kutta Hethods， Comn． ACH，10，（1967），102−104 3） Gruttke， W．B．： Pseudo Runge−Kutta Methods of the Fifth Order， J． Assoc． Comput． Mach．， 17，（1970），613−628 4）田中正次：Runge−Kutta法の特性について情報処理 Vol．18 NO．9，（1967），870−878 5）田中正次：Pseudo−Runge−Kutta法とその2nd および3rd order Runge−Kuttaの打ち切り誤差の情報処理 Vo1．10 NO．6，（1969），406−417 6） Costabile，F．： Metodi Pseudo−Runge−Kutta di seconda specie， calcolo，7，（1970），305−322 7） Shintani， H．： Two−step methods for ordinary differential equations， Hiroshi皿a Hath．J．，14 （1984），471−478 8） Nakashi皿a，H．： Pseudo−Runge−Kutta Processes， Pub1．RIMS，Kyoto Univ．，23，（1987），583−611 9）三井斌友：数値解析入門，朝倉書店（1985） 10）新海、清水：Pseudo−Runge−Kutta法の特性についてその2 山梨大学計算機科学科卒業論文（1988）

Byrne型Pseudo-Runge-Kutta法の特性について 利用統計を見る