Runge–Kutta法でLorenzアトラクタを描く

#0. 初めまして
初めての投稿なので温かい目で見てください．
3つの常微分方程式をRunge–Kutta法で連立させて解き，Lorenzアトラクタを描いてみました．

#1. Lorenz方程式
E N. LorenzによってNavier–Stokes方程式から導かれた3つの常微分方程式です．太陽熱で暖められる地表と，上層で冷やされる大気の対流を記述するモデルです．

\begin{align}
\frac{dx}{dt}&=-px+py \\
\frac{dy}{dt}&=-xz+rx-y \\
\frac{dz}{dt}&=xy-bz \\
\end{align}

パラメターを$p=10$，$r=28$，$b=8/3$として，3次元相空間の軌道をプロットします．初期座標は$(x, y, z)=(0.0, 1.0, 2.0)$としています．

#2. ソースコード
CでRunge-Kutta法のプログラムを書きました．特に捻りのないシンプルなコードです．

attLorenz.c

#include <stdio.h>
#include <math.h>

const double p = 10.0;
const double r = 28.0;
const double b = 8.0 / 3.0;

double f(double t, double x, double y, double z) {
	return (p * y - p * x);	//dx/dt = -px + py
}

double g(double t, double x, double y, double z) {
	return (r * x - x * z - y);	//dy/dt = -xz + rx - y
}

double h(double t, double x, double y, double z) {
	return (x * y - b * z);	//dz/dt = xy -bz
}

int main(void) {
double t;
double x = 0.0;
double y = 1.0;
double z = 2.0;
double dt = 0.02;
double tmax = 100.0;
double k1, k2, k3, k4;
double l1, l2, l3, l4;
double m1, m2, m3, m4;

	for (t = 0.0; t <= tmax; t += dt) {
		k1 = dt * f(t, x, y, z);
		l1 = dt * g(t ,x, y, z);
		m1 = dt * h(t ,x, y, z);		
		k2 = dt * f(t + dt / 2.0, x + k1 / 2.0, y + l1 / 2.0, z + m1 / 2.0);
		l2 = dt * g(t + dt / 2.0, x + k1 / 2.0, y + l1 / 2.0, z + m1 / 2.0);
		m2 = dt * h(t + dt / 2.0, x + k1 / 2.0, y + l1 / 2.0, z + m1 / 2.0);
		k3 = dt * f(t + dt / 2.0, x + k2 / 2.0, y + l2 / 2.0, z + m2 / 2.0);
		l3 = dt * g(t + dt / 2.0, x + k2 / 2.0, y + l2 / 2.0, z + m2 / 2.0);
		m3 = dt * h(t + dt / 2.0, x + k2 / 2.0, y + l2 / 2.0, z + m2 / 2.0);
		k4 = dt * f(t + dt, x + k3, y + l3, z + m3);
		l4 = dt * g(t + dt, x + k3, y + l3, z + m3);
		m4 = dt * h(t + dt, x + k3, y + l3, z + m3);
		x += (k1 + 2.0 * k2 + 2.0 * k3 + k4) / 6.0;
		y += (l1 + 2.0 * l2 + 2.0 * l3 + l4) / 6.0;
		z += (m1 + 2.0 * m2 + 2.0 * m3 + m4) / 6.0;
	}
	
return(0);
}

#3. Lorenzアトラクタ
gnuplotでプロットしてみました．
蝶のように綺麗な軌道が描けました．3次元の相空間内で軌道同士は交差せず，自己相似な構造（フラクタル構造）を示します．
Lorenzアトラクタは曲線でも曲面でもない，2と3の間の非整数のフラクタル次元をもつフラクタル集合です．

#4. カオス理論
Lorenz方程式は決定論的な常微分方程式ですが，その軌道は非周期的で初期値によって結果が大きく変化するカオス的な性質を示します．この初期値鋭敏性によって，気象の正確な長期予測は非常に困難であると結論されます．これが有名なバタフライ効果です．

#参考文献
●決定論的非周期な流れ（Deterministic Nonperiodic Flow）：Edward N. Lorenz, Journal of Atmospheric Sciences, Vol.20, pp.130-141, 1963.
●ストロガッツ 非線形ダイナミクスとカオス　数学的基礎から物理・生物・化学・工学への応用まで（丸善出版）