陽的Runge-Kutta法の評価ルーチンとその特性について(常微分方程式系の数値解析とその周辺)

陽的

Runge-Kutta

法の評価ルーチンとその特性について

山梨大学 工学部 矢崎 寛 (Hiroshi Yasaki)

日本大学 工学部 田中正次 (Masatsugu Tanaka)

山梨大学 工学部 山下 茂 (Shigeru Yamashita)

\S

1.

まえがき

Runge-Ku$t$ta 法は, C.Runge,W.

Kutta

らにより 1900 年ごろ創始された常微分方程式の初期値

問題 $y$ $=$ $f$ ($x$ , y) $y(x_{0})=y_{0}$ (1.1) を解く代表的な数値解法で, 1 段階法で最も普通に使われる方法であり, 他の多段階法の出 発値の計算をはじめとして, その安定性とアルゴリズムの簡単さのために長区間を連続的に

積分する場合にも科学者や技術者によって古くから愛好されてきた方法である.

ここで陽的

$Runge-Kutt$

a

法の一般形を示そう. 常微分方程式の初期値問題(1. 1) において, $X=x$

.

における $y$ の数値解 $y$

.

が得られているとき, $x$ ニ $X$ 。$+$ における $y$ の数値解 $y,+$ を $y_{n+}=y_{\mathfrak{n}}+\sum_{1=1}c_{1}k_{I}$ $(n=0,1,2, \ldots)$ $k_{1}=h_{\mathfrak{n}}f(x_{\mathfrak{n}}y_{\mathfrak{n}})$ $k_{I}=h$。$f$

$( x_{\mathfrak{n}}+a_{1}h_{\mathfrak{n}}, y_{n}+\sum_{\text{」}=1}^{I-1}b_{I}, k_{\text{」}})$ _{$(i=2,3, \ldots p)$}

$h_{\mathfrak{n}}=x_{\mathfrak{n}+1}-x_{\mathfrak{n}}$ (1. 2)

によって求める方法を$p$段数陽的

$Runge-Kutt$

a

法という. $a_{I},$$b_{Ij}$ ,c、は公式を特

徴づける係数であり

$a_{1}=$ 」

$1 \sum_{3}^{1- 1}b_{I}$ $(i=2,3, \ldots gp)$

とする. また $a_{1},$$b_{I}$ $c_{1}$ は, 次のような係数マトリックスによっても表現される. $p$ $q$次法を得るための最低必要段数$p$は, $q\leqq 4$ では _$p=q$ であるが, $q>4$ では $p>q$ である ことが分かっている. 表 1.1 に段数と到達可能次数の関係を示す. 表1.

1

段数と到達可能次数の関係 また, 以下の年表は

$Runge- Kutt$

a

法の歴史を簡略に表記したものである。

こあような歴史的背景を経て, 現在では低次から高次に至るおびただしい数の陽的 Runge-Kutta法が存在し,

いかなる公式を利用すべきかユーザを戸惑わせるような現状を招来するに

至っており, 現在公にされている多くの公式の特性の総合的な調査が強く求められている

.

そこで我々の研究室では, 石原, 高山, 渋谷, 大坪により昭和

58

年度から

60

年度にかけ

て, 陽的 Runge-Kutta公式の段数と係数al.$b_{Ij}$

.

_Clが入力されたとき, その公式について次数,

打ち切り精度, 丸め誤差特性, 安定特性などの重要な情報を提供する陽的 Runge-Kutta 法 診断ルーチン D R E RK$F$ (Diagnostic Rout ine for Explicit Runge-Kutta Formula) を開発した (文献 (6), (11))

.

このルーチンでは打ち切り精度判定基準を, 局所打ち切り誤 差の主誤差項 $h^{q+1}t_{q}$ _{に関して求めているが, 今回われわれは}DRERK$F$に改良を加え, 更に $h^{q+2}tn+\iota$ _{の項についても求め, より慎重を期すようにした}. _また, 新たに公式の単調 性についても診断できるようにするとともに, やや物足りなかった絶対安定領域の図形情報 の出力機能についても拡張を加え, 絶対安定領域の内外に等高線を出力できるようにし, さ らに同一画面上に安定多項式のグラフも出力させ, 両者の比較が容易にできるようにするな ど安定性のより正確な把握を期した. 本研究は, このようにして発展させた新しい診断ルーチン

DRERKF

を用いて, おび ただしい数の陽的 Runge-Kut

ta

公式それぞれにっいそ特性の調査を行い, それを総合的にま とめることによってユーザの意志決定を支援し, さらには多くの研究者に有用なデータを提 供するという希望のもとに進められた. \S 2. 陽的 Runge-Kutta 法の評価の観点と尺度

2.

1

