集積打ち切り誤差の推定ができる初期値問題の数学ソフトウェア利用統計を見る

(1)

集積打ち切り誤差の推定ができる初期値問題の

数学ソフトウェア

山下茂若林晴彦石川茂田中正次（昭和56年8月31日受理）

On a Mathematical Software with Global Truncation Error

Estimates for Initial Value Problems

ShigeruYAMASHITA HaruhikoWAKABAYASHI ShigeruISHIKAWA MasatsuguTANAKA Abstract Recently， P． Merluzzi and C． Brosilow contrived Runge−kutta algorithms with built−in estimates of the accumlated truncation error． On the basis of the above methods we make a program which can estimate the accum− ulated truncation erros of numerical solutions of ordinary differential equations． Secondly， we explain how to use our program and how to interpret the results of computation．

1．まえがき

局所打ち切り誤差が推定できるうめ込み型Runge− Kutta公式がR．H． Merson（1956）によって最初に提案されてから，この型の公式は極めて有効であるのでいろいろな研究者によって研究され，低次から高次に及ぶさまざまな公式が公にされている1）2）3）4）。これに対し，近年局所打ち切り誤差の推定とほぼ同程度の手数により集積打ち切り誤差が推定できるうめ込み型の公式がP．Merluzziら（1978）によって提案された5）。本論文において著者は，上記のP．Merluzzi らによる公式を用いて，集積打ち切り誤差が推定できる常微分方程式の初期値問題に対する数学ソフトウェアを作り，その使用法，有効性，結果の解釈などについて述べる。 2．集積打ち切り誤差が推定できるうめ込み型公式常微分方程式の初期値問題を yt＝f（x，N）， rv（Xo）＝ツo （2．1）とする◎ ここで，f（x， or），），（x）はm次元関数ベクトル， yo はm次元ベクトルであって，or（X）は十分な滑らかさをもつものとする。そのとき集積打ち切り誤差が推定できるうめ込み型の公式は次のように表わされる。ツo＊＝ツ（Xo）ゆ ki＝hnf（ッ＊n．1＋Σ鋤是プ）ゼ＝1，2，…，μ ゴ＝1 ルツ1＝ツ（Xo）十Σ（β汁rのkz ゴ＝1 ツ。＝ツ＊n−1＋Σゐ品〃＝2，3，… 輌＝1 ぽ εn＝Σア流 π＝1，2，3，… ‘＝1 ツ・・−Nn＋〔（ hn＋i hn）P−・〕e・−n−・，2・… （2．2）ここで，ornはX＝Xnにおける数値解，εnはその集積打ち切り誤差の推定値である。固定刻みのときは Yn＊＝仇公式の係数は， Runge−Kutta法の集積打ち切り誤差の漸近展開が En＝hpBpy（P）（x）十〇（hp＋1）（2．3）

(2)

によって与えられるように決定される。ここでBpは任意に与える誤差定数である。したがって，Bpを大きくとると数値解ornの打ち切り精度は悪くなるが推定誤差の精度が上がる。一方，Bp小さくとれば，数値解Ynの打ち切リ精度はよくなるが，推定誤差が不正確になる。 P．Merluzziらは，主誤差関数がある仮定を満足する場合についてP＝＝1∼4に対してBpをいろいろにかえたいくつかの公式を提案している。本論文で提案する数学ソフトウェアにおいては，それらの中で最も高精度で実用性の高いと考えられる次の公式を使用した。」yo＊＝ツ（XO） kl＝hf（エn， ryn）

k・−hf（Xn−1＋畜捗・＋音り

k・＝・hf（Xn−・＋｛．1−h， yn−i＋−iltE7ki＋音虎・）（2・4）ぽ∫（〔許・ツ・一・＋91，7k，＋謀・＋音ゐ・）緬∫（Xn−・＋h・ツ耐克・一誓ん・＋手虎・）是・一乃∫ト1＋丁み・ツ・−1＋誓ゐ1＋晶虎・

」謬・＋謬克・）

⊇＋÷ゐ1−19911−llk・＋豊＋☆ゐ・＋、易6虎・泥㌫・＋昔虎・一瀞＋皇為・＋針・＋蛋τゐ・・・獅薰P磯8ゐ・一晶ゐ・＋古虎・＋☆ゐ・この方法は古典的な4次法と同精度の公式で，推定誤差が正確になるようにBpを比較的大きく選んだものである5）。この方法はRunge−Kutta法の主誤差関数がある特殊な形をとることを仮定しており，また集積誤差の漸近展開の第1項のみを評価するものである。したがって，うまく処理できない問題があったり，集積誤差の推定精度が十分でないなどの限られた有用性をもつ。しかし，集積打ち切り誤差を簡単なアルゴリズムによっそ評価しようとするこの試みの意義は大きい。 3．集積打ち切り誤差の推定ができる初期値問題の数学ソフトウェア、ここに提案するプログラムは，2節において述べた 6段数4次のうめ込み型公式（2．4）を使用して，常微分方程式の初期値問題の数値解とその集積打ち切り誤差の推定値を求めるものである。

