エクセルを用いたルンゲクッタ法による微分方程式の計算

エクセルを用いたルンゲクッタ法による微分方程式の計算

Excelを用いた科学技術計算が第2版になりました　30年10月！

amazon：Excelで操る! ここまでできる科学技術計算第2版

楽天： Excelで操る! ここまでできる科学技術計算 [ 神足　史人 ]

　　　　　　　　　　　　　　図１　ルンゲクッタ法で解いた例題のグラフ

　解析的に微分方程式を解くのは非常に煩雑で、また、必ず解けるとはかぎらない。
微分方程式を１階連立微分方程式に変換し、数値積分法の一種であるルンゲクッタ法
を適用すると、ワンパターンの計算で数値解が求まります。

４次のルンゲクッタ法

　次の１階２元連立微分方程式で説明する。

　dφ/dt = F( t , φ, ω)
　dω/dt = G( t , φ, ω)

　F,Gは時間tと変数φ、ωの関数

初期値　t=t0　の時、φ=φ0, ω=ω0とすると、Δt後の変数値は下式となる。

　　k1 = Δt・F( t0, φ0, ω0)
　　m1 = Δt・G( t0, φ0, ω0)

　 k2 = Δt・F( t0+Δt/2, φ0+k1/2, ω0+m1/2)
　　m2 = Δt・G( t0+Δt/2, φ0+k1/2, ω0+m1/2)

　　k3 = Δt・F( t0+Δt/2, φ0+k2/2, ω0+m2/2)
　　m3 = Δt・G( t0+Δt/2, φ0+k2/2, ω0+m2/2)

　　k4 = Δt・F( t0+Δt, φ0+k3, ω0+m3)
　　m4 = Δt・G( t0+Δt, φ0+k3, ω0+m3)

　　k = (k1+2・K2+2・k3+k4)/6
　　m = (m1+2・m2+2・m3+m4)/6

Δt後の変数値φ1,ω1は

　　φ1 = φ0 + k
　　ω1 = ω0 + m

となる。

この手続きを繰り返せば、離散的だが、すべての時間での変数値が求まります。
ちなみに、Δtは0.0001秒のような微小値を設定する必要があります。
Δtが大きすぎると解が発散する場合があります。

計算量が多いので、手計算は困難です。EXCELなどを使って計算しましょう。

数値積分法の基礎と応用
参考文献：日本機械学会編

技術論文作成のための機械分野キーワード100解説集

独修微分積分学改訂増補

数値解析

たった1日で即戦力になるExcelの教科書

入門者のExcel VBA―初めての人にベストな学び方 (ブルーバックス)

Excel VBAのプログラミングのツボとコツがゼッタイにわかる本

エクセルを用いたルンゲクッタ法

　ルンゲクッタ法を用い(1)の微分方程式を解いて見ましょう。

　１階微分方程式

　　　dy/dx = f(x,y) 　f(x,y) = -y -10x + 5 　　　・・・・・　(1)

　　　初期値　x0 = 0 , y0 = 1
　　　刻み巾　h = 0.1

x0がh増加すると、ｙ0はｋ増加する。よって、

　( x1, y1) = (x0 + h, y0 + k ) となる。

ｋは下式で計算する。

k1 = h・f(x0, y0)
k2 = h・f(x0 + h/2, y0 + k1/2)
k3 = h・f(x0 + h/2, y0 + k2/2)
k4 = h・f(x0 + h, y0 + k3)
k = (k1 + 2k2 + 2k3 + k4 )/6

まず図２のように、エクセルで初期値( x0, y0) からｋを求めます。

　　図２　微分方程式の初期値からkを求めるエクセルシート

図２の青色のセルが初期値（x0,y0)の領域で、初期値と刻み巾より、k1,k2,k3,k4,kを計算します。
この初期値領域に(x1,y1) =(x0+h, y0+k)を設定すると、y1の増分kが計算できます。

それでは、図３のように計算結果のエリアに初期値を設定しましょう。

　　　　　　　　　　　　　　　　図3　計算結果エリアに初期値を設定したエクセルシート

計算の繰り返し処理のため、ＶＢＡをつくります。

　　　　　　　　　　　　　　　　図４　ルンゲクッタ法のための計算繰り返し用VBA

計算実行ボタンを作り、VBAをマクロとして登録し、VBAマクロを実行します。

　　　　　　　　　　　　　　　　　　　　図５　VBAマクロ実行結果のエクセルシート

計算結果を図６のグラフにしました。

　　　　　　　　　　　　　図６　ルンゲクッタ法による微分方程式の計算結果グラフ

参考文献
数値積分法の基礎と応用
日本機械学会編

Runge-Kutta-Fehlberg法

　新潟大学工学部化学システム工学科からの依頼により、エクセルを用いたRunge-Kutta-Fehlberg法を製作しました。Runge-Kutta-Fehlberg法は誤差をチェックし、時間刻み幅を変更しながら計算するため、通常のルンゲクッタ法と比較して非常に精度が良い方法です。通常のルンゲクッタ法では発散してしまう温度境界層問題も解けます。Runge-Kutta-Fehlberg法によるVBAは化学工学(反応速度式などに適応)のみならず機械工学にもたいへん有用と思われます。この常微分方程式解法シート(ITOシート)は新潟大学の化学工学資料のページからダウンロードできます。

参考　：　化学工学資料のページ 　　Excelで気軽に化学工学

エクセルを用いたルンゲクッタ法による多質点系バネマスモデルの振動解析

上記振動解析において、４次のルンゲクッタ法とRunge-Kutta-Fehlberg法を比較しました。Runge-Kutta-Fehlberg法はエネルギー保存則(運動エネルギー＋バネのエネルギー)を満たしており、ひじょうに精度の良い解法であることがわかる。また４次のルンゲクッタ法より精度が良いにもかかわらずその演算時間が格段に小さく、性能が優れていることが解る。

エクセルを用いたルンゲクッタフェールベルグ法によるクロソイド曲線(緩和曲線)

計算が難しいといわれているクロソイド曲線も容易に算出できる。

ルンゲクッタフェールベルグ法により電気回路方程式を解く

電気回路にも適用できる。

Excelを用いてラグランジュポイントにおける物体の運動方程式を解く
(決して宇宙人には教えないで下さい)

ラグランジュポイントでは物体は奇妙な運動を行う。

ルンゲクッタ法を用いた単振り子