［山梨大学工学部研究報告第42号1991年12月］
論 文
6段数5次陽的Runge−Kutta法の安定性について
山下茂 矢崎寛 田中正次
（平成3年8月31日受理）
On the Stability of the Six-stage Fifth-order Explicit
Runge-Kutta Method
SigeruYAMASHITA HiroshiYAZAKI MasatuguTANAKA Abstract In this paper we propose three six−stage fifth−order explicit Runge−Kutta methods which have excellent stbility in a sense． 2．安定多項式のグラフ
1．はじめに
Runge−Kutta法は，常微分方程式の初期値問題 y’＝f（x，y） y（x。）＝y。 （1．1） を解く代表的な数値解法である。ここでfは十分滑ら かとする。特に6段数5次陽的Runge−Kutta法にお いては，種々の特徴をもつ数多くの公式が提案されて いる。 6段数5次陽的Runge−Kutta法の一般式は次のよ うに表すことができる。 f診1＝hnf（Xn，ヱ！n） （1．2） i 1 為一力♂（エn＋aihn，Yn＋Σ6品）（i−2，3，…6） 」＝1 6 Y。＋、＝9n＋Σ Cifei hn−Xn＋一工n iニ1 ここで，ai， b，， Ciは公式を特徴づける係数で，公式 （1．2）が5次の精度を持つように選ばれる。 我々は陽的Runge−Kutta法の3つの特性（打ち切 り誤差・丸め誤差・安定性）のうち特に安定性に着目 して，夙（ここでhは刻み幅，λはテス・ト方程式y・＝ λy（1．3）の右辺の複素定数である。）を実軸上で変化さ せたときの誤差伝播率（安定多項式の値）の挙動につ いて調査する。ついで絶対安定区間の長さ・誤差の拡 大率の二つについて優れている公式を提案し，同次の 既知公式と数値例や特性値を通して比較し，得られた 公式の有効性を示す。 ＊電子情報工学科，Department of Electrical Engineering and Computer Science． 安定性の諸概念についての説明は文献（1）に譲る。 6段数5次陽的Runge−Kutta法の安定多項式は P（hλ，γ・）一ゑ（努！）h＋76（響 （2．1） と表される。ここで76は安定性を支配するパラメータ で，公式の係数のみの関数である。γ6を既知公式のそれ γ6＝100 Shanks Lawson γ6 ＝10P
（hλ） 10 5 1 一5 0 5 一1 γ6 ＝−P0 γ6 ＝−P0 γ6＝−100 図一2．1各公式の安定多項式のグラフ hλ 一27一平成3年12月 山梨大学工学部研究報告 第42号 を含む様々な値に固定したときの，誤差伝播率（安定多 項式の値）のグラフを図一2．1に示す。この図から， hA≦0の領域で様々な特徴を持つようにγ6を定めて も，γ6が正の値である限り夙＞0の領域ではほぼ同じ ような挙動を示し，正の方向に発散してしまうことが わかる。また，76が負の値のときはこれとは逆にhA＞ 0，hA≦0どちらの領域でも負の方向に発散するが， γ6（hλ）6／6！が他の項を打ち消すほど小さくなるま でには76の値により多少の差がある。しかしこの場合 は，絶対安定区間が非常に短くなって実用にたえない。 3．今回提案する公式 今回は，文献（1）の考察をもとにして，絶対安定区間 ができるだけ長い公式，集積誤差の伝播率が小さい公 式，これらの中間的な性質を持つ公式の3種の公式の 導出を試みた。 絶対安定区間が最長になるのは図一3．1から，安定多 項式のグラフがy＝−1．0に接するときであることが わかるが，この場合絶対安定区間内にlP（hA，γ6）1＝1．0 となる点が存在し，そのような刻み幅に対しては誤差 が縮小せず，区間内においてどこでも安定であるとは 言えない。そこで次の公式HYO9を提案する（図一3．1 参照）。公式HYO9は， hAが負の実軸上で変化すると き，安定多項式のグラフがy＝0．9に接するように作ら れたもので，絶対安定区間がほぼ最長に近く，かつ安 定区間内に1P（ hλ r6）1＝1．0となる点が存在しないとい う長所をもっている。 また公式HYOOは，安定多項式のグラフh： y＝0．0に 接するように作られたもので，図一3．1から絶対安定区 間内の大部分において，他の公式に比べて安定多項式 の絶対値が小さいことがわかる。 最後に公式HYO5は，この2つの公式の中間的な性 質を持つものとして，安定多項式のグラフが〃＝−0．5 に接するようにして作られたものである。 【新規に提案する公式】 以下に示す3公式の係数は，4倍精度で計算して有 効桁16桁をとったものである。 《
HYO9》 （γ6：0．5268502962
0。0000000000000000D十〇〇 〇．1429260000000000D十〇〇 〇．7655560000000000D十〇〇 〇．2883760000000000D十〇〇 〇．8144880000000000D十〇〇 〇．1000000000000000D十〇1 0．1429260000000000D十〇〇 al a2 a3 a4 a5 a6 b21 b3ユ b32 ： b4ユ ニ b42 ： b43 ニ b51 ： 一〇．1284721728111050D十〇1 0．2050277728111050D十〇1 0．9572989505952828D−01 0．1700866049271763D十〇〇 〇．2255950001329544D−01 −0．4305347706123246D十〇1 K−T．1 HYO9 冊05 @ Lawson HYOO Butcher lShanks @ Fehユberg P（h λ， hλ 一 一4 一2 0 1 一1 ） 図一3．1 今回提案する公式と既知公式の安定多項式のグラフ 一28一6段数5次陽的Runge−Kutta法の安定性にっいて 《 《 b52 ニ
HYOO
al a2 a3 a4 a5 a6 b21 b31 b32 b41 b42 b43 b51 b52 0．7484269993452506D十〇1》 （γ6＝0．6495536165
0．0000000000000000D＋00
0．1006400000000000D十〇〇 〇．8055560000000000D十〇〇 〇．2982890000000000D十〇〇 〇．8527890000000000D十〇〇 〇．1000000000000000D十〇1 0．1006400000000000D十〇〇 −0．2418412944435612D十〇1 0．3223968944435612D十〇10．3062668016336113D−01
0．2427648456806534D十〇〇 〇．2489747415598549D−01 ： −0．5493826523510360D十〇1 0．7662584546262387D十〇1HYO5》 （γ6ニO．5510532079
al a2 a3 a4 a5 a6 b21 b3ユ b32 b41 b42 b43 b51 b52 b53 ： b54 ＝ b61 ニ b62 ： b63 ： b64 ＝ b65 ニ Cl C2 C3 C4 C5 C6 0．0000000000000000D十〇〇 〇．1348760000000000D十〇〇 〇．7655560000000000D十〇〇 〇．2885560000000000D十〇〇 〇．8221740000000000D十〇〇 〇．1000000000000000D十〇1 0．1348760000000000D十〇〇 −0．1407091428512115D十〇1 0．2172647428512115D十〇10．9024662689686832D−01
0．1747076930730974D十〇〇 〇．2360168003003427D−Ol −0．5357500814644698D十〇1 0．8993550993195941D十〇1 一〇．1176816295094726D十〇〇 −0．2246752657819788D十〇1 −0．2864408456991706D十〇1 0．6430679140120729D十〇1 0．1291503116724074D十〇1 −0．3296410764033206D十〇1 −0．5613630358198898D十〇〇 〇．8755923256701552D−01 0．0000000000000000D十〇〇 〇．6583239395491788D十〇〇 〇．4319226575989059D十〇〇 −0．2655669689391899D十〇〇 〇．8776113922408970D−Ol ） ） b53 b54 ニ b61 ： b62 ＝ b63 ＝ b64 ニ b65 ニ Cl C2 C3 C4 C5 C6 ニ ー0．2956802955379996D−Ol −0．1286400993198226D十〇1 −0．7903344524844271D十〇10．1225780820463179D十〇2
0．1368864872866470D十〇1
−0．3952532919473902D十〇1 −0．7707956331800891D十〇〇〇．9062910344267865D−01
0．0000000000000000D十〇〇 〇．9397691251710854D十〇〇 〇．4466972160935955D十〇〇 −0．5897244188253473D十〇〇 〇．1126289741179877D十〇〇 b53 ニ b54 ： b6ユ ＝ b62 ニ b63 ： b64 ニ b65 ＝ Cユ C2 C3 C4 C5 C6 一〇．1322186290194199D十〇〇 −0．2681657549531823D十〇1 −0．3571697981951037D十〇1 0．7388973042982103D十〇1 0．1181831766169313D十〇1 −0．3535503166722029D十〇1 −0．4636036604783496D十〇〇 〇．8758317113842049D−01 0．0000000000000000D十〇〇 〇．6231033927104228D十〇〇 〇．4323598682654962D十〇〇 −0．2320575003503454D十〇〇 〇．8901106823600596D−Ol 4．公式の特性 既知公式及び今回提案する3公式の各特性値を表一 4．1に示す。各特性値の正確な定義については，たとえ ば文献（4）を参照されたい。 読者の便宜のために簡単な定1生的な説明を付け加え よう。表一4．1の見出しの中で，A，2及び．4，3は打ち切り 誤差の大小を判定するための尺度で，その有効1生につ いて十分な実績があるものである。A（Se）びαは安定 性（誤差伝播特1生）の優劣を評価するための尺度で， 前者は絶対安定領域の面積，後者は絶対安定区間の長 さを表わし，これらの数値が大きいほど安定性がよい。 最後にRは公式の丸め誤差特性を判定するための尺度 で，公式の各kiのg部の係数及びY。＋1の係数の絶対値 の和である。もともと，公式の係数の最適化の際に丸 め誤差特性が著しく劣化することがあるので，それを 防止するために導入された数量である。Rは小さい程 よいことは勿論であるが，多少の差は問題にならない。 一29一平成3年12月 山梨大学工学部研究報告 第42号 5次法 、452 、453 R A（s．） α Butcherl ．Q．343d−03 7．262d−07 10．5 17．05 3．39 Lawson 4．861d−03 2．192d−06 7．8 25．84 5．60 Luther 1 1．597d−02 3．552d−05 5．4 10．10 2．52 Nystrom 1．049d−02 1．475d−05 12．3 15．17 3．21 Fehlberg i．172d−03 4．477d−07 27．9 15．44 3．19 Shanks 1．870d−04 5．522d−09 12570．3 18．37 3．55 HYO9 2．712d−03 1．098d−06 33．89 24．07 6．18 HYO5 2．606d−03 9．997d−07 38．77 25．76 5．77 HYOO 2．227d−03 6。406d−07 48．94 24．65 4．71 図表一4．1 各公式の特性値 5．数値実験 下記の弱stiff問題をいろいろな方法と刻み幅で解 き，エ＝3．0における誤差を図一5．1に示した。この図に おいて横軸には刻み幅が，縦軸にはx＝3．0における誤 差の絶対値の常用対数がとられる。 霜 舗一2
