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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

ああ。もちろん。

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

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

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

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

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

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

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

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

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

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

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

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