打ち切り誤差の大小とその評価 p段数q次法の $x=x_{\mathfrak{n}+1}=x_{n}+h$ における近似解 yn+’ の局所打ち切り誤差は $T_{\mathfrak{n}+1}=g_{q}h^{q+1}+O(h^{q+2})$ (2.1.1) と表すことができる. ここで $g_{q}$ は, 次式に示すように公式の係数のみの関数 $a_{q}$

,

と, (1. 1)

の右辺の関数に依存

1

する関数ベクトル

$g_{l}(x_{\mathfrak{n}}, y_{\mathfrak{n}})$ との積和の形で表現される. $g_{q}=\sum_{J<1}a_{\mathfrak{n}J}g_{1}(x_{\mathfrak{n}}. y_{n})$ (2.

1.

2) ここでlq は, 微分方程式(1.1) が単一か連立かと, 公式の次数 q に依存して定まる自然数で ある 我々は, p段数q次法の打ち切り誤差の大小を判定するために, その有効性につい て十分な実績をもっ次の二っの数量$A_{q2}$及び$A_{q3}$を使用する.

Aq2 $= \sum_{/=1}^{1q}|a_{n}\downarrow|$ (2.

1.

3) $A_{03}=\sum_{\approx 1}^{1q}a_{q}\downarrow 2$ (2.

1.

4)

2.

2

安定性とその優劣の評価

一般の常微分方程式(1. 1)を線形微分方程式

$y’=\lambda y$ ($\lambda$ は複素定数) (2.

2.

1)

で近似しよう. いま $P$段数$q$次陽的

$Runge-Kutt$

a

法 (以下$p$段数$q$次法と略称す

る) (1. 2)を (2.

2.

1)に適用すると, 簡単な計算から

$e_{\mathfrak{n}+1}=P_{p.q}$(hh$\lambda$, _{$7_{n+1}\ldots..\gamma_{p}$})$e_{n}+E_{n}$ (2.

2.

2)

が得られる. ここ て$=y_{n}.y(x_{\mathfrak{n}})$ はそれぞれ $x=x_{\mathfrak{n}}$ における数値解及び理論解,

$e_{n}=y_{n}-y(x_{n})$ (2.2.3)

$P_{n.u}$$( h\lambda, \gamma_{\pi+\iota}\ldots., 7)=\sum_{k=0}^{q}(\underline{h_{k}\lambda_{!})}^{k}+$

$\sum_{1=n+1}7I$ ( $\underline{h_{i}\lambda_{!})}^{1}$ (2.

2.

4) また, $E_{n}$ はこのステップで発生する局所誤差である. $P_{p.q}$$(h\lambda, 7_{q+1}\ldots..7_{0})$は$P$段数$q$次法の安定多項式と呼ばれる. (2.

2.

4)式の各 71 は, 公式の係数のみの関数である. (2.2.2)式から, 安定多項式は集積誤差の次のステップへの伝 播率になっていることがわかる. $p$段数$q$次法(1. 2) は, $|P,$_$.q$$(h\lambda, 7_{q+1}\ldots..\gamma,)|\leqq 1$ (2.

2.

5) であるとき, その h\mbox{\boldmath$\lambda$}に対して絶対安定であるという. また次式によって定義される領域 S を $P$段数$q$次法の絶対安定領域という. $S=$

{

$h\lambda|$ $|P_{p.q}$($h\lambda$ , $7_{B}+1,$$\ldots$, 7p) $|\leqq 1,$ $h\lambda$は復素数} (2.

2.

6) 安定性に有効な原点を含む絶対安定領域の単連結な閉部分領域$Se$を, $P$段数$q$次法の有効 絶対安定領域と呼ぶことにしよう. そのとき次式によって定義される閉区間$S1$_を, $p$段数$q$ 次法の絶対安定区間という.

$S1=R\cap Se=$ $[-\alpha, 0.0]$ (2.

_2.

7)

ただしR は実数全体の集合である. ここでは有界な絶対安定領域をもつ

p

段数q 次法の安定 性の優劣を, (1)有効絶対安定領域$Se$ _の面積

A

$(Se)$ _の大小 (2)絶対安定区間の長さ$\alpha$

(3) 有効絶対安定領域の包含関係 (4)同じ$h\lambda$ に対する安定多項式の絶対値め大小によって評 価する.

2. 3

丸め誤差に関する性質の評価とその他の観点 公式の丸め誤差に関する性質の良否を判定するために

,

やや粗雑な尺度ではあるが, 次式 によって定義される数量

Rl

を使用する. Rl $= \sum_{1=1}^{p}|c1|+\sum_{I=2}^{p}\sum_{\downarrow=1}^{I- 1}|b1j|$ (2.

3.

1)

また, すべての係数が非負で, $0\leqq a1\leqq 1.0$ および

a

$I\leqq a$

,

$(i\langle j)$ が満足されると

き, 陽的Runge-Kutta法(1.2)は単調であるという.

公式の係数の簡単さは (1) 零係数の数, $-(2)$ _{係数が分数であるときの最大分母の桁数,}

