数値計算を使って常微分方程式を解く~ルンゲクッタ法の解説~

資料請求番号:TS33 TS38

数値計算による微分方程式解法の基本~runge-kutta法~

ルンゲ・クッタ(runge-kutta)法は、コンピュータを使った微分方程式解法の基本として工学部の基礎教育に登場します。
中でも4次のルンゲ・クッタ法はアルゴリズムが簡単な上に高い精度を持つことから、実際の研究や企業におけるR&Dの現場にも利用されており、多くの工学者からの支持を受けています。

今回は、コンピュータを使った微分方程式解法(=数値計算)の基本的な考え方
非常に簡単なアルゴリズムの微分方程式解法であるオイラー法
そしてその進化系的存在である4次精度ルンゲ・クッタ法について解説します。

微分方程式をコンピュータで解くということ

ストーク

今からやろうとしていることは、こんな常微分方程式をコンピュータを使って解くということだ。

ストーク

微分方程式をコンピュータで解くにあたって、要となる考え方は?

シママ

方程式を離散化させること!コンピュータは基本、四則演算しかできないから、微分積分、微分方程式のような高等数学をコンピュータで取り扱いたいなら四則演算に落とし込む作業が必要なのよね。

ストーク

せやな。連続である微分方程式を、こんな風に離散化させて、四則演算にする。

ストーク

そしたら、前のyの値であるyiから、次のyi+1を求められるよな。
もちろん、初期条件x0y0が与えられていること前提やけどな。

シママ

hxの変化分よね。hを極限まで小さくすれば、微分の定義そのものになるわね!

ストーク

せやな。これを繰り返せば最終的に微分方程式の解y=f(x)を構成する点の集合が得られるたいね。

ストーク

本来の解が下の図Bのような場合、上の図Aの集合体でいかに下の図Bに近い形を描くか?これが数値計算の精度たい。

シママ

なるほど・・・。

シママ

・・・ところで、この中にあるφって何?

ストーク

このφの関数形が微分方程式の数値解の精度を握っとるたい。

オイラー法

ストーク

φの関数形がf(xi,yi)ならば、オイラー法になる。

シママ

あ、この式なら、解きたい方程式がdy/dx = f(x,y)の形であって、初期条件f(x0,y0)が与えられていれば、逐次的にyi+1を求められるわね。

ストーク

せやな。実際に図に起こすとこんな形になるたい。

シママ

ああ、そっか。xi,yiにおける青線の傾きをf(xi,yi)として、その傾きのままxi+1まで進んだとしたら、yi+1はどうなるのか計算するのがオイラー法なのね・・・。

シママ

・・・これじゃあ、凄く精度が悪いじゃない。

ストーク

せや。この方法を繰り返してたらもはや何を計算しとるのか分からんくなるたい。

シママ

hを小さくしたらどうにかなるかしら・・?

ストーク

どうにかならんこともないが、計算負荷が上がる。
しかも急激に値が変化するような関数や、曲率の大きい非線形な関数をこれで取り扱おうとするとなおさらたい。

シママ

そうね・・・しかも急激に値が変化するような関数は微分方程式の解の形として、よくあるわよね。例えば指数関数とか。
それをキチンと再現するためにhを小さくしていったら・・・。キリないわ。

ストーク

せやねん。そこで、f(xi,yi)だけじゃなくて、色々な値φを使って次のyi+1を予想しようってのが、ルンゲ・クッタ法たい。

ルンゲ・クッタ法

シママ

ワタシ、ルンゲ・クッタ法を実装して解いたことあるけど・・・。

シママ

この数式の通りにプログラムを書いただけで、実際のところ、ルンゲ・クッタ法のこと、よく分かってないのよね・・・。

ストーク

まぁ、確かにこの数式だけあれば、ルンゲ・クッタ法で微分方程式は解けるけどな。俺もExcel使って解いたことあるたい。

ストーク

でもな、この数式の意味をしっかり理解しながら計算を進めていくことができれば、
もしエラーが出たときにエラーの処理をしやすいし、計算結果を誤解することも防げるたいね。

シママ

計算結果の・・・誤解?

ストーク

数値解析の答えを実際の答えだと信じて疑わないことたい。
「数値解析ソフトやプログラムが出している結果なんやから、その答えは間違ってなくて現象もその通りになるはずや」と考える、一種の思考停止たい。
数値解析のアルゴリズムの意味や解析ソフトが解いている方程式の意味を分かっていないとこういう恥ずかしいミスをする。

シママ

(うっ、グサグサくる・・・。)

ストーク

そういう背景もあるから、ルンゲ・クッタ法の表している数式がどういう考えで作られているのか?考えてみる。
そして可視化、すなわちグラフにしてみたらどうなるのか知りたいという欲望を持つことはとても大事たいね。

シママ

まずね、(k1+2k2+2k3+k4)/6が分からないのよね・・・。
どういう考え方なのかしら・・・?

ストーク

