My Blog Logo
Blog
数値シミュレーション －微分方程式－

数値シミュレーション －微分方程式－

数値シミュレーション解説 第2回

カテゴリ:

シミュレーション

タグ:

シミュレーション数学

Published: 7/17/2018

Updated: 12/19/2025

約15分

数値シミュレーション解説の前回の記事

数値シミュレーション －数理モデル－

数値シミュレーション解説の次回の記事

数値シミュレーション －オイラー法－

更新履歴
  • 2025/12/19: プログラムを C# 14.0 に修正
こんにちは、Rafkaです。
今回は、多くの数理モデルに使用される微分方程式についてのさらなる解説と、本シリーズでの共通の表記法を定義します。

微分方程式

微分方程式 は、未知の関数とその導関数との関係を表した式で、 微分方程式を解くことは、未知関数の具体的な形を決めることを意味します。 解く際には積分が必要不可欠となるので、解には積分定数が含まれ、 積分定数が含まれた状態で表される解を 一般解 、 初期値等の情報を与えることで積分定数の値まで定めることで得られる解を 特殊解 と言います。 コンピュータで数値的に解く方法は、式変形を行って解く方法ではないので、 コンピュータによって求まる解は特殊解となります。
微分方程式には大きく分けて、 未知関数が本質的に 1 変数関数である 常微分方程式 （ Ordinary Differential Equation; O.D.E. ）と、 未知関数が多変数関数であり、 その偏導関数を用いて表される 偏微分方程式 （ Partial Differential Equation; P.D.E. ） の 2 種類があります。
常微分方程式 - Wikipedia のプレビュー画像
常微分方程式 - Wikipedia
ja.wikipedia.org のファビコン
ja.wikipedia.org
偏微分方程式 - Wikipedia のプレビュー画像
偏微分方程式 - Wikipedia
ja.wikipedia.org のファビコン
ja.wikipedia.org
それぞれについて解説していきます。

常微分方程式（O.D.E.）

nn 階の常微分方程式は 1 階の連立常微分方程式に変形することができますので、 常微分方程式は一般に以下のような形式で表すことができると言えます。
dx(t)dt=f(x(t),t)\frac{d \boldsymbol{x}(t)}{d t} = \boldsymbol{f} \bigl( \boldsymbol{x}(t), t \bigr)
ここで tt は変数、x(t)\boldsymbol{x}(t) は未知関数、 f(x,t)\boldsymbol{f}(\boldsymbol{x}, t) は右辺の式をそれぞれ表していて、 未知関数とf \boldsymbol{f} はそれぞれ以下のようなベクトルになります。
x(t)=[x1(t),x2(t),⋯ ,xn(t)]T\boldsymbol{x}(t) = [x_1(t), x_2(t), \cdots, x_n(t)]^T f(x(t),t)=[f1(x(t),t),f2(x(t),t),⋯ ,fn(x(t),t)]T\boldsymbol{f}\bigl(\boldsymbol{x}(t), t\bigr) = \Bigl[ f_1\bigl(\boldsymbol{x}(t), t\bigr), f_2\bigl(\boldsymbol{x}(t), t\bigr), \cdots, f_n\bigl(\boldsymbol{x}(t), t\bigr) \Bigr]^T
常微分方程式の具体例として最も有名なのはニュートンの運動方程式でしょう。
ニュートンの運動方程式は、 位置を時間に対する未知関数 x(t)x(t) とした時の以下のような 2 階の常微分方程式であると言えます。
md2x(t)dt2=Fm \frac{d^2 x(t)}{d t^2} = F
また、位置 x(t)x(t) 、速度 v(t)v(t) 、加速度 a(t)a(t) には以下のような関係があるので、
dx(t)dt=v(t),    dv(t)dt=a(t)\frac{d x(t)}{d t} = v(t), \;\; \frac{d v(t)}{d t} = a(t)
運動方程式の例を先程の式の表記に置き換えるとすると、 x(t)\boldsymbol{x}(t) と f(x(t),t)\boldsymbol{f}\bigl(\boldsymbol{x}(t), t\bigr) はそれぞれ以下のようになる。
x(t)=[x(t),v(t)]T\boldsymbol{x}(t) = [x(t), v(t)]^T f(x(t),t)=[v(t),F/m]T\boldsymbol{f}\bigl(\boldsymbol{x}(t), t\bigr) = [v (t), F / m]^T
この様にして、 2 階の常微分方程式であるニュートンの運動方程式を、 1 階の連立微分方程式に変換することができたので、 同様にして、高次の導関数を 1 次の導関数の連結で表現することで、 高階の常微分方程式を 1 階の連立常微分方程式に変換できることが分かりますね。

偏微分方程式（P.D.E.）

