3． 常微分方程式の数値解法の諸問題
3.1　台形法の不安定
. . Milne 法は Runge K. 法に比べて精度が悪く、(1.1)の数値解を x=0 から h=0.05 で求めていくと x=5 以上では有効桁はなくなり、正負に激しく振動しながら振幅が増大していく。これは誤差の伝播では到底説明できないので数値的不安定現象と呼ばれ理論的にも証明できたと見なされている様であるが、著者は数値的不安定と定義すべき現象ではなく、常微分方程式の数値解法が低次の固定次数近似で、且つ、固定区分幅であることに起因する解法の適用限界を越えているのであり、この様な解法で正しい解が得られると考えることが本質的に間違っていると考える。
. . Cp+1[a, b]なる任意の関数を Newton の前進補間公式や Tayler 展開で近似する場合、関数の全定義域に渡って同一区分幅、同一次数で、且つ、同一精度で近似することは一般的には出来ない。数値積分においては被積分関数により精度が異なるだけでなく、同一関数でも場所により著しく精度が異なり、１ステップの積分で有効桁の無くなる場合すらあることは前節に示した。常微分方程式の数値解法も数値積分であるから同様の誤差を生じ、単なる数値積分とは異なり、その誤差は次のステップに出発値の誤差として影響を与える。低次の固定次数、固定区分幅の解法ではこの誤差の累積を抑えられないことも適用限界を越える原因となり、数値的不安定現象と呼ばれる現象を生ずるのである。
. . 常微分方程式 dy/dx=f(x, y)を任意の点 x0からx1=x0+h 迄 Euler 法で予測し、台形法で修正して解くと、正確な f1 に収束できたとき相対誤差は(2.17)で与えられる。