これは平均しとるたい。係数の数字を足すと・・・・つまり1+2+2+1は6やろ?
1個分のk1,2個分のk2,2個分のk3,1個分のk4、全て足し上げて6で割れば平均になる。
特に、k2k3の影響が強く出るように設定しとるたい。こういうのを「重み付き平均」という。

シママ

kの種類が4つで6で割っているから平均と分かりにくいけど・・・言われてみれば、そうね。

シママ

と、いうことは・・・4種類のを求めて、それを平均したものを足してyi+1を求めているということになるのね。

ストーク

せや。単にφ=f(xi,yi)としていたオイラー法に比べたら確実にバージョンアップしとるやろ?

シママ

うん。

ストーク

それで4つのkの求め方って言うのが以下の通り。下の図と何度か見比べてしっかりイメージ掴むんや。

ルンゲクッタ法 手順

①点(xi,yi)でf(xi,yi)を計算して1回目の傾きとする(青色の直線)。これをxi+1まで伸ばすとyk1だけ変化する。

xiからh/2だけx方向へ、yiからk1/2だけy方向に進んだ点でf(xi+h/2,yi+k1/2)を計算し2回目の傾きとする(赤色の矢印)。この傾きで点(xi, yi)からx方向にhだけ進むとyk2だけ変化する。
※色の同じ直線は互いに平行

xiからh/2だけx方向へ、yiからk2/2だけy方向に進んだ点でf(xi+h/2,yi+k2/2)を計算し3回目の傾きとする(緑色の矢印)。この傾きで点(xi, yi)からx方向にhだけ進むとyはk3だけ変化する。

xiからhだけx方向へ、yiからk3だけy方向に進んだ点でf(xi+h,yi+k3)を計算し4回目の傾きとする(橙色の矢印)。この傾きで点(xi, yi)からx方向にhだけ進むとyk4だけ変化する。

⑤上記4ステップで求めたk1, k2, k3, k4にそれぞれ1,2,2,1の重みを付けて平均した量を求め、点(xi,yi)からx方向にhだけ進んだ時のy方向変化量を求め、yi+1を決定する。


シママ

あ、これ、k1だけでyi+1を決定したら、オイラー法と同じになるのね。

ストーク

お、そこ気づいたか?せやせや。1点だけで次の点を決めるオイラー法とは違って、4つの点を使って、その平均で次の点を決める。しかも、2個目、3個目はxixi+1の中間で傾きをとるたい。こういった工夫によって、関数の急激な変化に強くなるたいね。

シママ

確かに、微分方程式の解曲線y = f(x)に近い値が予測出来てるわね。さっきのオイラー法とは大違い。

ストーク

まぁな。本当はもっと点が近いんだが、見やすくするため、少し距離置いとる。それくらい精度ええたい。

例題1 基本問題

ストーク

じゃあ、ルンゲ・クッタ法を理解したところで、実践たい。この問題をやってみよう。

放射線汚染の原因となる核種の一つにヨウ素131(131I)がある。131Iは放射線を出しながら安定状態であるキセノン131(131Xe)へ変化する。131Iの数をNとすると、壊変定数λを用いて、核の崩壊速度dN/dtは以下の式で表される。

131Iの壊変定数λは0.08664 day-1であることが知られている。t=0において100個の131Iが存在する場合、それが半分になるまでの時間を4次精度ルンゲ・クッタ法による数値計算により求めよ。

シママ

・・・あれ?この微分方程式、変数分離で解けるんじゃ・・・?

ストーク

ああ。解ける。解けるんやが、ルンゲ・クッタ法の練習として、一ついい例題かな思て。

シママ

何の言語でプログラム書こうかしら?

ストーク

これ、Excelで計算できるたい。せやからExcelでやろ。
R&Dの現場によっては、プログラミングに触れたことのない人に仕事の引継ぎをすることが多いし、
プログラミング言語のコンパイラをダウンロードしたりと色々と面倒な手続きが要る場合が多いから、Excelで計算することが多い。
実際、俺の部署では計算はExcelとProⅡとFLUENTを組み合わせて行う。

シママ

え~・・・。CやFortranが使えないの・・・・?

ストーク

プログラムというプログラムは「VBAマクロ」と言ったとこたいね。この問題はマクロ組むまでもない簡単なものだから、Excelの表計算で十分やと思うたい。

ストーク

それで、ルンゲ・クッタ法で解く微分方程式の基本形が左で今回解きたい方程式が右や。

シママ

yNxtに相当するのね。

ストーク

せや。そしたら、ルンゲ・クッタ法の式は次の様になる。

ストーク

そして、こんな感じでExcelシートを作っていくたい。

シママ

それぞれのkの値を計算する時、絶対参照のB2セルがhに相当して、その後のカッコの中で-λNを表現しているのね・・・。

シママ

あ、確かに、k2を見るとh*(-λ*(N+k1/2))を計算しているわね!

ストーク

せや。この微分方程式の右辺には時間項がないが、
それぞれのセルでN+k1/2, N+k2/2, N+k3がしっかり計算されているのが分かるやろ?
そしてE3セルでそれぞれのkの値の重み付き平均をとって次のNの値を求めるたい。

シママ

