年末である。
職場の廊下に人気はない。（「にんき」ではない。そう私は任期制教員では・・・辞めておこう。）
雑事から逃れ、自由気ままに思考するのに最適な条件・・・であるはずなのだが、年内にしたかった、そしてすべき仕事が残っていて、思考に制約がかかってしまう。
時間もないので、辞めておけば良いのに、紙に標題のことに関する思考過程を書いていて、ふと思った。ここに書いておけば良いじゃないか。
推敲もせず、思ったまま書いていくので、おそらく整理されておらず、中途半端で終わるだろうが、ご容赦願いたい。
Runge-Kutta法、そしてRunge-Kutta-Fehlberg法で初期値問題を解く際に、従属変数（値を求めたいもの。例えば濃度）の値が0のものが変化し始めるときの誤差が大きく、独立変数（たとえば時間t）の変化幅が小さくなりすぎるという問題が生じる場合が多い。
その解決策を模索したい。
定常状態近似法のようなものを使うと良いのではないかと思った。
具体的に考察しよう。
たとえば、（式をあるアプリ風に書いて）
例）{ dn sub {i} } over dt = {k over 2} sum from {j=1} to {i-1} {n sub {j} n sub {i-j}} - k n sub {i} sum from {j=1} to {infinity} {n sub j}
のようなSmoluchowskiのRapid coagulation問題のようなことを考える場合を扱う。
iは１以上の整数。
独立変数は時間t。
t=0でn sub {1} = n sub{10}、i>1のすべてのiについて、n sub {i} = 0とする。
時間変化の幅Δt=hとして、
・・・なぁんて、紙に書いたものを写して書いていると、式を書くだけで疲れてきた。やはり、紙とペンの方が楽だ。
GoogleのJamboardにでも書こうか。
肝心のところにたどり着くのに、とても時間がかかる。
なので、途中は跳ばす。定常状態近似法についての考察を丸跳ばしにする。
Runge-Kutta方でのt=0からの最初のステップ(t=0→h)では、 n sub {i}が0ではないような最大の重合度iの値は、
1+1→2、2+2→4、4+4→8、8+8→16、16+16→32、32+32→64
で、k次について、2 ^ k なので、5次で32、6次で64になる。
したがって、Runge-Kutta-Fehlberg法の1ステップで、 ・・・
一部省略・・・
hが小さくなり、6次の項が桁落ち（C言語など）して0と誤って認識されるまでhが小さくなる。小さくなりすぎる。C言語でlong double型だと、10の-4000乗程度以下になるまでである。doubleだと、そんなに小さい数を扱えないので、返ってhの最小値が大きく、計算が早い。
でも、そんな、10の-4000乗程度の濃度の値って必要だろうか。
・・・
時間に余裕がなく、全部は書いてられません。
・・・
実際の現象では、有限の数の分子がランダムに反応（衝突・凝集でも等価）する。微分方程式が表す濃度n sub {i}およびその時間変化は、無限に大きい系、すなわち、含まれる分子の個数が無限大の状況での、濃度の値。
実際の系では、有限なので、モンテカルロで多くの分子を取り扱った方が、実際の結果に近いような気がする。
モンテカルロでは、多数の場合（分子と考えても良い）を取り扱い、その平均値が全体の平均値と等しいとするが、その多数とは、有限である。多くの分子を含む系から、ランダムに一つを取り出して、性質を乱数で発生させて決定する。これを同様に多くを取り出す（サンプル、抽出）する。
実際に実験で取り扱う系でも有限の個数の分子を含む。もちろん、１モル=6.02x10^23個の分子をモンテカルロシュミレーションで取り出す（計算する）のは大変長時間かかる。一生で計算が終わらないかもしれない。
ただし、現実は、無限に大きい数の分子が起こしているのではなく、有限である。
パラレルワールドというか、無限に起こる多くの可能性のうちの一つが現実である。
モンテカルロシミュレーションの結果は、同じ数の分子を取り出す場合でも、それを繰り返すと、値が異なる。特に、確率が低い事象で生成する分子、つまり低濃度の分子の濃度は、一回毎の、同じ値の大きな数の分子を取り出すシミュレーション毎に異なる。
これは、不確定性原理にも関係するような気がする。つまり、すべてを正確な値で測定・観測することは不可能で、ある程度の誤差が含まれる。
低確率の事象で生成する低濃度の分子が、発生するかどうかによって、系の状態が大きく異なる場合がある。そのような場合にカオスがクローズアップされるのだろうか。
ちょっと話がそれた。
Runge-Kutta法やRunge-Kutta-Fehlberg法で解いたとしても、特に短時間での重合度iの大きな分子の濃度は、間違いを含みやすい。上述のことが起こって、無限大の個数の分子を含む系では起こる事象が、起こらないことになってしまう。
ただ、それほど小さい確率の事象で生成する分子を正確に取り扱う必要は、上述の理由で、不要かもしれない。
ちょっと違った観点から考察すると、衝突も平均自由行程程度動ける程度の短い時間で考えると、衝突が起こらないという時間帯があっても良い。このような場合でも、微分方程式を解くと、衝突がわずかな程度に起こっていることになる。こうなると、むしろ、微分方程式を数学的に厳密に解くということは、現実に即していない、間違いのようにも思える。蛇足ながら、我々が暮らしている温度では気体中の窒素分子は、他の分子と衝突しなければ、一秒間におよそ500m程度の速度で運動している。平均自由行程はうろ覚えなので間違っているかもしれないが、一気圧なら分子直径の200倍程度だったような。一分子の直径をおよそ1nmとすると、平均自由行程は200nm程度。これを500m/sで割ると、平均的に一回の衝突を起こす時間間隔はおよそ10のマイナス9乗秒になる。もちろん、液体や固体中の分子の数密度は気体中よりもはるかに大きいのだが。
私がRunge-Kutta-Fehlberg法である系で衝突を計算させたら、long double型変数で計算させたところ、時間変化幅hが一時、
先の話に戻るが、C言語で扱える絶対値が最小の数は10のマイナス4000乗程度だが、これは、アボガドロ数の10の乗数、23の2乗以上である。その程度というのは、現実の系でも起こっていない、すなわち、反応などの速度が0とみなしても良いと思われる。
この0とみなしても良いだろう値、それをどう、C言語のプログラムに与えるか。
たまたま私が実行している環境なら、
double型なら10のマイナス307乗
long double型なら、10のマイナス4931乗
程度以下で、0とみなされている。double型で十分な気がしてきた。いや、もっと大きい値で0とみなしても充分だろう。