常微分方程式が、1 変数関数の導関数を用いて構成される方程式であったのに対して、 偏微分方程式は多変数関数の偏導関数を用いて構成される方程式になります。
例としては、前回の記事で示した以下のナビエ－ストークス方程式や、シュレディンガーの波動方程式があります。
∂v∂t+(v⋅∇)v=−1ρ∇p+ν∇2v+g\frac{\partial \boldsymbol{v}}{\partial t} + (\boldsymbol{v} \cdot \nabla) \boldsymbol{v} = - \frac{1}{\rho} \nabla p + \nu \nabla^2 \boldsymbol{v} + \boldsymbol{g}
本シリーズの目標である脳シミュレーションに使用する数理モデルは常微分方程式で表されるので、 偏微分方程式は紹介程度に留めさせてもらいます。

プログラム

数値計算を行うプログラムも、このシリーズでは示していこうと思います。
使用する言語は C# 14 になります。
数理モデルや数値計算の手法を簡単に入れ替えて試せるように、以下のような構成にします。
  • ModelFunction Delegate: 数理モデルの関数を表すデリゲート
  • SolverFunction Delegate: 数値計算手法を表すデリゲート
  • Simulate Method: 数理モデルと数値計算手法を受け取り、シミュレーションを実行するメソッド

Function Delegates

ModelFunction Delegate

ModelFunction デリゲートは、数理モデルの関数を表すデリゲートです。 常微分方程式は一般に、以下の数式で表せることは既に説明した通りですので、
dx(t)dt=f(x(t),t)\frac{d \boldsymbol{x}(t)}{d t} = \boldsymbol{f} \bigl( \boldsymbol{x}(t), t \bigr)
ModelFunction デリゲートはこの右辺の関数 f\boldsymbol{f} を表すものとして、以下のように定義しました。
delegate void ModelFunction(ReadOnlySpan<double> x, double t, Span<double> dx);
計算に必要なメモリは呼び出し元で確保して渡す形にして、 ReadOnlySpan<double> x と double t が関数の入力、Span<double> dx が関数の出力を表しています。

SolverFunction Delegate

SolverFunction デリゲートは、数値計算手法を表すデリゲートです。 数値計算手法は、数理モデルの関数を使用して数理モデルを解くものなので、 以下のように定義しました。
delegate void SolverFunction(ModelFunction model, Span<double> x, double t, double dt);
model で渡された数理モデルの関数に対して、入力 Span<double> x と double t を使用して計算を行い、 刻み幅 double dt を使用して数理モデルを解き、結果を Span<double> x に上書きします。
例えば、オイラー法を実装すると以下のようになります。
void Euler(ModelFunction f, Span<double> x, double t, double dt) {
  int n = x.Length;
  Span<double> dx = stackalloc double[n];
 
  f(x, t, dx);
  ScaleAdd(x, dx, dt, x);
}
ここで、ScaleAdd メソッドは、ベクトルのスケーリングと加算を行うユーティリティメソッドです。
void ScaleAdd(ReadOnlySpan<double> a, ReadOnlySpan<double> b, double s, Span<double> result) {
  for (int i = 0; i < result.Length; i++) {
    result[i] = a[i] + b[i] * s;
  }
}

Simulate Method

指定された数理モデルと数値計算手法を使用してシミュレーションを実行し、結果を数値の配列のリストとして出力します。 以下のように定義しました。
List<double[]> Simulate(
  ModelFunction model,
  Span<double> state,
  SolverFunction solver,
  double dt,
  int steps
) {
  var n = state.Length;
  var list = new List<double[]>();
 
  list.Add([0, ..state]);
  for (var i = 0; i < steps; i++) {
    var t = i * dt;
    solver(model, state, t, dt);
    list.Add([t + dt, ..state]);
  }
 
  return list;
}
state は初期状態を表すベクトルで、dt は刻み幅、steps はシミュレーションのステップ数を表しています。 state にはシミュレーションの途中経過が上書きされていきます。 各ステップの結果は、15 行目の処理でリストに新しい配列を都度生成して保存しています。

使用例