o注意

（1）プPグラムは固定刻みと可変刻みの両方が使用できるようになっているので，どちらか一方を選択する。（2）可変刻みにおける刻み幅調整法このアルゴリズムでは，刻みを，集積打ち切り誤差の推定値のノルムが指定した最大相対許容誤差 Emaxに等しくなるように調整している。集積打ち切り誤差の推定値をε。とする。εnのノルムをllεnl1で表わすと，P次法の時11 e。 1］の値はhP にほぼ比例するので，次の刻み幅Hは次式で計算できる。 H＝h。（Em。。／II S・ H）i／P そのとき（イ）もし1ε，ノッ杣＞2Emaxならば，第n番目の刻みを捨て，hn＝Hとおいて新しい刻み幅で第 n番目のステップを繰り返す。（ロ）もしEm。x＜H en／ッ。 il≦2Em・xならば，第n 番目の刻みを受け入れて，次のステップでは下記のように刻み幅を減らす。 hn＋1＝H ◎もしO． 5Emax＜llεπ／yn H〈Emaxならば，第n 番目の刻みを受け入れて，次のステップはhn＋1 ＝・ hnとおく。 ←）もし1］εn／yn ll≦90 ． SEmaxならば，第n番目の刻みを受け入れて，次のステップでは次のように刻み幅を増す。 hn＋1＝min｛H，2hn｝（3）パラメータ（入力）初期値Xox｛ n−1 （出力）エn＝Xo＋Σh， nニ1，2，…，e k・＝O （入力）初期値Ylo， Y20，…， orNOをY（1）＝ Ylo，_Y（ Y（2）＝ Y20，…， Y（N）＝ or ivoの順に与える。 n−1 （出力）Xn＝Xo＋Σh，における解 k＝O n＝1，2，…e Ft＞（入力）初期値問題におけるfi（i＝1，2，…， N）を計算する副プログラム名。副プログラムの使用法 SUBROUTINE F（X，Y，YP） X⇒（入力）独立変数

(3)

Y⇒（入力） Y（1）＝Yl， Y（2）＝Y2”…Y（N）＝Wなる対応を持つ大きさNの1次元配列。 YP⇒（出力）YP（1）＝fi（x，ッ1， Y2，…ツN），…， YP（N）＝fN（x，ッ1，ッ2，…， w）なる対応をもつ大きさNの1次元配列。 NEgN⇒（入力）連立1階常微分方程式の元数。 XEND⇒（入力）独立変数の最終値。 H⇒（入力）刻み幅。ただし，可変刻みにおいては，刻み幅の初期値になる。 EPS2V⇒（出力）集積打ち切り誤差の推定値。 Eハ4AX⇒（入力）指定する最大許容誤差。 NFt＞（入力） NF』1 固定刻み。 IVF』2 可変刻み。（4）サブルーチン ATERKl and ATERK2 与えられた初期値問題を，ATERK1が固定刻みで， ATERK2が可変刻みで解き，最終ステップまでの結果，つまり各ステップの独立変数の値，数値解，集積打ち切り誤差の推定値をこの順に印刷する。（5）サブルーチンRK4 RK4は， Runge−Kutta法の増分関数ん（i＝ 1，6）を計算して数値解Yと集積打ち切り誤差の推定値EPSNを求めるサブルーチンである。（6）使用例

C

連立1階常微分方程式ツ、t＝ツ2，ツ、（0）＝0 ｛ツ2∫＝一ツ1元ツ2（0）−1 を可変刻みを使って，最大相対許容誤差EMAX ＝10−4として初期点X＝＝O．0から最終点Xend ＝・5．0まで積分する。 10 20 ＊＊EXAMPLE＊＊ REAL Y（2）

NF＝2

NEQN＝2

X・＝0．O Y（1）＝0．O Y（2）＝1．O H＝0．1

XEND＝5．O

EMAX＝1．O E−4 GO TO（10，20）， NF CALL ATERK1（NEQN， X， Y， H， XEND）

ST6P

CALL ATERK2（NEQN， X， Y， H， XEND，

EMAX）

STσP

END