(3) _{無理係数であるか否かのチェック, によって判断する}. \S 3. 陽的Runge-Kutta法の評価のためのプログラムDRERK$F$

上のような観点に基づいて, 陽的

$Runge- Kutt$

$a$法の諸特性を評価するシステム

がDRERK$F$_であり, 以下はその診断機能の概要である. 公式名 段数 公式のパラメータ (有理数\rangle 無理数、分数、小数) $c_{\iota}=1arrow 1\Sigma^{1-}b_{I}1$ $1=\Sigma^{p}a_{t}1=1$ 残差による判定 (残差く$10^{-\epsilon}$ ならば $=$ ) $0$係数の数 分数の最大分母の桁数

1

次から

9

次までの次数条件式に よる次数のチェック 到達可能次数が得られているかの チェック 打ち切り誤差の主項の出力 打ち切り精度判定基準 $A_{q2}$

.

A.

3の出力 打ち切り精度判定基準

A

$(q+1)2$

.

A

$(n+1)3$の出力 丸め誤差特性判定基準 Rの出力 安定多項式の出力 有効絶対安定領域の面積の出力 絶対安定区間の出力, 絶対安定領域のプロット図 DRERK$F$は主に

INPU

$T$, $ORDER$ , $ERR0R$_{, $ST$}

ABILIT

$Y$ _{の 4} っ

のサブルーチンから構成されている.

3.

1

SUBROUTINE

I N P U$T$ 機能 診断に必要なデータ, すなわち公式の段数, 公式の係数の種類 (小数か, 分数か, 無理数 を使用しているか) , 公式の係数a1, $btJ$, $c1$ を読み込む. そして係数マトリックスの 出力, さらに係数のチェックを行い, その結果も同時に出力する

.

もしこの時点でエラーが 検出されると, メッセージを出力し診断を停止する. また, 零係数の数と分数係数の最大桁

数も出力する. 対応できる公式の係数のタイプは次の 4 タイプである.

.

小数型 (無理数未使用) 分数型 (無理数未使用)

.

小数型 (無理数使用)

.

分数型 (無理数使用)

3.

$\frac{2SUBR0UTlNEORDER}{\text{機能}}$ 入力された公式が何次法で, 連立微分方程式と, 単一微分方程式のどちらを対象に作られ た公式であるかを次数条件式より判定する. ここで次数条件式は

,

打ち切り誤差および打ち 切り精度判定基準の計算の準備として, 達成次数より

2

次上の次数条件式群の左辺の値まで すべて計算しておく. また, 異常公式を指摘する. すなわち, 入力段数の達成可能次数より 小さい次数法の指摘を行う.

.

アルゴリズム

次数条件式の血汲の正確な値は

,

DATA文を用いて1次から9次までの連立微分方程式, 単一微分方程式にっいてそれぞれ R$R$ (I) $(1=1.483)$

.

$TR$ (J) $(j=1.104)$ に入れてお く. ただしT$R$は連立と重複していないもののみとする. 入力された係数パラメータにより,

次数条件式群の 7i の値

R $L$, $TL$ を計算する. 次 数条件式群が成立しているかどうかの判定基準は, (1) 連立については $|$

RL

$(I)-RR(1)|<1.0\cross 10^{-8}$ (2) 単一については $|$

TL

$(J)-TR$$(J)$ $|<1.0\cross 10^{-\epsilon}$ とする.

3.

$\frac{3SUBR0UTINEERROR}{\text{機}6\text{能}g}$ このサブルーチンではサブルーチン

ORDE

$R$

において計算しておいた達成次数より 2 次

上までの次数条件式群の左辺の値により打ち切り誤差の主項

$g_{n}$ とその次の項 $g_{n+1}$ , それ による打ち切り誤差特性判定基準 Aq2.

A

$(q+1)2$ , $A_{q3}$ ,

A

_{$(q+1)$ $} ($q$ は達成次数) を求め る. さらに丸め誤差特性判定基準と, 単調性についての評価を行う

.

アルゴリズム 連立, 単一の微分方程式の右辺の関数に依存する部分を, それぞれRF(I)(I$=1,482$),$TF(J)$ $(J=1.200)$に, それぞれしまっておく. また, 各次数に対する打ち切り誤差項の項数, 連立の 場合と重ならない単一のものの項数, さらにどの次数条件式と対応しているかなどがわかる ように先頭番号などをにしまっておく. サブルーチンORDE $R$から引き渡されたD $E$により, 連立, 単一のどちらに対応するか 判断し, 連立と単一のいずれかについて打ち切り誤差の主項の値を計算する

.

こうして求め た打ち切り誤差の主項の $a_{QI}$ より打ち切り精度判定基準 Aq2, $A_{q3}$ を計算する.

A

$rz=\sum_{I=}^{nu}]$

a

$q1|$ A$q3=\sum_{1=1}^{num}$(a

$n1$)

