久しぶりに Runge-Kutta法でローレンツ方程式を解いたよ
Runge-Kutta法っていうのは、常微分方程式を数値計算するための、一解法です。
ルンゲ＝クッタ法 - Wikipedia
ローレンツ方程式というのは、カオス理論にとって
極めて重要な方程式であります。
ローレンツ方程式 - Wikipedia
で、これを解いて何をしようとしているか、というと
プロフィール画像をつくろう！というわけです。
kei-os2007という idは、学生時代の
chaos理論との素晴らしい出会いから、きています。
自分の人生は、どんなアトラクタにひきこまれていくのだろう、
というワクワクする気持ちも入っています。
というわけで、解きますー。
Cソース
Cで書いてみます。
Wikipedia には y' = f(t,y) ってなってますが
ここでは、より一般的（と思われる） f(x,y) を使って書いてます。
きざみ幅の h は step にしています。
lorenz.c ということで。
$ cat lorenz.c #include <stdio.h> #include <stdlib.h> #include <string.h> #define OLD (0) #define NEW (1) // 3D data type typedef struct { double x; double y; double z; } Point; // parameter double Pr = 10.0; double b = 8.0/3.0; double r = 28.0; static double dxdt(double x, double y, double z) { return Pr*(y - x); } static double dydt(double x, double y, double z) { return (-x*z + r*x - y); } static double dzdt(double x, double y, double z) { return (x*y - b*z); } static double rk4_3d(double (*proc)(double x, double y, double z), Point *p, double step) { double k0, k1, k2, k3; k0 = step*proc(p->x, p->y, p->z); k1 = step*proc(p->x*k0/2, p->y+k0/2, p->z+k0/2); k2 = step*proc(p->x+k1/2, p->y+k1/2, p->z+k1/2); k3 = step*proc(p->x+k2, p->y+k2, p->z+k2); return ((k0 + 2*k1 + 2*k2 + k3)/6); } int main(int argc, char **argv) { double step = 0.001; int i, iter = 100000; Point points[2]; FILE *fp; if (NULL == (fp = fopen("./lorenz.dat", "w"))) { fprintf(stderr, "cannot open output file\n"); exit(-1); } // initial condition points[OLD].x = 0.1; points[OLD].y = 0.1; points[OLD].z = 0.1; for (i = 0; i < iter; i++) { points[NEW].x = points[OLD].x + rk4_3d(dxdt, &points[OLD], step); points[NEW].y = points[OLD].y + rk4_3d(dydt, &points[OLD], step); points[NEW].z = points[OLD].z + rk4_3d(dzdt, &points[OLD], step); fprintf(fp, "%lf\t%lf\t%lf\n", points[NEW].x, points[NEW].y, points[NEW].z); // for next calc. memcpy((void *)&points[OLD], (void *)&points[NEW], sizeof(Point)); } fclose(fp); return 0; }
