方程式と微分方程式（常微分方程式）は同じ関数の別の表示法である。しかし微分すると絶対値の情報が欠けるので，微分方程式には初期値が必要となる。

【例題1】微分方程式で方程式を再現する<rk01.xlsm>

自然現象や物理現象の数式モデルの多くは微分方程式となる。微分方程式が１階の が与えられれば，初期値から出発して数値積分できる。 このような数値積分の手法としてはルンゲ・クッタRunge-Kutta法が標準的である。

また，多階の常微分方程式も１階の正規形常微分方程式の連立式に書き直せるので，初期値がすべて与えられていればRunge-Kutta法により解くことができる。

ここでは正規形の連立常微分方程式を解くための，「微分方程式解法シート」を提供する。シート内にはRunge-Kutta法 およびRunge-Kutta-Fehlberg法のVBAプログラムを組み込んである。しかし，VBAプログラム中の変数とシート上のセルの値を対応させておくことで，式，定数をVBAプログラム中に書く必要がないようにした。微分方程式はシート上のセル内にセル座標で記述する形式なので，いろいろな常微分方程式の解法に使える。

＞ 参考 （同じようにExcelのセルに記述する形式で常微分方程式を解く方法を考案している例）
舞鶴高専　伊藤教授　講義資料　微分方程式｜熊本大学　趙教授　Excel VBAによる数値計算法についての研究｜POLYMATH ODE_Solver｜ 北陸電力　ミックスオイラー法｜Excelで使える数値計算ツールSUITEXL

【例題2】生態系のモデル<ｒｋ02.xlsm>

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

ここではx をウサギ（被食者）の個体数，y をキツネ（捕食者）の個体数とする。t は時間で，両式の第２項(xy )が個体群間の相互作用を表す「交差項」である。ウサギが1000匹，キツネが100匹の状態から，このモデルにより個体数の変化を予測せよ。a =0.01, b =0.05, c =0.0001とする。

図が「微分方程式解法シート」である。セルB1に微分方程式の数2を書く。定数をG2:G4に記しておく。B5に式(1)を，C5に式(2)を記述する。その際，式中の変数x, y はセルB3, C3を用いる。時間t の積分区間と刻み幅をB7:B9に設定する。x, y の初期値1000, 100をB12,C12に記入して，ボタン”Runge-Kutta”をクリックすると連立常微分方程式の積分が実行される。結果がシートの13行以下に出力された。 計算結果を図中のグラフで示す。初期はキツネが増え，ウサギが減少するが，ウサギが減少しすぎるとキツネも減少し，その後ウサギは増加する。このような生態系の周期変化が，簡単なモデル式により表現されていることがわかる。

【例題3】逐次反応 <rk03.xlsm>

反応物AがRを経てSになる簡単な逐次反応：

では，各成分の濃度を表わす連立常微分方程式が次のように書ける。

ここでC_i は各成分の濃度，r は反応速度，k₁, k₂ は反応速度定数である。k₁ = 0.2，k₂ = 0.2，初期濃度をC_A= 1, C_R = 0, C_S = 0 として濃度の時間変化を示せ。

図が「微分方程式解法シート」である。定数をG2:G3に，上式をB5:D5に記述する。積分区間・刻み幅，初期値(B12:D12)を設定し，ボタンクリックで積分を実行する。図中のグラフに各成分の濃度変化の様子を示す。

【例題】Schrödinger 方程式 <schroed.xlsm>

量子化学の基礎式，シュレーディンガ−方程式Schrödinger equationは波動関数wave function φ に関する２階の常微分方程式である。（定常状態。左辺の定数 2m/h²は1とする。）

E が固有値であり，波動関数の解は固有値E のとびとびの値で存在する。ここではV に井戸型ポテンシャルwell-type potential：

を考えて解いてみる。 （参考資料：物理のかぎしっぽ：井戸型ポテンシャルの数値解）基礎式を次式の正規形の連立常微分方程式 ：