$2$ (3.

3.

1) (num ; 次数条件式の個数) また,

.

その次\emptyset項 $g$$,+1$ についても同様}\llcorner \check求める.

っいで丸め誤差特性判定基準Rlを計算し, 単調性の検証もする.

Rl

$= \sum_{1=1}p|c1|$ $+ \sum_{1=2}\sum_{\downarrow=1}^{I-1}\Phi|b1j|$ (3.

3.

2)

3.

4

SUBROUTINE

S

TA

B I L I T$Y$ 機能 安定多項式を出力し安定性を支配するパラメータ 71 を求める. 有効絶対安定領域の面積と絶対安定区間を求める.

$h_{\mathfrak{n}}\lambda$を実軸に沿って動かしたときの安定多項式 $P_{p}$

.

$q(h_{n}\lambda, 7_{q+1}. , , 7_{\Phi})$ のグラフを出力 し絶対安定領域をプロットし, その内外に等高線も画く. アルゴリズム 1 段数 1 次から 11 段数 8 次までの安定多項式を

PP

(]), $(1=1,11)$ にしまっておく. サブ ルーチンERRO $R$から引き渡された段数$P$に応じて安定多項式を出力する. 安定多項式の各項の係数を, 段数$P$, 次数$Q$に応じて計算し, $P$次の項の係数を

ALPHA

(1) に,

$(P-1)$

次の項を

ALPHA

(2) に.

.

.

.

.

1 次の項を

ALPHA

(P) _{に入れる.} ただし

ALPHA

(1)$=0$ であるときは, 後で $D-K-A$法が使えないので, 添字を1ずっずらしておく

(

ALPHA

(1) $arrow ALPHA(2)arrow$ ).

このようにして安定多項式を求めたのち, $h_{\mathfrak{n}}\lambda$を実軸に沿って動かしながらそのときの安 定多項式の値をプロットする.

$|P_{p}$

.

$Q(h_{n}\lambda, 7_{q+l}. , . 70)|=1)$ とおいた方程式を解く. すなわち$\theta$ を$0$ から2 $\pi$ まで $(\pi/45)$ ラジアンずっ変化させ, 毎回この$p$次方程式を解く. これにより得られた $(90*P)$ _{個の解について, その実数部を REP}$(J, [)$ に, 虚数部を IMP $(J, 1)$ $( J=0,89. 1=1, P )$ にしまう. また,

解をソートして複素平面上にプロットする.

この ソーティングは, 第1番目の解の集合 (90 個) から第$P$番目の解の集合 (90 個) の入れ 替えを行い, 原点から一連の動作で単連結部をプロ ッ トする為に行う. また, ソーティ ング

された有効な解について順次台形則を適用することにより面積を近似する.

面積は, 単連結 部全体, 有効部分, 無効部分に分けて計算し, 同時に絶対安定区間も求める. また絶対安定領域の中と外に等高線を画くため, $|P_{p.q}\langle h_{n}\lambda$

.

_{$7_{q+1},$}

. .

$7_{0}$) $|=0.9.0.8$

0.

1

, $|P_{p}$

.

$q(h_{\mathfrak{n}}\lambda, 7_{n+1}, . . 7,)|=3.0$ _,

5.

$0$ ,

10.

$0$ とおいた方程式を解き, 得 られた解をソートしてプロットする. (i) _{安定性を支配するパラメータ 71 の求め方} 安定性を支配するパラメータ

71

は次のような式から求めることができる

.

71 $=i$! $\lambda b\Downarrow$ $Zh_{4}x$ $(q+1\leqq i\leqq p)$ (3.

4.

1)

$i$ 個の積和

安定多項式(2.

_2.

4) からもわかるように 71 は$(p-q)$_個ある.

(1 I)

_{絶対安定領域の面積の有効無効の判定基準}

【定義3. 1 】 AREA, AREAI, AREA2,AREA3, (AREA4, AREA5)の定義

AREA

;

絶対安定領域のうち原点を含む単連結の閉領域の面積.

AREAI ;AREA の中で, 第一, 第四象限にはみ出した部分の面積. 無効とみなす. AREA2 ;AREAの中で有効とみなされる部分の面積.

AREA3

; 絶対安定領域のうち原点を含む単連結の部分にくびれがあった場合

,

その くびれより左側にある部分の面積. 有効な場合と無効な場合がある

.

AREA4,AREA5は2 っめ,

3

っめのくびれがあったときのためのものである. 面積の有効・無効の判断は, 図3.5 のようにくびれがあった場合, そのくびれの部分の Imaginary Part の大きさを LMIN, 第2象限の Imaginary

Part

の最大のものを

LMAX

とす るとき, $\frac{LMIN}{LMAX}<\frac{1}{10}$ _{ならば無効とみなし,} $\geqq\frac{1}{10}$ ならば有効とみなす. これは図3.

6