Eq3_1___________(3.1)
の任意の x0 における１ステップの解の相対誤差は、(2.17), (2.18)より、
Eq3_2___________(3.2)
従って、相対誤差は h のみで決まり、全区間に渡り一定であるから１ステップの解は x0 を何処にとっても h=0.1 では前節表 2.3 により有効桁４桁である。
. . １ステップの相対誤差を ε と置き、x0 から xn=x0+nh 迄、逐次に解 yk の近似値 zk を計算したときの累積誤差は次のようになる。
x0 において y0=a の解は、__________ Eq3_2a
(x1, z1)を初期値とする解は、________ Eq3_2b
(x2, z2)を初期値とする解は、________ Eq3_2c
以下同様にして、(xn−1, zn−1)を初期値とする解は、
Eq3_2d
新値はy(xn)=ae−nh であるから累積誤差の相対値は、
Eq3_3___________(3.3)
故に、x0=0 から xn=10 迄 h=0.1 で解くと n=100 であるから有効桁は２桁失われる。しかし、最初に y1 の近似値に生じた誤差 εy1 は表 2.3 から約 −7.9×10−4 であるが、最終誤差 nεyn は表 2.3 より約 −4.0×10−8 であるから最終誤差は最初の誤差に比べて遥かに小さいものであり、誤差は伝播しないことを示している。初期の誤差は急速に減衰するが累積誤差の減衰が積分値の減少に追い着かないのである。実際に、x0=0, y0=10, h=0.1 とし Euler 法で予測し台形法で修正する方法により微分方程式(3.1)を解くと、x=0.1 では y=9.04762、x=10 では y=4.50226×10−4 である。表 2.3 の真値と比べれば、前者は有効桁４桁、後者は２桁であるから上記の理論誤差とぴったり一致している。
. . この解の精度を改善するには h を小さくするよりないが、n が h に逆比例して増加するので ε がそれ以上に改善されなければ効果がない。この場合は、h=0.05 にすると(3.2)よりε=−1.0681×10−5 となるので１桁程度の精度改善が望め、さらに h を小さくすれば全域で単精度の解を得ることもできる。しかし、理論計算と異なり、h を小さくしすぎると x+h の計算で h が情報落ちすることと n の増加により、x の大きいところでの精度は理論値程は改善されない。
Fig3_1
図 3.1
逆に、h を大きくすると１ステップの誤差は大きくなるが総ステップ数 n が h に反比例して少なくなり累積誤差の影響は少なくなる。しかし、新たな問題が生ずる。積分では被積分関数 f が既知であるから問題はないが、微分方程式の場合は Euler 法による予測値 yk+fkh が負側に飛び出し修正のための正しい fk+1 が求まらない事が起こる。実際の数値解法では図 3.1 に示すように、この負の解を用いた xk+1 における微係数を fk+1 として１回目の修正を行うことになる。繰返し修正の結果 zk+1>0 に収束すれば解 yk+1 の近似値であり、(3.2)、(3.3)で１ステップの誤差と累積誤差を評価してよい。その限界は予測値 yk+fkh が −yk となったときである。このとき微分方程式(3.1)より fk+1=−(−yk)=−fk であるから平均勾配は零となり１回目の修正値は yk に等しく p 点になる。p 点における微係数は(3.1)より fk+1=−yk=fk であるから平均勾配も fk となり、２回目の修正値は −yk で q 点となる。従って、２安定となり解は永久に収束しない。予測値が −yk より大きいときは fk+1=−yk+1<−(−yk)=−fk であるから１回目の修正値は p 点より下になり、２回目の修正では fk+1=−yk+1>−yk=fk となるから修正値は q 点より上になる。従って、p, q 間の１点に収束すると考えてよい。しかし、この様な領域では修正回数の最大値設定による打ち切りを必要とする。
. . 故に、(3.2)、(3.3)の使用し得る限界は x に関係なく次のようになる。
Eq3_3a
又、予測値が零又は負になるのは yk+fkh≤0、即ち、h≥1 であるから、この微分方程式の場合は h<1 では不安定現象は生じない。例えば、h=0.5 にすると ε=−1.3459×10−2 であるから x=10 では nε=−2.6918×10−1 となり有効桁１桁の解が得られる。
. . 微分方程式が、___________ Eq3_4___________(3.4)
の場合には、任意の x0 から１ステップの解の相対誤差は(2.20)より、
Eq3_5___________(3.5)
この場合は。(2.20)と同じで１ステップの相対誤差が出発値の x 座標の値により異なるので xn−1 から出発する n 番目のステップの相対誤差を εn で表せば(3.3)を導いたのと同様にして n ステップの累積誤差を求めることができる。
x0 において y0=a の解は、______ Eq3_5a
(x1, z1)を初期値とする解は、____ Eq3_5b
Eq3_5c____
(x2, z2)を初期値とする解は、____ Eq3_5d
Eq3_5e
同様にして、(xn−1, zn−1)を初期値とする解は、____ Eq3_5f
Eq3_5g
Eq3_5h ____により{ }の中を整理すると、
Eq3_5i
Eq3_5j
真値はEq3_5kであるから累積誤差の相対値は、
Eq3_6___(3.6)
εk は(3.5)の x1 を xk、x0 を xk−1 に換えたものであり、２次の項が省略できる範囲では累積誤差の相対値は総和となり数式計算でそれを求めることも出来るが、省略できない範囲では(3.6)の省略のない式で２倍精度計算により求めなければならない。
. . この場合には、h<1 でも Euler 法による予測値が零又は負になることが起こる。
Eq3_6a
h=0.1 では xk≥10、h=0.5 では xk≥2 において予測値は零又は負になるが、表 3.1 に示すように、その初期においてはまだ解は有効桁を持っている。予測値が −yk となったとき、
Eq3_6b
従って、１回目の修正値は図 3.1 の p 点より上、即ち、yk+1>yk となり、２回目は、
Eq3_6c
であるから修正値は q 点より下になり、繰返し修正により発散してしまう。そのような領域の概略値は yk+fkh≤−yk より xk≥2/h である。h=0.1 のときこの領域は xk≥20 であり、２進法の計算機では通常指数部の制限からここまで計算できないので不安定現象は起きないことになる。しかし、それは１ステップの相対誤差が十分小さい領域の場合であり、そうでない領域では yk に代えて近似値 zk を用いなければならない。従って、１ステップで有効桁が失われるような xk=10 の近辺では不安定現象は起こる。h=0.5 のときはこの領域は xk≥4 であり、表 3.1 の実際の数値解が示すように不安定現象が起きる。表 3.1 の累積相対誤差の絶対値が xk=4 に於いて１より大であるのがそれを示し、この後は(3.5)は使用出来ないので表 3.1 の赤字表示より後の(3.6)による累積相対誤差の理論値は意味を持たない。
表 3.1 x0=0 より xn 迄の累積相対誤差と不安定
xn累積相対誤差の理論値h=0.5 の数値解真 値
h=0.1h=0.5修正打切20回修正打切21回
h
1.0
2.0
3.0
4.0
4.5
5.0
5.5
6.0
6.5
7.0
7.5
8.0
9.0
10.0 		1.25208×10−5
1.07736×10−3
1.74273×10−3
−6.56812×10−3
−3.90923×10−2
−7.06485×10−2
−1.15783×10−1
−1.76409×10−1
−2.53271×10−1
−3.45328×10−1
−4.49270×10−1
−5.59424×10−1
−6.68277×10−1
−8.50824×10−1
−9.55597×10−1		 8.14845×10−3
3.14543×10−2
5.77183×10−2
−2.96036×10−1
−1.05177×10 0
−9.41751×10−1
−1.15108×10 0
−2.70822×10−1
−6.94063×10 0
7.62576×10 1
−1.54732×10 3
4.64174×10 4
−2.04983×10 6
−1.20181×1010
−2.85004×1014		 8.88889×10 0
6.22222×10 0
1.41414×10 0
9.25918×10−2
−6.21529×10−3
−6.55409×10−2
−6.79378×10 0
5.52986×103
−2.94219×10 7
8.77619×1011
−1.28826×1017
8.31974×1022
−2.14462×1029
Over flow
:
8.888899×10 0
6.222229×10 0
1.414149×10 0
9.372549×10−2
1.408009×10−2
1.670349×10−1
2.162209×10 1
2.419389×10 4
1.930789×10 8
9.358699×1012
2.404099×1018
2.911109×1024
1.500829×1031
Over flow
:
8.82496×10 0
6.06530×10 0
1.35335×10 0
1.11089×10−1
3.35462×10−3
4.00652×10−4
3.72665×10−5
2.69957×10−6
1.52299×10−7
6.69158×10−9
2.28973×10−10
6.10193×10−12
1.26641×10−13
2.57675×10−17
1.92874×10−21
. . 実際の数値解法に於いては各ステップの発散は修正計算の繰返し回数制限により打ち切られ停止し、そのときの値を解として次のステップに進む。従って、制限回数が奇数であれば修正計算は常に図 3.1 の p 点側で停止するから解は正側に異常増大し、いわゆる幻影解となる。制限回数が偶数のときは q 点側、即ち、出発値と反対符号の側で停止するので解は負になり、次のステップは負側の解を求めることとなるが、修正計算の発散は反対側の正側で停止する。以後同様にしてステップ毎に符号を変え、激しく振動しながら異常増大する解となる。何れの場合も異常増大の速度は修正計算の繰返し制限回数に依存し、制限回数が多い程急速に増大する。表 3.1 には２通りの修正回数の例を示す。
. . 従来、これは数値的不安定現象と呼ばれているが、上記に述べたように、Euler 法で予測し台形法で修正する解法を適用限界を越えて強行した結果であり、表 3.1 の赤字表示以降の数値は解として何の意味も持っていないし、幻影解とか数値的不安定現象と呼ぶべきものでもない。この適用限界の正確な値は図 3.1 に示すように２回目の修正点が q 点、即ち、予測値に一致するときの x の値であるが、h だけでなく微分方程式の形に依存する。
. . このような領域に於いても新たに初期値を設定し、h を小さくして計算すれば正しい解を求めることが出来るが、それは h に反比例して解法の使用可能範囲が広がるからである。又、よく知られているように、x の大きい方から零に向かって解を求めた場合には不安定現象が生じないが、微分方程(3.4)の場合にはこの範囲では Euler 法による予測値が負になることが起きないし、１ステップの相対誤差(3.5)が正であるので近似解が負になることも起きないからである。

