拡散方程式をDormand-Prince(ルンゲクッタ8次)で計算(Excel VBAソース付き)
久々に偏微分方程式をExcel VBAで計算するシリーズ。まずは拡散方程式
∂u/∂t = D ∂^2u / ∂x^2
だ。
まずは空間部分を差分化して、
du(n, t)/dt = D * (u(n+1,t) + u(n-1,t) - 2*u(n,t)) / (Δx*Δx)
のような常微分方程式をそのままDormand&Prince(14段8次ルンゲクッタ法)で計算してみた。
さて、こういう常微分方程式をオイラー法で解くと、
D * Δt /(Δx*Δx) < 0.5
が必要で、結構刻み幅に気を使う。しかし、このDormand-Princeならば
1.2くらいにまで大きくしても全然安定。こりゃすごいな。
計算結果はこんな感じ。
ソースはこんな感じ(実用性はないけど)。
« RFワールドNo.11を買った。 | トップページ | 高速道路無料化-利用証明書はちゃんと受け取ろう。 »
「学問・資格」カテゴリの記事
- 高周波(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)
こんにちは。数値計算を興味深く読んでいます。
Dormand-Prince8次と書かれていますが、Fehlberg7(8)次と思っていいのでしょうか?
投稿: | 2010年10月 2日 (土) 01時17分
おっとそこに気付く人がいらっしゃるとは、、、
正式なDormand/Prince8次の係数はいわゆるdop853に記載されていますので
そちらをご参照ください。リンクはこちらです。
http://www.unige.ch/~hairer/prog/nonstiff/dop853.f
うちのは、、、途中で出典が違うことに気がついたけどめんどくさくなって
そのまま放置しているという、、、
投稿: tonagai | 2010年10月 2日 (土) 05時21分
こんばんは。早速のコメントありがとうございます。
次数の違いが気になってコメントさせて頂きました。
今後の数値計算も楽しみにしています。
投稿: | 2010年10月 2日 (土) 22時54分