こんばんは
またまた小難しいテーマになりますが、
興味のあるかたはお読みくださいませ。
解析的に解けないが、数値積分すれば解が求められるという
非線形微分方程式の解法には昔から興味がありました。
10年くらい前に、一度、表計算ソフトのエクセルを使って
解いてみようとトライしたことがあります。
計算式をチマチマとセルに入力して、コピペを繰り返して、
なんとか解らしきものはグラフの形になったのですが、
解の質＝有効桁数・・・
要するに、
この答えはどれくらいの桁まで正しいのかわからない・・・
一桁？、二桁？、・・・・・？
とりあえず計算の条件を変えて違いを調べるとなにかわかるかも？
ということでやり始めたところ、
エクセルの使用セルの量が爆発的に増えて、
ファイルの容量も数100メガバイトに膨らんでしまい、
開くだけでも10分近くかかるような状態でした。
エクセルは、セルの中の計算式はファイルを開いたときに
実行するのですね。
計算式の量が少ないときには開いたとたんに終了するのですが、
（ふつうは問題ない時間です、）
その量が膨大になると開くだけでも時間が掛かってしまいます。
そんなこんなで、中途半端なまま、自主的に終了していました。
まあ、誰に言われたことでもなく、
興味本位で始めたことなので、大勢に影響はないのですが・・・・・
それでもときどきファイルを開いては、
これなんとかしたいなあ・・・と思いつつ、
気が付けば10年ほど月日が経っていました。 (´・ω・`)
最終的には非線形微分方程式を解くのが目標ですが、
まずは厳密解がわかｔっている簡単な線形の方程式を解いて
その性質、
とくに推定残差（厳密解との差）を評価する計算式
（らしきものでもいいので）が見つかれば、
それを使って解析的な厳密解の残差を推定したいなあ
と思っています。
さて、数値積分の計算結果は一般に次の式になります。
計算結果の値 ＝ 真値 ＋ 残差
になるので、変形すると、
残差 ＝ 計算結果の値 ー 真値
両辺を真値で割って。
残差 / 真値＝ （計算結果の値 ー 真値） / 真値
になります。
具体的な数字では、
例えば
計算結果の値＝１１、真値＝１０ とすると
残差 / 真値＝ （計算結果の値 ー 真値） / 真値
＝ （１１ ー １０） / １０
＝ ０．１
です。
ここで、今、計算結果の値だけはわかっていて、
もし、
何らかの方法で 「残差 / 真値」 を知ることができれば
計算結果の値 ｘ 「残差 / 真値」
を計算すればおおよその残差を推定できます。
上の例では
推定残差
＝ 計算結果の値 ｘ 「残差 / 真値」
＝ １１ ｘ ０．１ ＝ １．１
です。
推定残差にも±の幅があると考えて、
計算結果は １１±１．１ 程度の範囲内にあることになります。
かなりザックリ言えば、
二けた目が怪しいから信用できるのは一けただけ
ということですね。
この 「残差 / 真値」に相当する式を探していきます。
実際に数値積分に使うツールは表計算ソフトのエクセルです。
とは言っても、セルに計算式を埋め込む普通の方法では、
あまりにも作業時間やファイルの容量が増えてしまうので、
やってられない・・・
これ過去に経験済みということで、
今回はエクセルのマクロ機能を使ってプログラムを組みました。
エクセルのマクロは2000年代くらいから10年ほど使っていましたが、
ここ10年くらいはとんとご無沙汰しています。
それで、復習の意味もあって、改めてプログラミングです。
とは言っても（この言葉多いですね。）、
セロから始めたわけではなく世の中に出回っている
数値積分のサンプルプログラムをベースに
アレンジしてみました。
前置きが長くなりましたね。
では、ここからが本筋です。
まず、
サンプルプログラムにも入っている1階の線形微分方程式が
こちらです。
dy/dx=-２ｘｙ
この微分方程式は変数分離すると簡単に解けてしまいます。
式を変形して
dy/ｙ＝ー２ｘｄｘ
これを積分して
∫（１/ｙ）ｄｙ＝-∫（２ｘ）ｄｘ
∴ loge(y)＝-ｘ＾２ ＋ｃ1
（ｃ1は任意の定数）
さらに変形して
ｙ＝ｃ・exp（-ｘ＾２）
（ｃは任意の定数で、ｃ＝exp(c1) )
ですね。
このように方程式の厳密解がわかっていますから、
数値積分で解いた結果と比較できますね。
積分にはサンプルプログラムに使っている
ルンゲクッタ法で行いました
計算するにあたって、決めないといけないのは、
ｘとｙの初期値とｘの範囲、それにｘの刻み幅ｈですね。
特に刻み幅ｈの値は重要です。
ｈを決めれば、淡々と計算すればそれで終わりですが、
今回はｈをあるひとつの値（これをｈの初期値とします。）にして
決めたｘの範囲の各々のｙを計算します。
次にｈを半分にして同じように計算し、さらにその半分・・・
というように複数回同じことを繰り返します。
つまり、刻み幅ｈを
ｈ、h/2、h/4、・・・、ｈ/（２のべき乗） のように繰り返すわけです。
ケースバイケースですが、
今回は、1/（２の10乗）から1/(2の2０乗）くらいで
様子をみながら試行錯誤でやってみました。
同じ計算範囲で刻み幅ｈを小さくするので、
当然計算量は増えます。
ｈを半分にすれば2倍、・・・・1/1024にすれば１０２４倍・・・
1/(2の2０乗）にすればとてつもなく膨大な量になります。
まあ、私が計算するわけでもなく、
PCの解散が終わるまで
数時間程度待てばいいだけのことですが。
マクロ内で定義した繰り返し数の設定を変えるだけで済むように
プログラムを作ったので、
デバッグするときはこれを小さくしたりと
いろいろ工夫したつもりです。
同じ処理（計算）が何度でも簡単にできるのも
マクロの強みですね。
ただ、長時間（例えば数時間）マクロを走らせっぱなしだと
windowsから応答なしのメッセージが返ってきて
気持ち悪い状態になります。
マクロ自体は動作しているので問題はないのですが・・・
これを回避するためにマク処理のきりのいいところに
「Doevents」コマンドを入れて、
一旦windousに返すという処理を入れたりしています。
さて、出てきた計算結果の料理方法です。
刻み幅ｈを設定して変数ｘを変えながらｙを計算していきます。
ひとつのｈについて計算が終わると
さらにｈを半分にして同じことを繰り返します。
このとき、
真値をy_ture=f(x)、積分結果をy_cal=g(x,h)とおくと、
積分結果は真値と残差ε(x,h)の和になるので
y_cal = y_ture + ε(x,h)
または
g(x,h) ＝ f(x) + ε(x,h)
と表すことができます。
残差ε(x,h)は未知数で、
真値f(x)も不明ならばこれも未知数になります。
刻み幅ｈを変えたときの計算結果について考えると
変数ｘが同じ値ならば
真値は変わらず残差だけが異なることになります。
つまり
x=x1、h=h1及びh=h2(h1≠h2）とおくと
g(x1,h1) ＝ f(x1) + ε(x1,h1)
及び
g(x1,h2) ＝ f(x1) + ε(x1,h2)
なので、辺々の差を取れば
g(x1,h1) - g(x1,h2) ＝ ε(x1,h1) - ε(x1,h2)
となります。
ただ、このままでは「残差 / 真値」の形にならないので、
若干発想を転換して、式を変形します。
y_cal = y_ture + ε(x,h)
= y_ture ・ ｛ 1 + ε(x,h) / y_ture ｝
または
g(x,h) = f(x) + ε(x,h)
= f(x) ・｛ 1 + ε(x,h) / f(x) ｝
と変形しておいて、 先に計算した g(x1,h1) - g(x1,h2) を
y_cal またはg(x,h) で割り算します。
すると、一応「残差 / 真値」の形が出てきます。
まず、g(x1,h1) で割ると
｛g(x1,h1) - g(x1,h2)｝ / ｛ g(x1,h1) ]
= { ε(x1,h1) - ε(x1,h2) }
/ [ f(x1) ・｛ 1 + ε(x1,h1) / f(x1) ｝ ]
= { ε(x1,h1) - ε(x1,h2) } / f(x1)
/ ｛ 1 + ε(x1,h1) / f(x1) ]
≒ { ε(x1,h1) - ε(x1,h2) } / f(x1)
・ ｛ 1 - ε(x1,h1) / f(x1) }
ただし、
ε(x1,h1) / f(x1) << 1 が成立するものとします。
次に g(x1,h2) で割ると
｛g(x1,h1) - g(x1,h2)｝ / ｛ g(x1,h2) ]
= { ε(x1,h1) - ε(x1,h2) }
/ [ f(x1) ・｛ 1 + ε(x1,h2) / f(x1) ｝ ]
= { ε(x1,h1) - ε(x1,h2) } / f(x1)
/ ｛ 1 + ε(x1,h2) / f(x1) ]
≒ { ε(x1,h1) - ε(x1,h2) } / f(x1)
・ ｛ 1 - ε(x1,h2) / f(x1) }
ただし、
ε(x1,h２) / f(x1) << 1 が成立するものとします。
二つの式を辺々引き算すると
左辺 = ｛g(x1,h1) - g(x1,h2)｝
・ [ ｛１ / g(x1,h1) ｝ - ｛１ / g(x1,h2) ｝ ]
右辺 = { ε(x1,h1) - ε(x1,h2) } / f(x1)
・ [ 1 - ε(x1,h1)/f(x1) - 1 + ε(x1,h2)/f(x1) ]
= -1 ・ ｛ { ε(x1,h1) - ε(x1,h2) } / f(x1) ｝＾２
=-1・｛ （ふたつの刻み幅の残差の差/ 真値）の2乗 }
となります。
右辺のマイナスは両辺に-1を掛けて正の値にしておきます。
左辺は
刻み幅の異なるふたつの数値積分した計算結果から求められ
これが「残差 / 真値」の形の2乗になることがわかります。
したがって。左辺のルート
（これを評価式と呼ぶことにします）
を計算すれば
「残差 / 真値」の形となります。
これを「推定残差（比率）」と呼ぶことにします。
問題は
この「ふたつの刻み幅の残差の差 { ε(x1,h1) - ε(x1,h2) } 」
と現実の残差の関係です。
直観的には何らかの関係がありそうに思えますが
この挙動が見えないので、
実際に計算して、「真値残差」と比較していきます。
ではまた
部分方程式 dy/dx=-２ｘｙ に話題を戻します。
まず、図1がdy/dx=-２ｘｙの厳密解のグラフです。
初期値はｘ＝０のときｙ＝１としました。
ｙ＝exp(-x^2)をｘが0から10までプロットしました。
数値積分もこの範囲で行います。
普通に見るグラフですね。
図 １
このままでは縦軸ｙの範囲が狭くて見ずらいので、
図2のように縦軸を対数表示に変えて
さらに数値積分した結果をふたつ重ねて比較します。
ｘが４より小さいときは3本のカーブはほぼ重なっていますが、
ｘがそれより大きくなると次第にずれが大きくなりますね。
図 ２
図2の丸印辺りを拡大したものが図3です。
拡大すると文字通り桁違いの差がみえてきます。
図 ３
もっとも、その大きさ自体十分小さくてほぼゼロになりますが・・・
では、数値積分した結果について。
まず、相対刻み幅と真値残差（比率）の関係が図４です。
ここで、横軸の刻み幅は初回を0.1として
これを１とした相対値で示しています、
右方向で刻み幅が小さくなります。
縦軸は真値残差（比率）を対数表示で表しています。
上方向で残差が大きくなります。
なお、青いカーブは残差の最大値を、
オレンジ色は最小値を示しています。
なお、真値残差（比率）は絶対値としています。
図 ４
直観的には刻み幅が小さくなれば残差も小さくなるはずで、
図4の赤色と緑色の矢印の刻み幅の位置までは
確かにそうなっています。
が、
さらに刻み幅が小さくなるとそれ以後は逆に増加に転じます。
これは数値積分が本質的に積算処理の連続であることに
起因していて、
積算回数は刻み幅に反比例して増加し、
わずかな残差が回数分積算された結果です。
刻み幅は
やみくもに小さくすればいいものではないということです。
まさに「塵も積もれば山となる」ですね。
次に、
評価式を使った推定残差(比率）と相対刻み幅の関係が図５です。
図４と同じ条件にしています。
二つの図を比べると同じ傾向なのがわかります。
図 ５
図6が刻み幅をパラメータにして
変数ｘと推定残差（比率）の関係です。
横軸がｘで右方向で増加します。
また、縦軸は推定残差（比率）で対数表示しています。
パラメータは刻み幅で下方向で小さくなります。
図 ６
変数ｘが増加すると推定残差（比率）も増加していて、
数値積分の積算の結果ですね。
図７は、
図5に重ねて推定残差(比率）の漸近カーブを加えたもので、
緑色と水色の破線がそれぞれ減衰方向の最大値、最小値の
漸近カーブになります。
水色の矢印方向の傾きは刻み幅の4乗です。
また紫色とオレンジ色の破線は増加方向の漸近カーブで
色の矢印方向の傾きは刻み幅の－１乗になっています。
推定残差（比率）の最大値、最小値の極小点は
それぞれ赤色と緑色の矢印で示しています。
図 ７
では、図４と図5を重ねてみましょう。
図８になります。
対数表示ですが、ほぼ一致しているといえますね。
図 ８
次に、
変数ｘに対する二つの残差（比率）の違いを調べてみました。
刻み幅は図８の赤色の矢印に対応した値としました。
図９では
真値残差(比率）は絶対値ではないので、±の極性があります。
図 ９
ｘの増加とともに単調増加とおもいきや、
どちらも増加したり減少したりやや複雑な動きをしています。
次に、
図９の絶対値をとって、
最大値と最小値にマークしたものが図１０です。
赤い丸印が最大値、緑色の丸印が最小値で、
オレンジ色の直線が推定残差（比率）の最大ラインです。
このラインより下なら問題無しです。
図 １０
図１１では、
横軸に真値残差（比率）、縦軸には推定残差（比率）として
相関図としてみました。
赤と緑の丸印は図１０と同じでオレンジ色も同じ意味です。
図 １１
というわけで
今回の例の結論は、
図８のように真値残差（比率）と評価式の相対残差(比率）が
ほぼ一致したというこどですね。
次は、
やはり厳密解がわかっている別の微分方程式について
調べてみます、
最後までお読みくださったかた、
大変お疲れ様でした。
＆ ありがとうございました。。