3.2　Runge K. 法の不安定
. . 微分方程式 dy/dx=f(x, y)を Runge K. 法で解くとき、区分幅 h の２区間を１ステップで解くと公式は(3.7)で表される。②式は中間点 x1=x0+h の関数値 y1 を Euler 法で予測し、その点の f を求めることを意味し、③式は後退 Euler 法 y0=y1−f1h により修正して再度その f を求めることを意味する。④式は修正された f により x2=x0+2h の関数値を中点法で予測して f を求め、⑤式は中間点 x1 の予測値及び修正値に対する f、即ち、②③の平均値をその点の f として Simpson の公式により積分することを意味する。従って、これら３点の f の正確値が求まったとき、解の精度は Simpson の積分公式の精度に等しくなり、⑤式の y2 を z2 と置けば Runge K. 法の１ステップの相対誤差は(2.21)で与えられる。
. .

Eq3_7___________(3.7)
本来の Runge K. 法は１区間を１ステップで解くから(3.7)の h を h/2 に置き換えたものであるが、これは Milne 法の修正式の積分区間の半分を Simpson の積分公式で積分することになり著しく Runge K. 法に有利になるので本著では(3.7)を用いる。これにより、３点の f の正確値が求まったとき１ステップの解の精度は同じになる。
. . 微分方程式(3.1)の場合には、任意の x0 における１ステップ２区間の数値解の相対誤差は正確な f0, (f1+f2)/2, f3 が求まったとして(3.1)の解を(2.21)に代入すれば(2.22)と同じになる。従って、全域に渡って相対誤差は同じで、h のみにより定まり、h=0.1 では有効桁７桁(−1.22855×10−7)、h=0.05 では有効桁９桁(−3.65068×10−9)である。従来の比較の仕方では Runge K. 法は Milne 法に比べ有効桁で２桁有利となる。
. . n を偶数とし、x0 から x0+nh 迄逐次に解を求めた場合の累積相対誤差は台形法の場合と同様にして求めることが出来、１ステップで２区間進行するから(3.3)の導出過程でステップ数を n/2、区分幅を 2h に置き換えればよい。故に、
Eq3_8___________(3.8)
従って、h=0.1 で x=0 から 10 迄逐次に解を求めると、略 1.5 桁有効桁を失う。
. . 実際には、Runge K. 法の x1, x2 における f の値は正確には求まらないので１ステップの相対誤差は(3.7)から求めなければならない。
Eq3_8a
これらを(2.21)に f1=0.5(f1+f2), f2=f3 として代入すれば、
Eq3_9___________(3.9)
h=0.1 とすれば、１ステップの相対誤差は 3.15153×10−6 となるから、(2.22)で計算される理想的な場合に比べて有効桁で１桁精度が悪くなる。h=0.05 なら 0.905842×10−7 となり、h=0.1 の理想的な精度に略一致する。従って、Runge K. 法は Simpson の 1/3 則の１積分区間、即ち、２区間を２ステップで解いたとき略対等の精度となるのであり、被積分関数を２次関数で近似した積分精度より悪い。2.2 節において Tayler 展開の精度について述べたように、Runge K. 法が Tayler 展開で h4 迄一致するといっても、それは h→0 のとき意味を持つのであり、h=0.05、即ち、本来の Runge K. 法の区分幅で 0.1 程度では差分法に基づく４次式の近似精度は得られない。
. . この微分方程式の場合、(3.7)において中点 x0+h の予測値と修正値及び x0+2h の予測値のどれかが零又は負となったとき正常な結果は得られなくなる。x0+h における予測値は y0+f0h=(1−h)y0 であるから h<1 では正、修正値は y0+f1h=(1−h+h2)y0>0 であり、最後の予測値 y0+2f2h=(1−2h+2h2−2h3)y0 は h<0.6…では正である。又、累積相対誤差が正であるから近似解は常に正である。従って、１ステップの区分幅が 2h<1 であればいわゆる幻影解は生じない。
. . 微分方程式(3.4)の場合には、任意の x0における１ステップの数値解の相対誤差は正確な f0, (f1+f2)/2, f3 が求まったとして(3.4)の解を(2.21)に代入すれば、Simpson の 1/3 則の相対誤差(2.23)と同じになる。従って、精度は x0 により異なり、h=0.1 のとき x0=0 においては有効桁７桁であるが(−1.67501×10−7)、x0=10 においては真値から出発して１ステップで有効桁２桁(−3.03927×10−2)しかない。これは誤差の伝播ではない。h=0.05 ではそれぞれ有効桁９桁(−2.60742×10−9)と４桁(−5.37379×10−4)となり、従来の比較の仕方では Runge K. 法は Milne 法 に比べて有効桁で２桁有利となる。
. . n を偶数とし、x0 から x0+nh 迄逐次に解を求めた場合の累積相対誤差は台形法の場合と同様にして求めることが出来､１ステップで２区間進行するから(3.6)の導出過程でステップ数を n/2、区分幅を 2h に置き換え、εk は k を偶数とし(2.23)の添え字 2 を k に、1 を k−1 に、0 を k−2 に置き換えればよい。故に、
Eq3_10___________(3.10)
h=0.1 で x=0 から 10 迄逐次に解を求めた場合の累積相対誤差は(3.10)の省略のない式を２倍精度で計算すれば−1.79450×10−1 となるので、数値解の有効桁は１桁となってしまう。
. . この場合も実際の Runge K. 法の x1, x2 における f の値は正確には求まらないので１ステップの相対誤差は(3.7)を用いて(3.9)と同様に求めなければならない。
Eq3_10a
これらを(2.21)に代入すれば、
Eq3_11
この結果の添字 2 を k に、1 を k−1 に、0 を k−2 に置き換えたものが(3.10)の εk になる。h=0.1 及び 0.25 に対して、x0=0 から xn=10 迄、逐次に数値解を求めた場合について(3.10)の省略のない式により２倍精度計算で求めた累積相対誤差の理論値及び実際の数値解を表 3.2 に示す。この理論相対誤差は各ステップの初期値が真値である場合の(3.11)で計算される εk を用いて(3.10)で計算したものであり、実際の数値解では各ステップの誤差を持つ解を次のステップの初期値とするので εk の大きい所では累積相対誤差は更に大きくなる。
表 3.2 dy/dx=−xy ,___y0=10 (ステップ幅=2h)
xnh=0.1h=0.25 (*の真値=8.82496)真 値
理論相対誤差数 値 解理論相対誤差数 値 解
2h
1.0
2.0
3.0
4.0
5.0
6.0
7.0
8.0
9.0
10.0		 −6.77422E−9
1.15426E−6
1.75014E−4
2.43666E−3
1.61592E−2
7.35364E−2
2.76247E−1
1.01209E00
4.65713E00
4.29526E01
1.57251E03		 9.80199____
6.06531____
1.35359____
1.11361E−01
3.40883E−03
4.00069E−05
1.94372E−07
4.60717E−10
7.16428E−13
1.13255E−15
3.03491E−18		 −1.12447E−5
−5.98691E−5
8.54474E−3
1.77018E−1
2.11227E00
5.59433E01
1.23652E04
3.11624E07
7.77049E11
1.64140E17
2.58049E23		 8.82487 *
6.06494
1.36492
1.30755E−01
1.04405E−02
2.12208E−03
1.88338E−03
7.13537E−03
9.84067E−02
4.22950E00
4.97713E02 		9.801986
6.065306
1.353352
1.110899E−01
3.354626E−03
3.726653E−05
1.522998E−07
2.289734E−10
1.266416E−13
2.576757E−17
1.928749E−21
. . 赤字表示部に於いて理論相対誤差は１を越えているが、台形法とは異なりその後も数値解と真値から求めた相対誤差に一致している。しかし、h=0.1 の場合は赤字表示部以降も減少関数となっているが、x=10 では相対誤差は 103 にも達し、h=0.25 の場合には x=6 以降は増加関数となり x=10 ではその相対誤差は 1023 にも達している。このとき数値解の値は 497、即ち、初期値 10 の 50 倍である。従来、これらは幻影解と呼ばれているが、Runge K. 法の適用限界を越えているのであり、数値解としての意味は全くなく、幻影解と呼ぶべきものではない。
. . この微分方程式の場合の Runge K. 法の適用限界は中間点の予測値が零又は負になったときで、yk+f0h=(1−xkh)yk>0 により xk<1/h となる。h=0.25 では x<4 で良く合っているが、h=0.1 では x<10 となり限界が広すぎる。実際にはこれより前に xk+2h の点の予測値が零又は負となり(3.7)の f03が正になるので、解の精度は有効桁１桁にも満たず使用に耐えるものではない。その限界は、
Eq3_11a
. . 故に、h=0.1 では x<6 となるが累積相対誤差が正であることによりこの限界は x=7 迄広がる。h=0.25 の場合には限界が解の変曲点に近いことや区分が粗いことにより、この推定は狭すぎる。この限界を正確に決定することは困難であるが、この後も真の解に似て減少関数であることから、概略 x<1/h を以って限界と考えてよい。従来、h=0.1 でこの区間を１ステップで解く Runge K. 法とこの２区間を Simpson の 1/3 則で修正して解く Milne 法とを比較して Runge K. 法は安定であるとする議論を行っているのは誤りである。対等の比較をするためには Runge K. 法も１ステップで２区間を解く、即ち、従来の解法の表現では h=0.2 で解かなければならないのであり、この場合には Runge K. 法も不安定現象と言われる現象を生ずる。h=0.1 を１ステップで解く場合でも、第三章第 5.1 節の表 5.1 に示すように２階、３階等の微分方程式では１階の Milne 法と同様に暴走する。