SUBROUTINE F（X， Y， YP） REAL Y（2）， YP（2） YP（1）＝Y（2） YP（2）＝−Y（1）

RETURN

END

c C c C C c C

10

20

★’★★★★★★★★★★★★★★★★★★★★★iSr∧★★★★★★★★★★★★★★★★★★★★★★★★★★★止★★ ★ 去 ★ ACC．UMUしA†ED TRUNCATエON ERROR OF RUNGε＝1くUTTA ．METHOD ★ ★ ★ ★★★★・★★★★★★★★★★t★★★★★★★★★★★★★★★★★★★★★★＊★★★★★★★★★★★★★★★★ ★t EXAMPしε ☆★ Rf三AL／， Y’（2） NFt・；l

NEQN＝2

X＝0．．O Y（．1）＝0●O Y（2）＝1脅e H＝O●1 XEND＝5●O_E凹層 `X＝1●0ε■4

Gn TO（10・20）gNF

cALし ATERK1（NεQN．XtY・HwXEND）

STOP

CAしし ATERK2（NεQN’X，YgH，XEND●εHAX）

STOP

END

SUBROUT工NE F（X●Y・，YP） REAし Y（2）⑱YP（2） YP．（1）＝Y（2）・ YP（2）＝”Y（1）．

RETURN

END

(4)

C c c c C c c ★★★★★★★★★★★★★★★★★★★★★ ★ ★ ★ SU6ROUT工NE ATERK1 ★ ★ ★ ★★★★★★★★★】k★★★★★★★★★★★ SUBROUTI∼ε ATERκ1（NEQN●X●Y，H，XEND） REAし EPSN（5）●Y（NEQN） c c C c C C c c c c

5蹴｝ぞイ‖）脚：1き；｝瓢賠㈱言5｛言｝ll；喫え、），，（5、，・，PSN（’，11．

＋ 1）8，3X）） ’ 10 N＝N＋1 CAしL RK4（NEQN，XgYgH，EPSNgN） XFXi・H

2。羅辮5！iil言§ii！illξ1妾；麟§一・・（EPS−・一）

GO『「OIV εND ★★★★★★★★★★★★’★★★★★★★★★ ★ 一 ★ ★SU日ROU．TIN．E ATERK2★ ★ ★ ★★★★★★★★．★★★★★★’★★★★★★★ ミUR・E°8｛↓2ξN（；謡維10ε‖…lll；瑠呈↓講llREしEPS・（5） 1麗鼠1手｛稲），￥1郭…経碁｝諺謬段閻設号離ラli言手；1：？Sglll，，、x）．N（s・．・EPSN（・十・11曽1）・●3X）） DO 5 1＝1●NEQN ESAVε（1 ）＝◎・O YSAVE（1）＝Y（1）

5CONTINUE

セ｛SA．VE＝H 10 N＝N＋1 CONτ正NUE

20

CALし RK4（NεQKI曹X■Y●H●EPSN●N） EPSMAX＝1iOε一20 111eNEQN 40 DO RELEPSN（1）＝A8S（EPSN（1）！Y（工）） EPSMAX＝AMAX1（EPSMAXgRELEPSN（1））

㍑i鰻懸難：瓢1蹴66・・1。7。

H・BH ．P，∼、．99Y，lil∼織（、、！HSAVE）・★（1．・／4，・）．1．・）・ESAVE（1） 60 CONTINUε GO，．了O．．20 70 工F（εPSMAX・しE・EMAX）GO TO 80 H．N＝8H， 8。12（18S職．しE．〈。．5・Et・1・X））G。 T・．9・ HN＝H Go TO IQO 90 HN＝A凹IN1（BH●（2・0★H）） 100 xSAVE．＝X

］§8iill恕1灘i縮撚酬一1−一一・1・…i EQN）

rF（X●EQ・XεNO） RETURN 140 1：＝コ，NEQ卜▲ DO YsAvE〈工）r−Y（1） ll全V三｝｛1；・ξ（澗IH）、、4．。）−1．・）・EP脚） 140 CONT工N． UE HSAV”E ・H H，＝．HN GOTO 10’ 印D

(5)

C C C c C C C C c C c C c C ★★★★★★★★★★★★★★〉 ★ ．★ ★ RK4 ★ ★ ★ ★★★★★★★★★★★★★★

SUBROUTエNE RK4（NEQNgX，Y●H●EPSN●N）

iK5（6）●K6（6）●KY（6） REAし K1（6），K2（6）9K3（6）tk4（6） R、EAL εPSN（NεQ・N）gY（NεQN） ★★★★ K1 ★★★★ 、、 DO ↑O L＝1●N亘uN KY（L）＝Y（し）・ 10 CONT工NUε CAしL F（X●K「了gK1） ★★★★ K2 ★★★★ DO 20 L＝1，NEQN KY（L）＝Y（し）＋H★（4●0／’15■0）★K3（Lう 20 CO−NTINUE CALL F（X＋4●0／15●0★tt，KYgK2） ★★★★．κ3 ★★★★ bO 30 L＝1，lEQN KY（L）＝Y（L）＋ト｛★（1●0／10●0 一’kK1（L）・ト3宥O／10・Ot．K2（L））

30 CONT工NUE

CALL．F（X＋2，0／5・0★H●KY，K3） ★★★★ K4 ★★★★ DO 40 L＝1・tV EQN ’ KY（L）＝Y（L’）＋H窒（7．0／64・0★Kコ（L）＋15・0！64・0＊K2（L）＋5多0！32・U★K3（し）） 40 CONTINUE 『 ’− CALL F（X＋1・0／2・0★H，KY’K4） ★★★．★ K5 ★〕「★★ DO 50 し＝1・NEQN KY（L）＝Y（L）＋ト｛★（K1（し二）−1590／4◎0★j（2（L）＋15●0／490★K3（t、）） 50 C・ONTI卜8UE CAしL F（X＋H●KY，K’5） ★★★★ K6 ★★★★ DO 60 L＝1gNεQN ．，． KY（L）＝Y（L）＋H★（272●0／81●0★｝く1（’し）＋5．0／36・0★K2（L）．−623’5．●6！972●0★｝〈3（L）＆＋1276●0／243●∩善K5〈L））

60 CONTINUE

CALL F（X＋7●0／3・0★H●KY●’K6）工F（NN・EQ●1） GG TO 80 ＝1・NEQiO DO’70 L Y（L）＝ヤ〈し’）＋H★（7●0／24．0★K1（し）一・ 875・0！696eO・★K3（L）−i・64●0！33喰0★K4（し）＆

70 CONTヱNUE

＋1●0／48●0★K5〈し）＋27●0！5104●0★K6《L）．） GO’TO 100 80 DO 90 L＝1書tVEQ・N ＆’

・90 CONTINUE

1°° Q9S；｛9）と語巴1…‖！24．。．Kl（し）、575．0！2。88．0．、3（L）．8．0！33．O．K4（L）＆＋1●0／144●0★K5’（L）＋990／51U4●0★K6（L．）） 110 CONTLeJUE．

RETURN

END Y（L）＝Y（L）◆H★（1●0／4●0★K1（L）−1025●0！1044●0★K3（L）＋56●0．／33●0★K4（し）＋1．0／36、●0★K5〈L）＋9・0！1276●0．★K6（L））

STEP

1 2 5 4 5 6 7 8 9 10 11 12 13 14 X ・1000000E・ト60 ●3000000E＋00 ，7ξ）00f｝00εナニOC ．11’00000E＋01 ．1＼｛10000｛〕E十〇1 ・1900000ε＋01 ●23COOOOE＋01 ・270qOO〔〕EiO’1 ●3100000E＋01 ．3488453E＋01 ●3b”769r）5E＋01 ．4265358E＋01 ・4653811ε＋01 ・5000CO｛〕E＋01 Y（1）、．9983342E−01 ・2955215ε＋00 ，64．4266．1．EナOU ．8912795ε＋00 ●．9975827E＋00 ●9463913E＋OU ●7、457．855E＋、00 ●4274350ε＋00 94159878E’01 一命3399684E÷0〔｝一●6708808E＋00 −●9018303E・卜00 −99984027E＋Od ●玲9590225ε＋00 Y（2）．995GO42ε＋00 ■9S53406ε＋00 ・7649029E’←00 ．4536373E＋00 ●7074936E−01 口・3233124Eナ00 −●6663347E＋00 −・9｛〕41615E＋OO ■．Eggg244・2E＋00 −．9405493E十〇．0 −・7．417139E＋00 −●4323560E＋00 −・5857134E⇔01 ■2836705ε＋00 εPSN（1）・2777778E−07 ・1328158ε一〇5 ・4819226ε口04 ・6757171E−04 ，7628334E●04 ・7295164E−04 ・5810249E−04 ・3408012ε●04 ●4677007E−05 −・2191923ε■04 −・4461087E−04 −．6065538ε一〇4 ●・6766195ε一〇4 −．4064232E−04 EPSN（2）令27777ア8E−06 94333518E−05 ●5953344ξ・．04 ●3606730E−04 ．69067G5E“05 工・2・334Zi 62E−04 −．4991068E−04 −◆6859726E−04 巳●7645408E−04 −●6422216E−04 −●5113561E’−04 鵠．3042937E−04 一舎5188757E．−05 ．1↑45307ε一〇4 4．数値実験と結果の解釈 4．1数値実験ぜ前記のプログラムが，有効であるかどうかを検討するために，線形，非線形，単一，連立のいろいろな問

題について数値実験を行った。計算にはMELCOM

700皿を使用し，倍精度演算を用いた。（1）テスト問題。単一微分方程式

(6)

Ex・ 1 y＝cosx （解）ツ＝sinx Ex． 2 y＝：ycosx （解）or＝ e・iロ Ex・ 3 yry （解）Y ・eX Ex・ 4 yt＝一ツ（解）ツ＝e−x ツ3 Ex・ 5 y’＝− 2 （解）ッ＝1／1／x＋1 Pt（0）＝0 積分区間〔0，5〕ッ（0）＝1 積分区間〔0，5〕 Pt（0）＝1 積分区間〔0，5〕 Pt（0）＝1 積分区間〔0，5〕 Pt（0）＝1 Ex． 6 ヅ＝2ッ／（1十x）（解）ッ＝＠＋1）2 Ex・ 7 ツ’＝1一ツ2 （解）ツ＝tan hx Ex． 8 0r1＝一ッ十sin 2x 積分区間〔0，5〕 N（0）＝1 積分区間〔0，5〕 or（0）＝0 積分区間〔0，5〕 or（0）＝−0．4 （解）ツ＝（sin 2x−2cos 2x）／5 積分区間〔0，5〕 Ex． 9 yt・＝（y−−1）（xy −or−x） Pt（0）＝1．5 （解）ツー2芸裟1積分区間〔・・5〕 Ex．10 yt＝x十sin x十〇r or（0）＝0．5 （解）ツー2・・−x−・−c°s守1nエ．積分区間〔0，5〕 E・・・・yt−一

?ﾝ

ツ（2）一・（解）or＝9／（x3＋1）積分区間〔2，7〕 Ex．12 y＝or 一一 2x／ッ y（0）＝1 （解）rv＝1／互i＋1 積分区間〔0，5〕ノノ

E…3y−±一撒一y20・（・）一一・

（解）or＝一⊥ 積分区間〔1，6〕 x Ex．14 ヅ・＝（3x2十2x）／8Pt or（3）＝3 （解）ツー（号）石＋・積分区間〔3・8〕 Ex．15 yt＝レ／『’N− −2L ツ（1）＝16／9 x （解）ツー（丁＋174−tl・x−）2積分区間 Ex・16 rvt＝2エツ（解）ッ＝eX2 x十ツ Ex．17 yt＝（解）ツ＝x（log x十1） Ex．18 ヅ＝x2十⊥ x

（解）ツー芸号

．y（0）＝1 積分区間〔0，5〕 Pt（1）＝1 積分区間〔1，6〕 or（1）＝0 積分区間〔1，6〕。2元連立微分方程式

Ex

@19

o謬

（解）｛；1：二：．x

Ex

@2

試ﾇ、

（解）｛露蒜

Ex’ 21 _{o1惑、＋3。。，x} Ex’ 22 ツ1（0）＝1 ｛ Y2（0）＝−1 積分区間〔0，5〕ツ1（0）＝0 ｛ツ2（0）＝1 積分区間〔0，5〕ツ1（0）＝1 ｛ツ2（0）＝1 （解）

ﾉ：二撫霊〕

1 ツ1 （解）｛ Ex’ 23

_o

ヅ・＝了計1一2xor・

y・一一5芸、＋2xツ1

ツ1＝〆x＋1c・s（x・）ツ，＝〆針1sin＠・） y1＝−1．38ツー0．8ツ2 ｛ツ1（0）＝1 ツ2（0）＝0 積分区間〔0，5〕 Pt1（0）ニ＝−2．9995 ッ2’＝−2．16y，−1．920r2 0r2（0）＝＝4．001

（解）｛㌶：1麗㍊：lll㌔区間〔・・5〕

E苦2

氏F：≡ （；：：：：：1

ツ1

（解）｛繊積分区間〔…〕

5 ｽ霊ピ：：1；：1

（解）ツ1一c｛（・一㌃）2−21・g（・一吉）−1｝ツ・一c｛2（・一㌃）（言）積分曙5〕一（225）／（1一㌃）｝。3元連立微分方程式

6 _o妻i総

（解）

o欝i−8i−．

7 _o》≡裁

ツ1（0）＝3 ｛ツ2（0）＝4 ツ3（0）＝−8 積分区間〔0，5〕ツ、（0）＝2 ｛ツ2（0）＝0 ツ3（0）＝1

(7)

（解）

ﾛil三：：：：∵積分区㌦〕

（II）実験方法（1）固定刻みの場合結果は次のような数量Rnで要約した。

ふ虹讐欝難劃り誤差）

このRnをε，、の推定精度と呼び，その絶対値を次の三つの範囲に分けた。すなわち， 0．9≦二lR？b ≦91．1 0．5≦91Rnl≦2．0 0．1≦lRπ1三≦10．0 いま，各問題を区間〔0，5〕において固定刻みh＝ 0．1で積分すると，ステップ数は50になる。おのおののステップでの推定精度の絶対値が上の三つの範囲のいずれかに入れば，入ったそれぞれの範囲で1ずつカウントして，積分が終了したら，それらの三つの範囲に対してカウントした数を総ステップ数で割る。それに］00を掛けたものをパーセントとして表示した。表一4．1，表一4．2および表一4．3はそれぞれ，単一方程式，2元連立方程式，および3元連立方程式について上記の結果を示したものである。（2）可変刻みの場合 3節において述べた最大相対許容誤差Emax＝ 10−4 として，同じ問題を可変刻みで計算した結果について，（Dと同様な諸量を求め，表一4．4，表一4．5および表一4．6に示した。図一1∼図一7は，loglo lR。1を縦軸に，独立数エを横軸にとって数値実験結果を図示したものである。このloglo lRnlが零に近いほど推定精度が良いということになる。また，縦軸の値が1と一1のところに破線が引かれているが，これは4．1における 0．1≦lR∋≦10．0に相当している。図一1・図一2，図一3および図一4はそれぞれEx．7， Ex．8， Ex．21， Ex．27の固定刻みによる計算結果を示したものである。また図一5，図一6，図一7および図一8は，それぞれ，Ex．7， Ex．8， Ex．21， Ex．27の可変刻表一4．1 固定刻み単一微分方程式 Ex． 1 Ex． 2 Ex． 3 Ex． 4 Ex． 5 Ex． 6 Ex． 7 Ex． 8 Ex． 9 Ex．10 Ex．11 Ex．12 Ex．13 Ex．14 Ex．15 Ex．16 Ex．17 Ex．18 o．9≦二lRπi o．5≦二1Rn｝ o．1≦こ11∼nl ≦二1．1 ≦2．0 ≦二10．0 80．0％ 50．0 38．0 32．0 2．0 2．0 0．0 46．0 0．0 46．0 23．3 0．0 0．0 56．0 0．0 48．0 0．0 0．0 98．O％ 92．0 100．0 100．0 24．0 2．0 26．0 92．0 16．0 100．0 100．0 2．0 4．0 100．0 10．0 76．0 10．0 2．0 100．0％ 98．0 100．0 100．0 72．0 8．0 82．0 100．0 98．0 100．0 100．0 12．0 16．0 100．0 28．0 96．0 34．0 10．0 表一4．2 固定刻み2元連立微分方程式 Ex．19 Ex．20 Ex’21 Ex．22 Ex．23 Ex．24 Ex。25 0．9≦［Rnl≦91．1 Rn，i Rn，2 32．0％ 46．0 58．0 22．0 86．0 0．0 100．0 32．0％ 46．0 44．0 24．0 86．0 2．0 100．0 0．5≦lRni三≡92．0 Rn，1 Rn，2 100．O％ 100．0 100．0 72．0 98．0 18．0 100．0 100．O％ 98．0 98．0 68．0 100．0 8．0 100．0 0．1≦司R∋≦10．0 Rn，1 Rn，2 100．0 100．0 100．0 98．0 100．0 34．0 100．0 100．0％ 100．0 100．0 96．0 100．0 26．0 100．0 表一4．3 固定刻み3元連立微分方程式 Ex．26 Ex．27 0．9≦；＿lRnl≦91．1 Rn，1 Rn，2 Rn，3 5．0％ 5．0％ 5．O％ 2．0 10．0 14．0 0．5≦：IRn｛≦92．0 ＆，1 Rn，2 Rn，3 45．OhOo 45．0600 45．O％ 60．0 64．0 70．0 0．1≦9．IRnl≦：10．0 Rn，1 Rn，2 Rn，3 90．0％ 90．O％ 90．O％ 100．0 94．0 94．0

(8)

みによる計算結果を示す。 4．2 結果の解釈まず，固定刻みについていえぽ，Ex．6， Ex．12， Ex・13， Ex・15， Ex・17， Ex．18の6題を除いて，集積打ち切り誤差の推定値と真の集積誤差との比の絶対値が，ほとんど0．1から10．0以内に入っている。これは言い換えれば，集積打ち切り誤差の推定値が真の集積誤差と比較して，オーダーにして±1の範囲内（すなわち1／10から10倍の範囲）におさまっていることである。連立微分方程式についてもEx． 24を除い表一4．4 可変刻み単一微分方程式（Emax＝10−4） 1°・蚤！｛詞゜・箋！『・1°・芸ぽ1 ステップ数て同様なことがいえる。次に可変刻みの場合には，大きな刻み幅にもかかわらず良い推定誤差が得られ，数例を除いて，固定刻みの場合と同様に，推定誤差は大体オーダーにして真の誤差の±1の範囲内におさまっている。良い結果が得られていない問題についていえば，Ex．6Ex．18とは解関数がそれぞれッ（x）＝（x＋1）2およびrv（x）＝x3 ／2−x／2となるから，（2．3）の〉，（4）（x）が零となってしまい，誤差の推定が不可能になったと推測される。また，Ex．12， Ex．24について言えば， Ex．12 の一般解は Ex． 1 Ex． 2 Ex． 3 Ex． 4 Ex． 5 Ex． 6 Ex． 7 Ex． 8 EX． 9 Ex．10 Ex．11 Ex．12 Ex．13 Ex．14 Ex．15 Ex．16 Ex．17 Ex．18 35．7％0 23．5 35．7 21．4 11．1 0．0 9．1 28．0 22．2 26．7 25．0 0．0 0．0 16．7 0．0 14．1 0．0 0．0 71．4％ 64．7 85．7 42．9 66．7 28．6 54．6 76．0 88．9 86．7 58．3 7．7 8．3 50．0 33．3 70．3 42．9 16．7 100．0％ 88．2 100．0 85．7 100．0 71．4 100．0 88．0 88．9 93．3 100．0 30．8 50．0 83．3 88．9 93．8 85．7 50．0 14 17 14 14 9 7 11 25 9 15 12 13 12 6 9 64 7 6 含）−n．OO δ

9

−1．00 一2．00 Ex．7y’・・1−y2， y（0）ニ0 ヘー…−tSも 2．00 3．00 TIM−X 一3．00 図一1 Ex．7の数値触こおける集積打ち切り誤差の推定精度（固定刻み，刻み幅h＝0．1）・、表一4．5 可変刻み2元連立微分方程式（Emax＝10”4） Ex．19 Ex．20 Ex．21 Ex．22 Ex．23 Ex．24 Ex．25 0．9≦；＿lRnl≦！1．1 Rn，1 Rn，2 21．4％ 35．7 35．0 3．9 0．0 12．5 100．0 21．4％ 28．6 25．0 6．5 25．0 0．0 83．3 0．5≦91Rnl≦≡92．0 Rn，1 Rn，2 42．9％ 85．7 85．0 33．8 12．5 37．5 83．3 42．9％ 85．7 75．0 32．5 62．5 6．3 100．0 0．1≦91 Rπ1 sglo．0 Rn，i Rn，2 85．7％ 100．0 95．0 93．5 50．0 43．8 100．0 85．7％ 100．0 95．0 93．5 75．0 37．5 100．0 スプテ数ツ 14 14 20 77 8 16 6 表一4．6 可変刻み3元連立微分方程式（Emax＝10”4） Ex．26 Ex．27 0．9≦：lRnl≦91．1 Rn，1 Rn，2 Rn，3 6．5％0 6．5％ 6．5％ 8．3 8．3 8．3 0．5≦91Rnl≦92．0 Rn，1 ．R。，2 Rn，3 48．4％ 48．4％ 48．4％ 50．0 75．0 75．0 0．1≦lRπ1≦10．0 Rn，1 R。，2 Rn，3 90．3％ 90．3％ 90．3％ 88．3 100．0 100．0 スプテ数ツ 62 12

(9)

3．00 2．00 1．00 乞芭 6−0・00 日一1．00 一2．O．O 一3．00 Ex．8 y’ ＝−y十sin 2x， y（0）＝−0．4 2．00堂 3．00 TIME−X Ex．8の数値解における集積打ち切り誤差の推定精度（固定刻み，刻み幅hニ0．1） 3．00 2．00

2

色 6−0・00

9

一1．00 一2．00 Ex．7 y，＝1−y2， y（0）＝0

・ω∨

図一5Ex．7の数値解における集積＿3．00 打ち切り誤差の推定精度（可変刻み） 3．00 2．00 E・・21｛鱈。＋』，｛拙二1 1．00−一一一一一．一一一一一一一一一一．＿＿一一一一一一＿．一一＿一一．．＿一一＿一一一 2 za

訂∵∵≧4’°°5’°°

一2．00 図一3 Ex．21の数値解における集積打＿3．00 ち切り誤差の推定精度（固定刻み，刻み幅h＝0．1） 3．00 2．do 1．00

2

色 R−O．OO

8

一1．00 一2．00 一3．00 3．00 TIME−X 図一6 Ex．8の数値解における集積打ち切り誤差の推定精度（可変刻み）

；｛厭≡戴……1…≡1

∂ −O．00

8

・−P．00 一2・° 1 −3．OO」図一4 Ex．27の数値解における集積打ち切り誤差の推定精度（固定刻み，刻み幅h ＝O．1）

ll

ノ・x，21｛；㍍．、c。sx，・／・，2（・）一・fy・（0）＝1 1．00−一一一一一一一一一一一一一・一一一一・一一一・一一一一一一一一一一一一p・一一一一一一・一一≡一一…一

2

色 2 0−O．OO茜ti

9

TIME−X 2．00 3．00 5．00

》l

r／一1．00−一一…………一一一一一一一……一一一一…一一一一一……一「 ll：1 図一7 Ex．21の数値解における集積打ち切り誤差の推定精度（可変刻み）

(10)

y7＝−y、＋Y2 yZ＝y一2y、＋y、 Y3＝Y2｝Y3 ， y、（0）＝2 y，（0）＝O Y3（0）＝1

＿』

／／

2

色 TIME−・X

§llllヤミぎ

＼／ −2 ］図一8 Ex． 27の数値解における集積打一3．00」ち切り誤差の推定精度（可変刻み） or（x）・＝ン 1＋2・C＋ce2x（問題の初期値でc ・O） Ex．24の一般解は Yi（x）＝（1 −Cie2x）c2eX’（問題の初期値でCl＝0， c2 ・＝ 1） 1十c1θ2甜ツ2＠）＝ C2e−；v （1−Cie2x）2 である。上の一般解からわかるように，これらの問題は，初期条件によって本来零になるべき定数が丸め誤差などのために零にならず小さな値になると，それによって導入される無縁成分が数値的不安定現象を引き起こすたちの悪い問題である。したがってこの種の問題に対して失敗することはあながち悲観するに当たらない。この種の問題は，たとえば指数関数の部分が解曲線に現われるように初期値を決めた方が，もっと良い結果が得られるかもしれない。また次に，図一7を見ると，積分区間の最後のところで推定精度の値が急激に小さくなっている。これは積分区間の最後で刻み幅を強制的に推定値より小さくとったために起こったもので，積分区間の上端における集積打ち切り誤差を正当に評価するには，別の方法，たとえば，補間などによらなければならないだろう。一般に刻み幅を小さく変更することは過小評価になる要素を持っているので，この点刻み幅変更に用いるアルゴリズムはなお検討が必要であろう。

5．結

び以上の考察から知られるように，提案する数学ソフトウェアは，Merluzzi−Brosilow法の主誤差関数に関するある仮定が満たされない特殊な問題には適用できないが，きわめて有用であると考えられる。このプログラムの使用により，数値解のもつ真の誤差の範囲をほぼ推定することが可能である。このことは局所打ち切り誤差の推定による従来の方法では不可能であって，本プログラム存在の意義もまたここにあると考えられる。一一謝辞本研究を進めるにあたり，有益な御討論をいただいた山梨大学田口東講師に深謝する。文献 1） R．H． Merson：An operatioral method for the study of integration processes， Procedings of aSvmposium on Data processing， Weapons Research Establishment， Salisbury， South Australia，1957 2）田中正次：Rung−Kutta法の打ち切り誤差にっいて，情報処理Vo1．17， No．12．1976． 3） D．G． Bettis：Efficient embedded Runge−Kutta methods， Lecture note in mathematics 631， Numerical treatment ordinary differential equations， A． Dold and B． Eckmann（ed．）1976 4） J．H． Verner：Explicit Runge−Kutta methods with estimates of the local truncation error， SIAM Journal on Numerical Analysis Vo1．15， No．4，1978 5） P．Merluzzi， C． Brosilow：Runge−Kutta Integ． ration Algorithms with Built−ln Estimetes of the Accumulated Truncation Error， Comp． uting 20， 1−16（1978）

集積打ち切り誤差の推定ができる初期値問題の数学ソフトウェア 利用統計を見る