【C++】微分方程式の勉強1回目【甘口】

こんにちは、タカです。

最近雨続きだったせいか、チエちゃんがちょっと体調悪そうで心配です。明日はお休みなのでゆっくり休んでくれるといいなぁ

今日は微分方程式のお勉強です。甘口なので超簡単な内容です。

微分・・・方程式・・・とアレルギー反応が出るかもしれませんが、実はみんなも使っているし、よく見ています。

例えば経済の価格変動、人口変化、コロナの感染者数変化など様々な変化に使われています。

今回は人口モデルとロジスティックモデルについてプログラムにするまでを書いていこうと思います。

微分方程式ってなに？
人口モデル(マルサスモデル)
ロジスティックモデル
以上！

微分方程式ってなに？

微分方程式は、簡単に説明すると、とある現象の変化をモデル化したものです。

とある現象は、例えば上で書いたように株価変化や感染者数変化です。

モデル化とはこれまでの変化を数理モデルにする、つまり微分方程式を作ることです。変化を式にするのでその式は微分された状態というわけです。なので式を解くには積分が必要です。

普通の方程式は解くことによって数字が得られますが、微分方程式では関数が得られます。

式を解くことで何がいいかというと、今後の予測に使ったりできるということです。シュミレーションです。

それでは実際人口モデルを例にして見ていきましょう

人口モデル(マルサスモデル)

人口モデルは人口の増え方をモデル化したものです。

(単位時間での人口増加量) = (人口増加率) × (人口)

人口 y を時間 t の関数として y = y(t)、人口増加率を k と書くと

$$\displaystyle \frac{dy}{dt}=ry$$

となります。解くと

$$y(t)=y_{0}e^{rt}$$

となります。これをプログラムで計算して使うとこんな感じにできます。

fInitValue が $y_{0}$ で fIntercept が $r$、ループで回している i が $t$ です。

//------------------------------------------------------------------------------
// @brief	人口モデル
//------------------------------------------------------------------------------
void MyMath::PopulationModel( float fInitValue, float fIntercept )
{
	// 初期値
	float fPopulation = fInitValue;
	// ログ出力
	printf( "Population : %f\n", fPopulation );

	// iを時間と仮定して、とりあえず100回回す
	for( int i = 0; i < 100; ++i )
	{
		// 人口モデル計算
		fPopulation = fInitValue * exp( fIntercept * i );

		// ログ出力
		printf( "%d : %f\n", i, fPopulation );
	}
}

今回は計算だけですが、グラフにすると指数関数的なグラフになります。

ロジスティックモデル

人口モデルだと指数関数的に増えていきますが、現実はある一定の数で上限を迎えます。そこで増加の抑制として環境収容力Kを加えたのがロジスティックモデルです。

$$ \frac{dy}{dt} = ry(\frac{K-y}{K})$$

解くと

$$y(t) = \frac{y_{0}Ke^{rt}}{K-y_{0}+y_{0}e^{rt}}$$

あとは同じようにプログラムにします。

//-----------------------------------------------------------------------------
// @brief			ロジスティックモデル
//------------------------------------------------------------------------------
void MyMath::LogisticModel( float fInitValue, float fIntercept, float fCapacity )
{
	// 初期値
	float fPopulation = fInitValue;
	// ログ出力
	printf( "Population : %f\n", fPopulation );

	// iを時間と仮定して、とりあえず100回回す
	for( int i = 0; i < 100; ++i )
	{
		// ロジスティック方程式計算
		fPopulation = (float)(( fInitValue * fCapacity * 
						(float)exp(fIntercept * (float)i )) /
						( fCapacity - fInitValue + fInitValue * 
						(float)exp(fIntercept * (float)i)));
								
		// ログ出力
		printf( "%d : %f\n", i, fPopulation );
	}
}

グラフにするとKに収束するグラフになります