3.3　Milne 法の不安定
. . Milne 法の１ステップの予測子及び修正子は(3.12)で与えられる。予測子①は x0, x1 及びそれより前の２点 x−2, x−1 における関数値が既知であることを前提とし、被積分関数 f をそれらの点で既知である４点を通る３次関数で近似し、区間[x−2, x2]を積分することにより y2を求めるものである。しかし、修正子②は予測点を含む３点を用いて[x0, x2]を Simpson の 1/3 則により積分するものであり、修正の繰返しにより正しい f2 に収束したとき y2 の精度はSimpson の 1/3 則の精度に一致する。

Eq3_12___________(3.12)
従って、累積相対誤差は(3.7)で表す Runge K. 法の理想的な場合に一致する筈であるが、実際は Runge K. 法の精度には及ばない。それには二つの原因がある。
. . Milne 法は修正子で２区間積分して１区間しか出発点が進行しないのでステップ数が２倍になるのが一つの原因で、これにより累積相対誤差は２倍になるが、もっと重大な原因は修正子の積分区間の中間点の値を前のステップで求めた近似値に固定することにあり、これにより著しい誤差を生ずる。微分方程式(3.4)の場合、Runge K. 法は１ステップの計算で ε>0 の誤差が入っても次のステップの初期勾配が f=−x(y+ε)となり、Euler 法の計算値を小さくする、即ち、解の値を小さくし、誤差の補正効果がある。一方、Milne 法は誤差を増加する作用を持つので誤差の大きい所では急激に精度が悪化する。
. . 正しい値 y−2, y−1, y0, y1 を与えて y2 を予測し、y0, y1, y2 を用いた y2 の修正が収束したとする。この収束値と真値 y との誤差をε>0とすると、y2=y+ε である。次のステップは h だけ進行するからこの修正値は y1 になり、この点 x1 における勾配 f1 と真値 y に対する勾配 f の関係及び(3.12)①の予測値は次のようになる。yf は x1 における勾配に真値 f を用いた予測値、即ち、最初のステップと同様に真の出発値４点を与えた場合の予測値である。その結果、y2 の予測値は f1 の誤差に依存する誤差が加算される。
Eq3_12a
この予測値を(3.12)②へ代入した修正値は f1 の誤差が加算され次のようになる。
Eq3_12b
y(−x2yf)は f0, f, −x2yf を②へ代入した修正結果であり、最初のステップと同様に真の出発値４点を与えた場合の修正値である。y2 の修正値において、予測値の持っていた f1 の誤差に依存する誤差は −x2h/3 が掛けられ減少するが、新たに f1 の誤差に依存する誤差が加算され、以後、修正の度に前回迄の誤差は係数 −x2h/3 を掛けられ、新たに f1 の誤差に依存する誤差の加算があるのでこの誤差は収束しない。真の出発値４点を与えた場合の解が一桁も二桁もの有効桁を失う所ではこの誤差が著しい悪影響を与える。y2 はこの他に真の出発値４点を与えた場合の解と真の解との差による誤差を持つ。これ等の誤差はステップ数だけ積み上げられるので、解の精度が落ち始めると急速に精度は悪化し､暴走するのである。

