資料請求番号：TS33 TS38
ルンゲ・クッタ（runge-kutta）法は、コンピュータを使った微分方程式解法の基本として工学部の基礎教育に登場します。
中でも4次のルンゲ・クッタ法はアルゴリズムが簡単な上に高い精度を持つことから、実際の研究や企業におけるR＆Dの現場にも利用されており、多くの工学者からの支持を受けています。
今回は、コンピュータを使った微分方程式解法（＝数値計算）の基本的な考え方と
非常に簡単なアルゴリズムの微分方程式解法であるオイラー法
そしてその進化系的存在である4次精度ルンゲ・クッタ法について解説します。
今からやろうとしていることは、こんな常微分方程式をコンピュータを使って解くということだ。
微分方程式をコンピュータで解くにあたって、要となる考え方は？
方程式を離散化させること！コンピュータは基本、四則演算しかできないから、微分積分、微分方程式のような高等数学をコンピュータで取り扱いたいなら四則演算に落とし込む作業が必要なのよね。
せやな。連続である微分方程式を、こんな風に離散化させて、四則演算にする。
そしたら、前のyの値であるyiから、次のyi+1を求められるよな。
もちろん、初期条件x0とy0が与えられていること前提やけどな。
hはxの変化分よね。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で割れば平均になる。
特に、k2とk3の影響が強く出るように設定しとるたい。こういうのを「重み付き平均」という。
kの種類が4つで6で割っているから平均と分かりにくいけど・・・言われてみれば、そうね。
と、いうことは・・・4種類のhφを求めて、それを平均したものを足してyi+1を求めているということになるのね。
せや。単にφ=f(xi,yi)としていたオイラー法に比べたら確実にバージョンアップしとるやろ？
うん。
それで4つのkの求め方って言うのが以下の通り。下の図と何度か見比べてしっかりイメージ掴むんや。
①点(xi,yi)でf(xi,yi)を計算して1回目の傾きとする（青色の直線）。これをxi+1まで伸ばすとyはk1だけ変化する。
②xiからh/2だけx方向へ、yiからk1/2だけy方向に進んだ点でf(xi+h/2,yi+k1/2)を計算し2回目の傾きとする（赤色の矢印）。この傾きで点（xi, yi）からx方向にhだけ進むとyはk2だけ変化する。
※色の同じ直線は互いに平行
③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だけ進むとyはk4だけ変化する。
⑤上記4ステップで求めたk1, k2, k3, k4にそれぞれ1,2,2,1の重みを付けて平均した量を求め、点(xi,yi)からx方向にhだけ進んだ時のy方向変化量を求め、yi+1を決定する。
あ、これ、k1だけでyi+1を決定したら、オイラー法と同じになるのね。
お、そこ気づいたか？せやせや。1点だけで次の点を決めるオイラー法とは違って、4つの点を使って、その平均で次の点を決める。しかも、2個目、3個目はxiとxi+1の中間で傾きをとるたい。こういった工夫によって、関数の急激な変化に強くなるたいね。
確かに、微分方程式の解曲線y = f(x)に近い値が予測出来てるわね。さっきのオイラー法とは大違い。
まぁな。本当はもっと点が近いんだが、見やすくするため、少し距離置いとる。それくらい精度ええたい。
じゃあ、ルンゲ・クッタ法を理解したところで、実践たい。この問題をやってみよう。
放射線汚染の原因となる核種の一つにヨウ素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の表計算で十分やと思うたい。
それで、ルンゲ・クッタ法で解く微分方程式の基本形が左で今回解きたい方程式が右や。
yがN、xがtに相当するのね。
せや。そしたら、ルンゲ・クッタ法の式は次の様になる。
そして、こんな感じで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)を求めると以下の様になる。
赤文字の数式がグラフ上の実線、ルンゲ・クッタ法で求めた解を青の点線で表すと以下の様になる
あ～スゴイ！ピッタリ！！
ああ、でもな、今回は解曲線が単調減少という簡単な形だったから数値計算で再現しやすかったという考え方も必要たい。
解がどんな曲線になるのか想像しながら、精度と計算負荷両面から考えて絶妙な刻み幅を設定して計算するのもエンジニアとしての腕前たいね。
それから、今回は変数分離できる簡単な微分方程式を扱ったが、本来、紙面上で解くのが不可能だったり、高度な数学を用いる必要がある場合にこの方法を使うのが一般的たい。例えばロトカ・ヴォルテラの方程式のような問題にはルンゲ・クッタ法が有効たいね。
～ロトカ・ヴォルテラの方程式～
生態系における被食者と捕食者の生存闘争モデルとして有名なLotka-Volterraモデルは次式の連立常微分方程式である。
ここではxをウサギ（被食者）の個体数、yをキツネ（捕食者）の個体数とする。tは時間である。
ウサギが1000匹、キツネが100匹の状態から、このモデルにより個体数の変化を予測せよ。a = 0.01, b = 0.05, c = 0.0001とする。
出典：ロトカ・ヴォルテラの方程式
この問題は、tを独立変数、xとyを2つの従属変数とする連立微分方程式たい。
この問題の2つの方程式は以下の様に置くことができる。fxで一つの関数を表していることに注意だ。もちろんfyもや。
この連立微分方程式にルンゲ・クッタ法を適用すると、以下の様になる。
うわぁ～結構複雑になったね・・・。
まぁな。でも、独立変数ではh/2やhを、従属変数ではk/2やkを足したうえでfx(t,x,y)やfy(t,x,y)から傾きを求める。というところはさっきの核崩壊の問題と同じたい。
あと、kx1とky1の違いは使っている関数がfxかfyかの違いだけたい。
確かに・・・・。
これもExcelで解くことができる。
計算量は増えたけど、やっていることは核崩壊の問題と同じなのね。
I2セルのkx2を見ると、B6が絶対参照でh,その後のカッコの中でfxを表現している。元々fx = ax-cxyだから、これにkx2のルンゲクッタ法の式を適用すればいいのね。
せやな。そうやって一個一個考えていくのが理解への近道やな。方程式が2つになっただけで本質的には
最初に説明した①～④の手順は全く同じ。変数と方程式が一つずつ増えただけたい。
うん。
この方法で解いた、1000日分のウサギとキツネの個体数の予測がこんな感じだ。
ルンゲ・クッタ法を使えばあらゆる常微分方程式が精度よく解けるようになるのね～。
ああ、さっきも言ったが、解く前に解きたい方程式と実現象とのギャップはないか、ケアレスミスでおかしな数値計算結果を出していないか？しっかり吟味して使おうな。
※追記
ここで解いた連立微分方程式を解くソースコードをVBAで作りました！
こちらもよろしくお願いいたします！
ルンゲ・クッタ法はアルゴリズムが簡単な上に高い精度を持つことから、実際の研究や企業におけるR＆Dの現場にも利用されており、多くの工学者からの支持を受けています。
今回はコンピュータを使って微分方程式を解くという技術の考え方から始まり、
非常に簡単だが精度は低い微分方程式解法であるオイラー法の説明を行い、
そのオイラー法の発展形として、4次精度ルンゲ・クッタ法について説明しました。
数式とグラフを見比べながら考えると、ルンゲ・クッタ法の考え方がよりわかるようになると思います。
最後に例題を2つ示して、実際にルンゲ・クッタ法を用いて微分方程式を解いて実践による使い方を説明しました。
スポンサーリンク