として解く。固定刻み幅のRunge-Kutta法による。

【例題4】粒子の飛跡<rk04.xlsm>

重力gのみ考慮して，射出された物体の軌跡を求める。水平方向をx 座標，垂直方向をy 座標とし，各座標の粒子の位置[m]をx, y，粒子の速度[m/s]をu_x, u_yとし，時間をt [s]とすると粒子の運動方程式は次式である。

【参考】粒子の飛跡シミュレーション<fluid3.xlsm>
床面上0.2 mの高さの点から，密度ρ _p = 2500 kg/m³, 粒径D _p = 0.0003 mの球形粒子が初速度u ₀ = 15 m/s, 迎角 45度(=α₀)で静止空気中に放出された。このときの粒子の飛跡を求め，床面に落ちる位置を決定せよ。ただし，空気密度ρ _f = 1.2 kg/m³, 空気粘度μ _f = 0.000018 Pa･sである。

 水平方向をx 座標，垂直方向をy 座標とし，各座標の粒子の位置[m]をx, y，粒子の速度[m/s]をu_x, u_yとし，時間をt [s]とすると粒子の運動方程式は次式である。

また，u_x, u_yの初期値は共にu₀sinα₀=u₀cosα₀=10.61 m/sである。よって，x, y, u_x, u_y に関する連立常微分方程式(11.5)〜(11.8)を，初期条件t =0でx = 0, y =0.2, u_x=u_y= 10.61により解く問題となる。

図が「微分方程式解法シート」である。諸定数およびu, Re, C_Dの計算式をセルG2:G9に記入・記述する。このときu_x, u_yはD3, E3とする。方程式数(B1)を4とし，B5:E5に微分方程式(11.5)〜(11.8)を記述する。積分区間と刻み幅を設定し，初期値をB12:E12に書く。ボタンクリックでRunge-Kutta法による積分が実行される。その結果から得られた粒子の飛跡(x-y)を図中に示す。また床面に落ちる位置(y =0 )は計算結果よりx =0.1735 mである。

【例題】主人と犬問題 the master-and-dog problem<dog.xls>

主人はy 軸上を(0,0)から出発して速度v で歩く。犬は主人の方向に向かって速度w で走る。犬の位置を(x ,y )  として，犬の軌跡は次式の２階の微分方程式となる。

 ここでy ’ =z とすると，２階の微分方程式が以下の連立常微分方程式となる。 刻み幅可変のFehlberg法により解く。

単振動の微分方程式は重力とフックの法則より、

 である。x は位置，おもりの質量をM, kをバネ定数として，ω² =k/M である。

この２階の常微分方程式をルンゲ・クッタ法により解く。正規形とするため速度v (=dx/dt)を用いると，上式は

のような連立常微分方程式となる。初期条件は、
x =1, v =0 （距離１引っ張って離す場合） または
x =0, v =1 （静止状態からおもりに速度を与える場合）
などが与えられる。

 【例題5】単振動 <rk05.xlsm>

表計算+BASICで振動の微分方程式の解を求める。これはバネ振動のシミュレーションと言える。初期条件を変えて振動の様子をみる。

さらに、図のようにダンパーを加えた系を考える。

力のバランスは、
（質量Mによる慣性力）＋（ダンパーの制動力）＋（ばねの力）＝（外力） より、

である。これは速度v を位置ｘ で書き換えると、v =dx/dt なので，x に関する２階の常微分方程式：

 である。 (1)式は以下の正規形連立常微分方程式となる。

【例題6】バネマスダンパ系<rk06.xlsm>

図のようなR(抵抗), L(コイル), C(コンデンサ)からなる回路で，スイッチを入れたときの過渡現象を調べる。直列回路の関係式より，

である。これは電荷 q と電流 i に関する次式の連立常微分方程式となる。

【例題7】R-L-C直列回路の過渡現象<rk07.xlsm>