3.4　解の存在と一意性
. . 微分方程式(3.13)は従来の数値解法では解くことができないとされている。その理由は関数 √y, (y≥0)が y=0 で Lipschitz 条件を満たしていないので解が一意でないからであるといわれる。確かにこの解は(3.14)に示すように無数に存在するが、求積法では解くことができるのであるから、解が一意でないことは数値解法で解けない理由にならない。実際、台形法や著者の解法では正確に解くことができる。

Eq3_13___________(3.13)
Eq3_14___________(3.14)
. . 求積法においても自然な解は y=x2/4 または y=0 であり、微分方程式に照らして考えてみると、これらを組み合わせた(3.14)も全て解になっていると解釈できるのであって厳密には求積法で求まるとはいえない。これらは人工的な解であり、不自然な解というべきである。本来、微分方程式は自然な物理現象を数式的に記述することを目的とするものであり、これらアナログ現象は高次微分も連続であるべきで、これが不連続な解は自然界には定常的には存在し得ない不安定な解である。(3.14)は２次導関数が x=c で不連続であり、不確定な定数 c を含むことが安定でないことを意味する。
. . 数値解法においては、Euler 法では２次、Runge K. 法や Milne 法では３次の導関数まで連続でないと剰余項が意味を持たないから(3.14)の解を直接求めることはできない。しかし、y=x2/4 または y=0 は求まらなければならないが、これらの解法では後者は求まるが前者は求まらない。x=0 で Lipschitz 条件を満たさないが、固定点であるから差し支えないのであり、前者が求まらないのは単に初期値および初期勾配が共に零であることによる。Euler 法や Runge K. 法では初期勾配 f0=0 の場合には x1=x0+h の関数値は y1=y0 であり、y0=0 であるから y1=0、この点の勾配は(3.13)により y のみで決まるから f1=0 となる。従って、全ての点で y=0 となる。Milne 法では出発値をこれらの方法で求めるからやはり同じであるが、別の方法により初期値以外も零でない正しい出発値を与えれば y=x2/4 を正確に求めることができる。
. . Eq3_14a_____を与えれば、_____ Eq3_14b
予測値は、_______ Eq3_14c_______よって、_______ Eq3_14d
修正して、_______ Eq3_14e_______よって、収束。
. . 求積法においては　y=0　ではないという仮定を用いて解を一意に制限しているのであるから数値解法でも恒等的には零でないという仮定を用いればよい。上記の　Milne　法の出発値はこれを意味している。初期値のみから出発できる解法では　y1=ε, (>0)とすれば良いのであるが、反復修正のできない Euler 法や Runge K. 法では正しい結果を得ることはできない。従って、これらを用いる Milne 法等の P-C 法では正しく解くことができない。これを解く最も簡単な方法は台形法により反復修正することであり、収束は遅いが区分幅 h に関係なく正しい解に収束する。
. . y0=0, y1=ε であるから(3.13)により、f0=0, f1=ε1/2 である。故に、
Eq3_14f
. . εは任意であるが、ε=h2 とするのが収束は速い。f1=ε1/2 の後退 Euler 法により y0=0 を求めると、ε−hf1=0 となることによるものであるが、更に、f0=0 と f1 との平均値を用いると2ε−hf1=0 となり、正しい解 ε=h2/4 が求まってしまう。正しい f0, f1 を用いた台形法の積分値を近似解とし、１ステッブの解の相対誤差を(2.17)により計算すると零であり、台形法で真の解に収束するのは当然といえる。従って、台形法の反復修正により出発値を決める P-C 法では y=0 及び y=x2/4 を求めることができ、求積法と同様に微分方程式(3.13)に照らして考えれば(3.14)も解であることは容易に分かる。Euler 法や Runge K. 法及びこれらを用いる Milne 法等の P-C 法ではこれを解くことができないのは Lipschitz 条件と無関係であり、解法として問題があるというべきである。
. . 台形法で y≡0ではない解を求めるためには Euler 法による予測値に対し(3.15)の条件をプログラムすればよい。数値解の実例を表 3.3 に示す。表の真値が x に対する理論値と僅かにずれているのは x の計算機内部の２進値が単精度の丸めにより真の x に対し僅かにずれているが、外部十進変換の丸めにより x の外部表現にそれが現れないからで真値の誤りではない。
Eq3_15___________(3.15)
表 3.3____dy/dx=√y ,____y(0)=0____(h=0.2)
x数 値 解真　　値x数 値 解真　　値
0.2
0.4
0.6
0.8
1.0		 1.000001E−2
4.000003E−2
9.000005E−2
1.600000E−1
2.500000E−1		1.0000000E−2
4.0000001E−2
9.0000007E−2
1.6000000E−1
2.5000000E−1		 1.2
1.4
1.6
1.8
2.0		3.600001E−1
4.900001E−1
6.400001E−1
8.100001E−1
1.000000E00		3.6000002E−1
4.8999998E−1
6.4000001E−1
8.1000006E−1
1.0000000E00
. . 解が恒等的には定数でない微分方程式も、初期勾配が零であり右辺に x を含まない場合には Euler 法や Runge K. 法及びこれを用いる P-C 法では解くことができない。又、右辺に x が含まれていても Euler 法では正しい解は得られない。
. . dy/dx=x, y(0)=0 の解は y=x2/2 であるが、通常は、f0=0 であるから x1=h では y1=f0h=0 となってしまう。f は y を含まないからこれ以降の f は正しい値となるが、解は真の解を平行移動した y=x(x−h)/2 となる。即ち、dy/dx=x−h/2, y(0)=0 の解に化けてしまう。この解は表 3.4 に示すように、h を如何に小さくしても元の微分方程式の数値解とは言えない程の隔たりを持つ。これは桁落ち、情報落ち、丸め誤差やそれらの伝播ではなく、初期勾配が零である場合には Euler 法の[x0, x0+h]の積分の相対誤差が −1 であることに原因がある。次式及び表 3.4 に示すように x の大きい所では h の影響が小さく、数値解の精度はよくなる。この場合も、y1=ε>0 を仮定すれば f1=x1=h、後退 Euler 法ε−f1h=0 により ε=h2 となり、台形法による修正で正しい解に収束する。実際、f1 と f0=0 との平均値を使えば 2ε−f1h=0 により簡単に正しい解が得られる。
Eq3_15a
表 3.4____dy/dx=x ,____y(0)=0____(Euler)
x数 値 解y=x(x−h)/2δy=x2/2
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
:
10.0		 0
1.000000E−2
3.000000E−2
6.000000E−2
1.000000E−1
1.500000E−1
2.100000E−1
2.800000E−1
3.600000E−1
4.500000E−1
:
4.950000E01		 0
1.000000E−2
3.000002E−2
6.000001E−2
9.999999E−2
1.500001E−1
2.099999E−1
2.800000E−1
3.600003E−1
4.499999E−1
:
4.950000E01		 −1.0
−0.5
−0.333333
−0.25
−0.2
−0.166667
−0.142857
−0.125
−0.111111
−0.1
:
−0.01		 5.0000001E−3
2.0000000E−2
4.5000003E−2
8.0000002E−2
1.2500000E−1
1.8000001E−1
2.4499999E−1
3.2000001E−1
4.0500003E−1
5.0000000E−1
:
5.0000000E01
. . dy/dx=1/y, y(0)=0 も y=0 で Lipschitz 条件を満たさないが、固定点であるからこれが解けない理由にならない。これが解けない理由は f0=∞ にあり、f0 を用いる解法では解くことはできない。求積法でも同じであり、明言はしないが y≠0 を仮定して一般解 y2=2x+c を求め、これが x=0, y=0 でも成立するなら c=0 となるから y2=2x が解であるとしているのであり、数値解法でも x=0, y=0 を避けなければならない。
. . x=0 を避けて x0=h, x1=2h とすると Eq3_15b Eq3_15c であるから、正しい f0, f1 を用いた台形法の積分値を近似解として、１ステッブの解の相対誤差を(2.17)により計算すると次式のε1 となるので、h に関係なく有効桁２桁の精度しか得られない。
Eq3_16_______(3.16)
h=0.1 で有効桁７桁の初期値 y0 を与えて解いた数値解を表 3.5 の A に示す。最初の１ステップで解の有効桁は２桁に落ちてしまうが、その後はその精度を保ち、x の増加と共に精度は改善され、x=100 では有効桁５桁となる。最初の１ステップで有効桁２桁に落ちてしまうのは台形法の欠点ではあるが、同表の B に示すように初期値の精度は有効桁２桁あればよいということでもあり、x0=0 から解くには有利な点でもある。x の小さいところでの精度を上げるには h を小さくすれば良い。h=0.001, x0=h で解けば、最初の１ステップの精度は有効桁２桁であることに変りはないが、x の増加に伴う精度の改善は表 3.5 に比べて２桁増加し、x=1 で有効桁５桁の解が得られる。
表 3.5 dy/dx=1/y, h=0.1, x0=h表 3.6 h=0.1, x0=h×10−6
xAB真　値x化けた解真　値
0.10.4472130.450000(初期値)10−74.47213×10−4(初期値1)
0.2
0.3
0.4
0.5
:
1.0
2.0
3.0
4.0
5.0
:
:
10.0
:
:
100.0		 0.637454
0.779994
0.899673
1.00500
:
1.41819
2.00297
2.45196
2.83058
3.16421
:
:
4.47352
:
:
14.1426		 0.639319
0.781506
0.900981
1.00617
:
1.41902
2.00355
2.45243
2.83100
3.16459
:
:
4.47378
:
:
14.1427		 0.6324555
0.7745966
0.8944271
1.0000000
:
1.4142135
2.0000000
2.4494897
2.8284271
3.1622776
:
:
4.4721361
:
:
14.142135		 0.1
0.2
0.3
0.4
0.5
:
1.0
2.0
3.0
4.0
:
:
10.0
:
:
100.0		 111.804
111.805
111.806
111.807
111.808
:
111.812
111.821
111.830
111.839
:
:
111.893
:
:
112.694		 (初期値2)
111.8052
111.8060
111.8069
111.8078
:
111.8123
111.8212
111.8302
111.8391
:
:
111.8928
:
:
112.6942
x0>0 を任意の位置にとると、１ステップの相対誤差は、
Eq3_16a
となり、x の大きい所では x0 と x1 の比は略１であるからε1 は略 h/2x1 となるので１ステップの相対誤差が小さくなること、関数値が大きくなることにより初期の誤差は影響がなくなることによる。しかし、必要な精度を得るために h を小さくすれば h に反比例して著しいステップ数の増加を招く。
. . h より小さい点を x0 に選ぶと１ステップの相対誤差が著しく大きくなるのでこの微分方程式を正しく解くことは出来ない。x0=ah, (a<1)とすれば、
Eq3_16b
Eq3_16c であるからx≥x0, y≥y0 では Lipschitz 条件を満たしているが、a=10−6 にとると１ステップの相対誤差は約 250 となり正しい解は得られない。h=0.1 で解いた数値解を表 3.6 に示す。最初の１ステップで 0.447 となるべきものが 112 にもなり、その後は殆ど変化しない。実際には、これから後は x0=0.1, y0=111.804 の解を求めることに変わってしまったのである。その解は、 Eq3_16d であり、表 3.6 の数値解はこの真値に良く一致している。
. . 最初のステップの相対誤差が著しく大となる理由は台形法の平均勾配(f0+f1)/2 が平均値の定理による勾配と著しく異なるからであり、a→0 では f→∞ であるから台形法の平均勾配も無限大となるので、x0=0, y0=0 からは f0 を用いる限り解けない。従って、このとき解が得られない理由は Lipschitz 条件とは関係なく、１ステップの相対誤差が無限大になってしまうからである。区間[0, h]でも平均値の定理を満たす勾配は存在し、それが求まれば x1=h における解は求まる。この微分方程式の解の x の小さいところでの平均値の定理を満たす勾配は区分点の勾配の平均ではなく、y を独立変数にすれば dx/dy=y となり解は x=y2/2 であることから明らかなように、区分点の勾配の逆数の平均の逆数である。両者の差が著しく大きいので台形法の精度は著しく悪いのである。x0=0 では f0→∞ であるから後者の平均は 2f1 となり、y1=2f1h は f1 が正しければ正しい解である。しかし、正しい f1 を得るために y1 に適当な値を与えて反復修正すると２安定となり収束しない。この場合は y1 を仮定して後退 Euler 法で y0 を求め、それが零になるように y1 を修正しているのだから、x0=0 で Lipschitz 条件を満たさなければならない。
Eq3_16e
この２安定の幅は c が大きければ広く、c が真の y1 に近ければ狭くなり、真値を挟んでいるので中間点を近似値としてこの操作を繰り返せば後退 Euler 法ではなくなるので正しい値に収束する。更に簡単に求めるには、２つの安定値が等しくなる c を求めれば正しい解 y12=2h が求まる。この後は既に述べたように、x0=h, y0=(2h)1/2 で解けるので、この微分方程式は数値解法で解ける。
. . 一般に、dy/dx=1/yn, y0=0, (n>0)のときも同様にして概略の y1 を求めることができるが、正確な値を求めるには平均勾配は、f0 の逆数の n 倍と f1 の逆数との和掛ける 1/(n+1)の逆数、即ち、(n+1)f1 としなければならない。この平均勾配も次式に示すように f1 が正しければ平均値の定理の勾配と正確に一致する。
Eq3_16f
. . y1 を任意に与えた反復修正は(n+1)f(y)h=y の浮動点を求めることであるが、縮小写像になっていないので収束しない。しかし、連続した２回の反復修正値について先のy1 の n 乗と次の y1 の積を y1 の n+1 乗と見なすと正しい解を得ることができる。
Eq3_17
これは連続した２回の反復値が収束し、同じ値になったと仮定したことと同じであり、連続した２回の反復値を等しいとしても求まる。簡単には、最初の仮定 y1=c と次の y1 が等しいとすれば(3.17)の結論の y1 が求まる。
. . この解法は微分方程式の求積解が分かっているから平均値の定理の勾配が分かったのであり、x=0 からの数値解法として適当ではない。そこで、x0=0 における f0 を使用しないために後退 Euler 法の反復修正により y1 の概略値を求めると、(3.18)の浮動点を求めることになるが、次式に示すようにこれも縮小写像にならないので浮動点を求めることは出来ない。
Eq3_18_______(3.18) Eq3_18a
しかし、方程式(3.18)の一つの解は y=n+1√h であり、これは真の解(3.17)の近似値として満たすべき精度に少し足りない程度である。例えば、n=1 の場合には(3.16)の精度があればよいが、y=√h は有効桁１桁の精度である。これは微分方程式と後退 Euler 法という条件だけから求まる近似値であるから初期値として使うことができる。従って、これ等の微分方程式を x0=0, y0=0 から数値解法により解くには最初の区間[0, h]については h を更に細分し、例えば、h1=h×10−4 位にして x0=y0n+1=h1 として解く。台形法では x=h 以降も細分しないと得られた精度を保てないのでステップ数が多すぎるのが難点であるが、著者の解法では高次解法で解き、十分な精度が得られる所迄きたら自動的に区分幅を 2h1, 4h1, …… と広くすることにより少ないステップ数で十分な精度の解を得ることができる。
. . y' =1/y, y(0)=0 を第２章のプログラム1(出力は 80 行と 28 行の必要なものだけに変更)で解いた結果を表 3.7 に示す。x=0～1 を 10 等分した点の解を求めるため、区間 0.1 を８等分した最初の点 x0=0.0125 における y0=0.158114 を第２章 2.2 節に従って入力する。E は表の１区間 h を更に E 等分した E 次近似で収束することを示す。低次の固定次数、固定区分幅解法ではこれだけの精度は達成できない。
表 3.7 y' =1/y, y(0)=0 (x0=0.0125, y0=0.158114, h=0.0125～0.05)
IBExy真　　値IBExy真　　値
1