のようにくびれが大きいとき, $h_{n}\lambda$を図のようにとれば, 途中に不安定な部分が出てくる

からである.

AREAI

を無効とするのもこの理由からである. AREA3, AREA4,

AREA5

は, このくび れの場合の対応のためのもので, くびれが 3 つある場合まで対処できる. 原点 Ｉ坩堕蠅壁

3. 5

評価ルーチンDRERK$F$の使用法 まず, 公式を評価したいユ$-$_ザは, 次のような公式の係数のファイル

9

_を作る

.

仝 式の係数が分数で, 無理数を使用していない場合 (例.

Heun

の3段数3次法) ファイルの内容

:

説明

3:

公式の段数 $F$

:

係数の種類 (分数は traction) $N$

:

無理数を使っているか (地)

1.

$OdO$

.

:a2 の分子

3.

$OdO$ :a2 の分母

2.

$OdO$

.

:as

の分子

Heun

の3段数3次法の

3.

$OdO$

:aa

の分母 係数マトリックス

1.

$OdO$

.

:

$b_{21}$の分子

3.

$OdO$

:

$b_{21}$の分母

0.

$OdO$

.

:

$b_{31}$の分子

3.

$OdO$

:

$b_{31}$の分母

2.

$OdO$

.

:

$b_{32}$の分子

3.

$OdO$

:

$b_{32}$の分母 以下略

式の係数が分数で, 無理数を使用している場合 (例.

Gill

の4段数4次法) ファイルの内容

:

説明

4

:

公式の段数 $F$

:

係数の種類 (分数はひaction) $Y$

:

無理数を使っているか $(X^{es})$

1.

$dO$_, :a2 _の $d$ の部分 (右図参照)

0.

$dO$

.

:a2 _の $e$ の部分 (右図参照)

$\frac{d+e\sqrt{f}}{g}$

0.

$dO$

.

_{:a2 の} $f$ の部分 (右図参照) 公式の7tーマット

2.

$dO$ :a2 _の $g$ の部分 (右図参照)

1. dO.

:

$a_{3}$ の $d$ の部分

0.

$dO$

.

:

$a_{s}$ の $e$ の部分

0.

$dO$

.

:

$a_{3}$ の $f$ の部分

2.

$dO$

:

$a_{3}$ の $g$ の部分

$l^{\backslash }’\lambda$下$\Re$

８ 式の係数が小数で, 無理数を使用していない場合 (例. 田中の5段数4次法) ファイルの内容

:

説明

5:

公式の段数 $D$

:

係数の種類 (小数は _{$\lambda^{ecima1)}$} $N$

:

無理数を使っているか (蜘)

0.

$235d0$

:

_a2

0.

$44d0$

:

$a_{3}$

0.

$994d0$

:

$a_{4}$

1.

$0d0$

:

$a_{5}$

0.

$235d0$

:

b21

$-0.02727517047d0$

:

$b_{S1}$

0.

$.4672751705d0$

:

$b_{ 2}$ 以下略 じ 式の係数が小数で,$\cdot$ 無理数を使用している既知公式は今の所ないので省く. 以上のような公式の係数ファイル

PAR. IN

を作ったら, 後はDRERK$F$を実行させれば 良い. プログラム中で

OPEN

($5$

.

$FILE=$

PAR. I

N) となっているのですぐに公式の評価結果

が得られる.

3. 6

使用例 この節では使用例として4段数4次古典的 Runge-Kutta 法の評価結果を掲載する. $*************************THE$

DIAGNOSIS

$****************************$ –

I NPUT DATA

– 公式の入力のチェック

$*stages–4$

:

段数の出力

$*4-stag_{-}e$ explici$t$ Runge-Kutta

formula

is

:

4 段数陽的 R-K 法の一般形

Kl

$=$

Hn

$*F$(

Xn

, Yn)

$Ki=$

Hn

$*F$(

Xn

_$+Ai*Hn$ ,

Yn

$+Sum$( $J=1$

.

$i-1)$

_Bi

_$i*Ki$ )

$[ i=2. 4]$

$Yn+1=$

Yn

$+Sum$

$(i=1 , 4)$

Ci S Ki $Xn+1=$ _Xn $+Hn$

*formula –

:

係数マトリックスの表示 1/

₂

$|$ 1/

2

1/

2

$|$ 0/

1.

1/

2

1/

₁

$|$ 0/

1

0/

1

1/

1

$|$ 1/

6

2/

6

2/

6

1/

6

*check –

:

係数のチェック

Th$is$ formula $is$ perfect

1

I$t$ has

no

error]

–Coefficients

–

$0$ coeff$ici$

ents

–

3

:

$0$係数の数

The largest

cofficient

–

1 Keta

:

係数の最大桁

–

ORDER

– 次数のチェック

This formula

is 4

th order method

The resul$t$ of test, thi$s$ formula

sati

sfied cond

itional

equat