これでこのまま、セルを下まで引っ張っていけば、核の数が半分、つまり50になる時間を求められそうね!

ストーク

ああ。実際にE列の数字が50以下になるところをMATCH関数で探してみたところ、82番目のセルでギリ50以上、83番目のセルで50以下になることが分かった。
つまり、半減期はD82セルで返せるってことたい。
せやから、答えは8日たい。

シママ

なるほど・・・。これ、変数分離で解いてみて答えのチェックしてみてもいい?

ストーク

ああ。もちろん。

シママ

初期値のN0としたら、半減期って、こんな風に変数分離法で求めることができて・・・、e底のlog2が0.6931、壊変定数が0.08664だったから・・。

シママ

あ、ホントに半減期8日になった!

ストーク

せやせや。4次精度のルンゲ・クッタ法を使えば、かなり精度高く微分方程式をコンピュータで解くことができるたい。
今回は刻み幅hを0.1にして計算したが、刻み幅1や5でも精度高く解曲線を出せる。

シママ

え~ホントに?そんなに荒っぽくても大丈夫なの?

ストーク

問題の微分方程式を解いてy=f(x)にあたるN=N(t)を求めると以下の様になる。

ストーク

赤文字の数式がグラフ上の実線、ルンゲ・クッタ法で求めた解を青の点線で表すと以下の様になる

シママ

あ~スゴイ!ピッタリ!!

ストーク

ああ、でもな、今回は解曲線が単調減少という簡単な形だったから数値計算で再現しやすかったという考え方も必要たい。
解がどんな曲線になるのか想像しながら、精度と計算負荷両面から考えて絶妙な刻み幅を設定して計算するのもエンジニアとしての腕前たいね。

ストーク

それから、今回は変数分離できる簡単な微分方程式を扱ったが、本来、紙面上で解くのが不可能だったり、高度な数学を用いる必要がある場合にこの方法を使うのが一般的たい。例えばロトカ・ヴォルテラの方程式のような問題にはルンゲ・クッタ法が有効たいね。

例題2 応用問題

~ロトカ・ヴォルテラの方程式~

生態系における被食者と捕食者の生存闘争モデルとして有名なLotka-Volterraモデルは次式の連立常微分方程式である。

ここではxをウサギ(被食者)の個体数、yをキツネ(捕食者)の個体数とする。tは時間である。
ウサギが1000匹、キツネが100匹の状態から、このモデルにより個体数の変化を予測せよ。a = 0.01, b = 0.05, c = 0.0001とする。
出典:ロトカ・ヴォルテラの方程式

ストーク

この問題は、tを独立変数、xyを2つの従属変数とする連立微分方程式たい。
この問題の2つの方程式は以下の様に置くことができる。fxで一つの関数を表していることに注意だ。もちろんfyもや。

ストーク

この連立微分方程式にルンゲ・クッタ法を適用すると、以下の様になる。

シママ

うわぁ~結構複雑になったね・・・。

ストーク

まぁな。でも、独立変数ではh/2やhを、従属変数ではk/2やkを足したうえでfx(t,x,y)やfy(t,x,y)から傾きを求める。というところはさっきの核崩壊の問題と同じたい。
あと、kx1ky1の違いは使っている関数がfxfyかの違いだけたい。

シママ

確かに・・・・。

ストーク

これもExcelで解くことができる。

シママ

計算量は増えたけど、やっていることは核崩壊の問題と同じなのね。
I2セルのkx2を見ると、B6が絶対参照でh,その後のカッコの中でfxを表現している。元々fx = ax-cxyだから、これにkx2のルンゲクッタ法の式を適用すればいいのね。

ストーク

せやな。そうやって一個一個考えていくのが理解への近道やな。方程式が2つになっただけで本質的には
最初に説明した①~④の手順は全く同じ。変数と方程式が一つずつ増えただけたい。

シママ

うん。

ストーク

この方法で解いた、1000日分のウサギとキツネの個体数の予測がこんな感じだ。

シママ

ルンゲ・クッタ法を使えばあらゆる常微分方程式が精度よく解けるようになるのね~。

ストーク

ああ、さっきも言ったが、解く前に解きたい方程式と実現象とのギャップはないか、ケアレスミスでおかしな数値計算結果を出していないか?しっかり吟味して使おうな。

シママ

※追記
ここで解いた連立微分方程式を解くソースコードをVBAで作りました!
こちらもよろしくお願いいたします!

まとめ

ルンゲ・クッタ法はアルゴリズムが簡単な上に高い精度を持つことから、実際の研究や企業におけるR&Dの現場にも利用されており、多くの工学者からの支持を受けています。

今回はコンピュータを使って微分方程式を解くという技術の考え方から始まり、
非常に簡単だが精度は低い微分方程式解法であるオイラー法の説明を行い、
そのオイラー法の発展形として、4次精度ルンゲ・クッタ法について説明しました。

数式とグラフを見比べながら考えると、ルンゲ・クッタ法の考え方がよりわかるようになると思います。

最後に例題を2つ示して、実際にルンゲ・クッタ法を用いて微分方程式を解いて実践による使い方を説明しました。

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です