2







3





 7
6
5
5
4
4
4
4
4
4
4
4
4
4
4
4
4
3
3
4
4		 0.025
0.0375
0.05
0.0625
0.075
0.0875
0.1
0.1125
0.125
0.1375
0.15
0.1625
0.175
0.1875
0.2
0.2125
0.225
0.2375
0.25
0.275
0.3		 0.2236069
0.2738614
0.3162279
0.3535535
0.3872984
0.4183301
0.4472137
0.4743418
0.5000001
0.5244045
0.5477226
0.5700877
0.5916080
0.6123725
0.6324555
0.6519203
0.6708204
0.6892024
0.7071067
0.7416198
0.7745966		 0.22360679
0.27386128
0.31622776
0.35355339
0.38729834
0.41833001
0.44721359
0.47434165
0.50000000
0.52440443
0.54772256
0.57008772
0.59160799
0.61237243
0.63245553
0.65192025
0.67082040
0.68920245
0.70710678
0.74161985
0.77459668		 4



5



6


7

8

9

10



 4
4
4
4
4
4
4
4
3
3
4
4
4
4
4
4
4
3
3


 0.325
0.35
0.375
0.4
0.425
0.45
0.475
0.5
0.525
0.55
0.6
0.65
0.7
0.75
0.8
0.85
0.9
0.95
1.0


 0.8062257
0.8366600
0.8600254
0.8944272
0.9219545
0.9486833
0.9746795
0.9999999
1.024695
1.048809
1.095445
1.140175
1.183216
1.224745
1.264911
1.303840
1.341641
1.378405
1.414213


 0.80622579
0.83666005
0.86002543
0.89442719
0.92195445
0.94868331
0.97467945
1.00000000
1.02469505
1.04880886
1.09544514
1.14017546
1.18321595
1.22474487
1.26491107
1.30384050
1.34164081
1.37840491
1.41421356


目　次　へ