ions

of 4th order

:

この公式は 4 次の次数条件式

for sys

tems

of D.

E.

を満たしていた

Therefore, this

formula is

4-stage 4th

order

explicit Runge-Kutta

formula

for

:

この公式は

4

段数

4

次陽的

systems of

D.E.

Runge-Kutta法である

$*$ check –

Th$is$ formula $is$

normal

1

For

a

4-stage method, the best possible order

is

41

:4

段数の達成可能次数は

4

次 なのでこの公式はNormalである –

TRUNCATI0N ERROR

– 打ち切り誤差特性 $*$ for systems of

D. E.

$T4=g4*Hn**5+g5*Hn**6+0$$( Hn**7)$

:

打ち切り誤差の主項の表示 $g4=$

a

4 14

$f$ $(f. f, f. f)$

:a4

$1=-0.34722222222222222D-03$ $+a4$ $2*f$ $(f (f. f. ), f)$

:a4

$2=020833333333333333D-02$ $+a4$ $3*f$ $(f(f). f, f)$

:.

a4

$3=-0.20833333333333333D-02$ $+$

a

4 44

$f(f (f, f. f))$

;a4 $4=013888888888888889D-02$ $+$

a

4

5*

$f(f(f (f. f)))$

:a4

$5=-0.20833333333333333D-02$ $+a4$

$6*f(f(f(f(f))))$

:a4

$6=083333333333333333D-02$ $+$

a

4

7* $f$

$(f(f). f(f))$

:a4

$7=-0.62500000000000000D-02$ $+$

a

4 84

$f$

$(f(f(f)). f)$

:a4

$8=-083333333333333333D-02$ $+$

a

4 94

$f(f (f(f), f))$

:a4

$9=0.41666666666666667D-02$

–CRITERIA OF

THE

SIZE

OF

TRUNCATI0N ERROR

(1st term) – $*$ for systems of D.

E.

$*A42=$

Sum

$( i=1, 9)$ abs(

a 4

$i$ )

:

打ち切り精度判定基準A42,$A_{43}$

の定義 $*A43=$ Sum $(i=1, 9)$ $squ$(

a

$4i$ )

A

42

$=$ $0$

.3506

$9444444444444D-01$

:

打ち切り精度判定基準A42

A

43

$=$ $0$

.21038

$290895061728D-03$

:

打ち切り精度判定基準A43

:

打ち切り誤差の第

2

項の表示 $g5=$

a5

$1*f(5)(f. f. f. f. f)$

:

a

$1=-0.17361111111111111D-03$ $+a5$ $2*f$ $(f(f). f, f. f)$

:

a

$2=-0.17361111111111111D-02$ $+a5$ $3*f$ $(f (f, f)$

.

$f$

.

$f$)

:

a

$3=0S6S05555555555556D-03$ $+a5$ $4*f$ $(f (f. f, f). f)$ ;

a

$4=0.17361111111111111D-02$ $+a5$

$5*f(f (f. f, f. f))$

;

a

$5=0520S3333333333333D-03$ $+a5$

$6*f(f(f (f. f. f)))$

;

_a

$6=-0.34722222222222222D-03$ $+a5$

$7*f(f(f(f (f. f))))$

:

a

$7=013SSSSSS8SSS88889D-02$ $+a5$ $\mathfrak{g}*f(f’(f(f(f(f)))))$ ;

a

$8=013S88888888888889D-02$ $+a5$ $9*f$ $(f(f). f\cdot(f), f)$ ;

a

$9=-052083333333333333D-02$

$+a510*f$

$(f\cdot\cdot(r_{:}f), f(f))$ ;

a

$10=-017361111111111111D-02$

$+a511*f$

$(f (f(f). f). f)$

;

a

$11=011754943508222875D-37$

$+a512*f$

(

$f(f(f))$

.

_f.

f) ;

a

$12=-0.69444444444444444D-02$

$+a513*f(f$

$(f(f). f. f)$

}

;

a

$13=0.31250000000000000D-02$

$+a514*f’(f(f(f’(f))). f)$

;

a

$14=0.69444444444444444D-02$

$+a515*f$

$(f(f(f’(f))). f)$

;

a

$15=055555555555555556D-02$

$+a516*f(f(f (f(f). f)))$

;

a

$16=0.41666666666666667D-02$

$+a517*f(f (f(f). f(f)))$

:

a

$17=-0.10416666666666667D-02$

$+a518*f$

$(f(f(f)), f’(f))$

;

a

$18=-069444444444444444D-02$

$+a519*f$

$(f (f ’ (f. f)). f)$ ;

a

$19=-0.34722222222222222D-02$

$+a520*f(f (f (f, f),$

$f$))

:

a

$20=0.34722222222222222D-03$

–CRITERIA

OF

THE S I ZE OF TRUNCATI0N ERROR

(2nd term) —

$*for$ systems of

D. E.