今まで示してきたプログラムを使って、実際に簡単な例について解いてみましょう。
例題に使うのは以下のような連立微分方程式です。
{dx1(t)dt=x2(t)dx2(t)dt=−x1(t)\left\{ \begin{align*} \frac{d x_1(t)}{d t} &= &x_2(t) \\ \frac{d x_2(t)}{d t} &= &-x_1(t) \end{align*} \right.
これは単振動の数理モデルで、方程式の解は、x1(t)=sin⁡(t)x_1(t) = \sin(t) 、x2(t)=cos⁡(t)x_2(t) = \cos(t) となります。それは、
dx1(t)dt=dsin⁡(t)dt=cos⁡(t)=x2(t)\begin{split} \frac{d x_1(t)}{d t} &= \frac{d \sin(t)}{d t} \\ &= \cos(t) \\ &= x_2(t) \end{split} dx2(t)dt=dcos⁡(t)dt=−sin⁡(t)=−x1(t)\begin{split} \frac{d x_2(t)}{d t} &= \frac{d \cos(t)}{d t} \\ &= -\sin(t) \\ &= -x_1(t) \end{split}
であることから確認できます。
さて、次はこの問題をコンピュータで解いても同じ結果が得られるかどうかを確認してみましょう。
解くのに使用するのは、今まで示してきたプログラムになります。
私は今回、以下のようなプログラムを作成してみました。
#!/usr/bin/dotnet run
 
var result = Simulate(
  (x, _, dx) => { dx[0] = x[1]; dx[1] = -x[0]; },
  stackalloc double[] { 0.0, 1.0 },
  Euler,
  0.1,
  100
);
 
foreach (var item in result) {
  Console.WriteLine($"Time={item[0]:F2} Pos={item[1]:F4} Vel={item[2]:F4} Exact={Math.Sin(item[0]):F4}");
}
 
void ScaleAdd(ReadOnlySpan<double> a, ReadOnlySpan<double> b, double s, Span<double> result) {
  for (int i = 0; i < result.Length; i++) {
    result[i] = a[i] + b[i] * s;
  }
}
 
void Euler(ModelFunction f, Span<double> x, double t, double dt) {
  int n = x.Length;
  Span<double> dx = stackalloc double[n];
 
  f(x, t, dx);
  ScaleAdd(x, dx, dt, x);
}
 
List<double[]> Simulate(ModelFunction model, Span<double> state, SolverFunction solver, double dt, int steps) {
  var n = state.Length;
  var list = new List<double[]>();
 
  list.Add([0, ..state]);
  for (var i = 0; i < steps; i++) {
    var t = i * dt;
    solver(model, state, t, dt);
    list.Add([t + dt, ..state]);
  }
 
  return list;
}
 
delegate void ModelFunction(ReadOnlySpan<double> x, double t, Span<double> dx);
delegate void SolverFunction(ModelFunction model, Span<double> x, double t, double dt);
簡単にプログラムを解説すると、 4 行目の (x, _, dx) => { dx[0] = x[1]; dx[1] = -x[0]; } で今回の数理モデルの関数を定義していて、 5 行目の stackalloc double[] { 0.0, 1.0 } で初期状態を定義しています。 6 行目の Euler で数値計算手法としてオイラー法を指定し、 7 行目の 0.1 で刻み幅、8 行目の 100 でステップ数を指定しています。 得られた結果を 12 行目でコンソールに出力していますが、 その際に厳密解である sin⁡(t)\sin(t) も一緒に出力しています。
C# 14 を動かす環境があれば、.NET 10 を使用して上記のプログラムのファイルをそのまま実行できると思います。 上記のプログラムを simulation-diff-eq.cs というファイルで保存した場合、以下のコマンドで実行可能です。
dotnet simulation-diff-eq.cs
得られた計算結果をプロットしてみると、以下のようなグラフが得られました。
赤線が x1(t)x_1(t) 、青線が x2(t)x_2(t) を表しています。
グラフからは確かに、x1(t)=sin⁡(t)x_1(t) = \sin(t) で、x2(t)=cos⁡(t)x_2(t) = \cos(t) である様子が確認できます。
これで、コンピュータを使用して常微分方程式を数値的に解けるということが分かっていただけたと思います。

まとめ

今回は、数理モデルを記述するための微分方程式について説明し、最後には実際に数値計算によって解く様子をお見せしました。
次回は、今回微分方程式を解くために使用したオイラー法の解説を行います。
Rafka
Rafka

Software Engineer

最近はドラクエ 1 & 2 をプレイしてます。楽しい！

数値シミュレーション解説の前回の記事

数値シミュレーション －数理モデル－

数値シミュレーション解説の次回の記事

数値シミュレーション －オイラー法－

前回の記事

数値シミュレーション －数理モデル－

次回の記事

発火率計算機を作成しました

関連記事

数値シミュレーション －ローレンツ・アトラクタ
数値シミュレーション －ローレンツ・アトラクタ

シミュレーション

約4分

数値シミュレーション解説 第8回

シミュレーション
数学
カオス

カオスの代表例であるローレンツ・アトラクタについて解説し、3D ビジュアライゼーションで初期値鋭敏性を体験します。

数値シミュレーション －二重振り子－
数値シミュレーション －二重振り子－

シミュレーション

約5分

数値シミュレーション解説 第7回

シミュレーション
数学

ルンゲ＝クッタ法を用いて二重振り子のシミュレーションを行います。

数値シミュレーション －各解法の誤差比較－
数値シミュレーション －各解法の誤差比較－

シミュレーション

約10分

数値シミュレーション解説 第6回

シミュレーション
数学

今まで紹介してきた数値解法の誤差比較を、斜方投射を題材に行っていきます。

数値シミュレーション －ルンゲ＝クッタ法－
数値シミュレーション －ルンゲ＝クッタ法－

シミュレーション

約6分

数値シミュレーション解説 第5回

シミュレーション
数学

常微分方程式の数値解法の 1 つであるルンゲ＝クッタ法の導出と実装を行います。

数値シミュレーション －ホイン法－
数値シミュレーション －ホイン法－

シミュレーション

約11分

数値シミュレーション解説 第4回

シミュレーション
数学

常微分方程式の数値解法の 1 つであるホイン法の導出と実装を行います。