4次R-K法で微分方程式を解く。数値解なのでアルゴリズムはわかっていてコードを書く。コピペでもいい。
dy/dx = 2π・cos(2π・x) を解く。解析階は y = sin(2π・x)である。よって数値解と解析階との差を出してみよう。
刻み幅 = 0.001の場合である。誤差が10のマイナス12乗より小さいので、まーそれはそれでいいのだが、言いたいことはそこではない。
C言語で計算してデータをファイル化、それをgnuplotで表示した。
例題は横軸10ぐらいで⊿t = 0.1 として実行。ほとんど瞬時で表示する。
ところがデータ数が増えると・・・
plot ‘data.txt’のdata.txtが8MBぐらいになるとgnuplotはしばらく応答なしになった。
そこだ。この手法では描画をするのに適していない。簡単な、つまりはデータ数が少ない時の確認用としては活用できる。データ数が膨大になったら この手法はよくない。最良ではない。
そこで思い出すのがwxMaximaであり、LTSpiceである。大量のデータがあるのに、描画が速い。そこに工夫がある。上のgnuplotは大量のデータを地道に真面目に表示している。細工はないのだ。
wxMaxima, LTSpiceに限らずだが、ここらの工夫は描画でのデータ間引きだろうと想像する。上の場合は⊿t = 0.001なので、gnuplotは0.001ずつのデータを地道に順番に描画するのだ。
それは元のブログラムとは関係ないからなあ。
そこだで。簡易的に描画はデータをファイルに落としてgnuplotで見ればいい・・・・と思っていた。ただ、それはデータ数が少ない時に有効であって。
データファイルが上の場合には8MBぐらいだが、それでも時間がかかった。20秒とか50秒とかの・・・ PCがポンコツなんでごめんね(笑)
数値計算のツールは微分方程式の数値解を求めてグラフ化するであろう。その場合にどのように工夫して描画を高速化しているか。間引くのだろなあと想像するけどね。
上の場合。横軸で⊿t = 0.001り意味はない。ベターと塗り潰されてしまうわ。なので間引いてもよい。その、間引くアルゴリズムってものがありそうだな。それは各社が工夫しているのだろうと想像する。アプリによって違うアルゴリズムだろうかもかもだ。
dy/dx = 2π・cos(2π・x) を解く。解析階は y = sin(2π・x)である。その差を上のグラフは示している。10のマイナス12乗のオーダーだから問題はない。。・・。という考え方はできる。だが、そこで段差がいくつかある。それはとのようにして生ずるのか? そこを考えて理解できれば僕の数値計算に対する理解は相当に深まると信ずる。
G++でdoubleの計算だ。4次R-K法だ。コードはパクリだ(笑)
上の場合に、さらに計算を続けると発散するのか? 発散するとして どうして? だ。誤差の蓄積が原理的にある・・・のは知っている。だが、それがどうして階段状の誤差として現れるのだ?
こういうものは考えると疑問は拡大する。doubleは60bitの倍精度の浮動小数点の表現だしなあ・・・・から始まる。4次R-K法事態は単純な数式だ。が、それがどうして、このような誤差の結果を生むのか? うむうむ。そこだな。生むのか・・・ うむうむ・・・ 笑えよな
入門レベルの例題をサイト等から拾って来てVScodeで動かして、あ、動いた、・・・終わり。
というスタイルでは学習としては効率はちょ~低い(笑)
で、その、ちょ~低いレベルの例題はWebに膨大にある(笑) どうでもいいようなものは大量にある・・・インターネットはゴミの山である・・・・と 毎日、ゴミブログを書いているワシが言いましたよってに。