$*A52=$

Sum

$( i=1. 20)$ _abs(

_{a 51}

): 打ち切り精度判定基準$A_{62}$

.

$A_{63}$

の定義 $*A53=$

Sum

$(i=1. 20)$ $squ$(

a

$5i$ )

A

$52=$ $0.88715277777777778D-01$

:

打ち切り精度判定基準

A52

A

$53=$ $0.46751422646604938D-03$

:

打ち切り精度判定基準$A_{53}$

–

CRITERI0N

OF

ROUND-OFF

ERROR

– 丸め誤差特性 $*$ for systems of

D.

E.

$*R1=$

Sum

$( i=1,4)$ _abs(Ci) $+Sum$$( i=2,4)$_Sum$(j=1, i-1)$

_abs

(Bij)

$*R2=$

Sum

$(i=1.4)$ _abs(Ci)

_:

_{丸め誤差特性判定基準}$R_{1}$

.

$R_{2}$の定義

Rl

$=$ $0.30000000000000000D+01$

:

丸め誤差特性判定基準$R_{1}$

R2

$=$ $0.10000000000000000D+01$

:

丸め誤差特性判定基準$R_{z}$

$*This$

formula is monotone.

:

この公式は単調な公式である 安定性に関する性質

–POLYNOMIAL

OF

STABI

L ITY –

PP

( 44)$=1+Z+Z**2/2[+Z**3/3!+Z**4/4$!

:

安定多項式の出力

$Z=h\lambda$ _である

Where, $Z=h*LAMDA$

–THE AREA OF

THE REGION

OF

ABSOLUTE STABILITY

–

:

絶対安定領域の定義

$*AREA$

:

the

area

of the region of absolute stability

$*AREA1$ ; the

area

of

the

non-effective

region of absolute stability (in _the

first

and the forth quadrant)

$*AREA3$; the

area

of

the

non-effective

region

of absolute

stability (smaller than

LMIN

(1) )

:LMI

$N$ は. \langleびれた部分のReal Part.

AREA

$=$

12.70008252277239400

:

絶対安定領域の面積

AREA

1

$=$ $0.50399783329155307$

:

第 1,

第

4

象限にはみ出た部分

$*AREA2=$

12.

19608468948084094

: 有効絶対安定領域の面積

AREA

$3=$ $0$

. 00000000000000000

:

くびれた部分の面積

$\backslash -$ THE

INTARVAL

OF

ABSOLUTE STABILITY

–

:

絶対安定区間の長さ

$(-2.78516567121230508 0.00000000000000000 )$

以上のデータに加えて,

図

3.8

のような絶対安定領域の図と安定多項式のグラフが得られる

.

虚軸

図3.8

絶対安定領域と安定多項式のグラフ

\S 4.

既知公式の特性

われわれは,

低次から高次までのおびただしい数の陽的

$R$

un

$ge-K$ $utt$

a

法の既知公

式を, 評価ルーチン

DRERK

$F$

を用いて評価した.

_{この章では紙面の都合から,}

7

段数

6

次法についての評価結果のみを表にまとめて掲載するとともに

,

各公式間

の比較を試みる

.

_{なお今回掲載できなかった他の次数の公式の評価結果にっいては,}

近い将来まとめて発表する予定である

.

4.

1

7 段数 6 次公式

打ち切り精度判定基準

$*A62=$

Sum

$( i=1, 48)$

abs

(

_a

₆

$i$ ) $*$

A

72

$=$

Sum

$( i=1.115)$

abs

(

a

7

$i$ )

$*$

A

63

$=$

Sum

$( i=1, 48)$

sqr

(

a

6

$i$ ) $*$

A

73

$=$

Sum

$( i=1,115)$

sqr

(

a

7

$i$ )

丸め誤差特性 _安定特性

$*$ Rl _$=Suw$( $I–$I.7) abs(C1) $*$ Su$n$$( i=2.7)$Sum\langle$j=1$

.

$I-1$) abs{$B$Ii) _{* 安定多項式}

