陰的ルンゲクッタ法(ラダウIIA)をExcel VBAに移植した。
硬い(スティフな)微分方程式を数値的に効率的に計算するには、陰的ルンゲクッタ法がいいということはよく知られている。しかーーーーし、これをVBAで計算している例を見たことがない。普通のCやFortranだって結構少ない。
このブログでは14段8次の陽的ルンゲクッタ Dormand-Prince法や、35段14次のルンゲクッタなど、普通のネットでいくらでも見られる4段4次ルンゲクッタなどとは一線を画す(うそです。単に私の暇つぶし。。。)方法でいろいろな微分方程式を計算してきたが、さらに進んで陰的ルンゲクッタに進もうと。
で、先日、連立一次方程式をLU分解をExcel VBAで計算するライブラリを移植したのでできるはず。できるはずだが、結構めんどくさかった。
陰的に計算するのに多変数のニュートン・ラフソン法を使う。なので、ヤコビアンも手で与える必要がある。方法は「常微分方程式の数値解法II」に出てきた「ラダウIIA」(Radau IIA)法を使った。3段5次でA安定。この本にはラダウIIAとか、ロバットIIICとか、マジンガーZの機械獣のような名前の手法(古すぎる。。。)がいっぱい。
で、ようやく完成。プログラムは無茶苦茶だが、参考になるかと一応ソースごと公開。実用的に使うにはちゃんと収束判定とかしてね。
まずは硬くないけど、ローレンツ方程式。
ふつーに計算できた。さて次は硬い方程式に進もう(続く)。
« Excel VBAにNumerical RecipesのLU分解ライブラリを移植して連立方程式を計算 | トップページ | 「無限記憶」を読んだ。 »
「学問・資格」カテゴリの記事
- 高周波(RF・マイクロ波・ミリ波・5G)関連ニュース2021年2月16日 IEEE Microwave Magazineの特集はオールデジタルのRFID、Microwave JournalはEバンド ミリ波通信に衛星や気球を使う話、アメリカの半導体企業がバイデンに投資を迫る、(2021.02.17)
- カオスを生じる電気回路、Chua’s circuitをLTspiceで回路シミュレーションしてみる。(2021.02.19)
- Labyrinth Chaos(迷宮カオス)を生むThomas-Rössler方程式のパラメータbを色々変えて、Python+Scipyでルンゲクッタ8次のDOP853(Dormand&Prince)を使って計算してGIFアニメ(2021.02.16)
- フィッツヒュー・南雲 (FitzHugh-Nagumo) 方程式をPython+Scipyでルンゲクッタ8次のDOP853(Dormand Prince)で計算。(2021.02.23)
- 「水晶振動子の等価回路計算」をカシオの高精度計算サイトkeisan.casio.jpの自作式としてUP! インピーダンスの大きさと位相がグラフ化できる。(2021.02.12)
« Excel VBAにNumerical RecipesのLU分解ライブラリを移植して連立方程式を計算 | トップページ | 「無限記憶」を読んだ。 »
コメント