$FP(\eta b)=|*7.*X222/l$!$7.t3/3$!\dagger 7.2*\{/t!*X2$s/5!$CA\mbox{\boldmath $\lambda$}NA62z126/0!

.

$7=h\lambda$

$*$ R2 $=$ Sun $(1–1.1)$ abs$(Ci)$

公式$S$han$ks$_は5_{次法なので}, _{安定性を支配する}

パラメータは

76.

77

の2 っある.

各公式の比較

文献から抽出した

7

段数

6

次公式は

1

6個ある. 公式

Bu

$t$che$r-opt$

.

$L$

awson-op

$t$_,

S TA

B.

TRUN.

I Ml

$0$_,

I

M9.

I

M5.

[Ml _は,

_{ゎれわれが過去に提案した公式である.}

係数が簡 単な公式は

Bu

$t$ch$er(3)$

.

(4)

$,$ (5) である.

打ち切り精度判定基準を見ると

,

公式

TRUN.

$But$che_$r-op$$t$

.

$S$

TA

B.

Lawson-op$t$ が良い

値を示している.

丸め誤差特性判定基準では,

_{ここでも公式 TRUN, それに公式}

_Bu

$t$

che

$r(1)$_, (3), (5) が上位にある. しかし $R1$

_{は著しく大きい場合を除けば}

,

さほど大きな影響はない

ので, 公式[Ml

_{までは許容範囲内である.}

_なお,

単調な公式は無い.

安定特性に関しては

,-

有効絶対安定領域の面積では公式

_STAB.

_{絶対安定区間の長さ}

では公式

Lawson

と IM9,

_{誤差の伝播率が小さい}

[Ml_がよい.

すべての評価基準において上位に位置する公式としては,

_{若干打ち切り精度判定}

基準は落ちるが

,

公式IMl, I$M5$_,

I

$M9$_{があげられる.}

..

まとめると,

一般的なユーザには打ち切り精度判定基準と丸め誤差特性判定基準

が良い公式 TRUN,

_{安定性を重視するユーザには公式}

_{[Ml, IM5, [M9 が推薦できる.}

\S

5.

結語 従来

, 陽的

$Runge-Kut$

$t$

a

法の特性を調査するには膨大な計算量を必要とし,

そのためのプログラムの作成にも大変な時間と労力を必要としてきた

.

しかし,

この評

価ルーチン

DRERK

$F$ の開発によってその評価が容易となり

,

様々な特性値をすぐに

得られるようになったことは, 今後の

$Runge-Kutt$

a

法の研究を強力に支援で

きるものと思われる

.

また本研究では,

DRERK

$F$

を用いて低次から高次までの数多くの陽的

$Runge$

$-Kut$

$t$

a

法の公式を次数ごとに分類し,

様々な評価基準から各公式の評価を行った

.

特に打ち切り精度特性と安定特性の部分は貴重なデータが得られたと思われる

.

これに より公式の相対的な比較が可能となり, 陽的

$Runge-Kut$

$t$

a

法のユーザの意志

決定を支援できるであろう.

今回は紙面の都合から

,

7

段数

6

次法についての評価結果

のみを表にまとめて掲載したが

,

掲載できなかった他の次数の公式の評価結果について

は, 近い将来まとめて発表する予定であるので,

そちらを参照していただきたい.

なお, 陽的

$Runge-Kut$

$t$

a

法公式は数多く存在するため,

埋め込み型公式を

はじめ今回評価しきれなかった公式もあり, また

DRERK

$F$

の出力をより分かりやす

くするなど改良の余地もあるので

, 今後の課題としたい.

【参考文献】

(1)

_{J. C. Butcher : The Numerical}

Analysi$s$

of

0rdinary

Differential

Equat ions,

John

Wi

ley

&

Sons

(1986)

(2)

_J.

_{H. Verner}

_:

Explici$t$ Runge-Kutta

Methods wi th Est imates of The

Local

Truncat

ion

Error,

SIAM J. Numer. Anal. Vol.

15,

No.

4

pp.

772

$-790$ (1978)

(3)

_R.

_{F. Churchhouse}

_:

_Handbook

_of

_Applicable

_Mathemat

ics,

_Vol.

III,

merical

Methods,

pp.

321-335

(4)

G.

$J$

.

Cooper

&J.

H.

Verner

:

Some

Explici$t$ Runge-Kutta

Methods of

High Order,

SIAM J.

Numer.

Anal.

Vol.

9,

No.

3,

pp.

389-405

(5) 田中 正次

:

Runge-Kutta 法の打ち切り誤差に関する研究.

.

東京大学に提出

した学位論文 (1972)

(6) 田中 正次

:

Runge-Kutta

法の研究と評価のためのソフトウェア

.

$l^{\backslash }/$ピュ$-$トロール

No.

12

(コロナ社),

pp.

29-35

(1985)

(7)

_{M. K. Jain : Numerical Solution}

_of

_Differential

Equations,

A Halsted

Press

Book,

pp.

32-40

(8)

J. C. Butcher

:

On

Runge-Kut

ta Processes of

high order,

J.

Austral

Math. Soc.

,

Vol.

4,

pp.

179-194

(1964)

(9)

F. Ceschino

&J.

Kuntzmann

:

Numerical Solut ion of

Ini

$ti$

al

Value

Problems, Prentice-Hall,

pp.

37-77

(1966)

(10)

A. Ralston

&P.

Rabinowi tz :

a

Fi

rst Course in Numerical

Analys

_is.

McGRAW-Hall, pp.

208-233

(1965)

(11) 高山 尚文 :7 段数 6 次陽的 Runge-Kuut ta 法の最適化と陽的 Runge-Kut

ta

法の評価 システムにっいて. 山梨大学計算機科学科修士論文 (1987)