« Fedora 8 を使って気づいたこと | トップページ | 日経Linux 12月号 »

2007/11/10

ローレンツ方程式

ローレンツ方程式というのは、MITの気象学者ローレンツさんが提示した非線形方程式です。この方程式はカオス的なふるまいをすることで有名です。

なんか難しげな響きがしますが、数値計算プログラムは簡単にできます。

#include <stdio.h>
#include <stdlib.h>

//---parameter---
int    sgm = 10;
int    r   = 28;
double b   = 8.0/3.0;

//---Lorenz Houteishiki---
double f(double x, double y)
{
	return ( sgm*(-x + y) );
}
double g(double x, double y, double z)
{
	return ( r*x - y - x*z );
}
double h(double x, double y, double z)
{
	return ( -b*z + x*y );
}

int main(void)
{
  int i, n;
  double x0, y0, z0, xn, yn, zn;
  double dt, kf[4], kg[4], kh[4];
  FILE *fp, *gp;
  
  x0 = 1.0;  y0 = 1.0;  z0 = 1.0;
  dt = 0.001;  n  = 200000;

  fp = fopen("lorenz.dat", "w");
  gp = popen("gnuplot -persist", "w");
  
  //===Runge-Kutta===
  for(i = 0 ; i <= n ; i++){
    kf[0] = dt * f(x0, y0);
    kg[0] = dt * g(x0, y0, z0);
    kh[0] = dt * h(x0, y0, z0);
    
    kf[1] = dt * f(x0+kf[0]/2.0, y0+kg[0]/2.0);
    kg[1] = dt * g(x0+kf[0]/2.0, y0+kg[0]/2.0, z0+kh[0]/2.0);
    kh[1] = dt * h(x0+kf[0]/2.0, y0+kg[0]/2.0, z0+kh[0]/2.0);
    
    kf[2] = dt * f(x0+kf[1]/2.0, y0+kg[1]/2.0);
    kg[2] = dt * g(x0+kf[1]/2.0, y0+kg[1]/2.0, z0+kh[1]/2.0);
    kh[2] = dt * h(x0+kf[1]/2.0, y0+kg[1]/2.0, z0+kh[1]/2.0);
    
    kf[3] = dt * f(x0+kf[2], y0+kg[2]);
    kg[3] = dt * g(x0+kf[2], y0+kg[2], z0+kh[2]);
    kh[3] = dt * h(x0+kf[2], y0+kg[2], z0+kh[2]);
    
    xn = x0 + (kf[0] + 2.0*(kf[1]+kf[2]) + kf[3])/6.0;
    yn = y0 + (kg[0] + 2.0*(kg[1]+kg[2]) + kg[3])/6.0;
    zn = z0 + (kh[0] + 2.0*(kh[1]+kh[2]) + kh[3])/6.0;
    
    //---OUTPUT---
    fprintf(fp, "%lf\t%lf\t%lf\t\n", xn, yn, zn);
    
    x0 = xn; y0 = yn; z0 = zn;
  }
  
  fclose(fp);
  fprintf(gp, "splot \"lorenz.dat\" w l\n");
  pclose(gp);

  return 0;
}

このプログラムは常微分方程式の近似解を求める方法としてルンゲ=クッタ法を採用しています。アルゴリズムを理解してしまえば簡単にプログラムにすることができます。(結構、精度もいい)

さらに、このプログラムは出力したファイルをGnuplotに渡してグラフを作成するようにしています。例えばこんな風にグラフが出力されます。

Screenshotgnuplot

これをローレンツアトラクタと言います。

|

« Fedora 8 を使って気づいたこと | トップページ | 日経Linux 12月号 »

C」カテゴリの記事

Numerical Computing」カテゴリの記事

コメント

コメントを書く



(ウェブ上には掲載しません)


コメントは記事投稿者が公開するまで表示されません。



トラックバック

この記事のトラックバックURL:
http://app.cocolog-nifty.com/t/trackback/222903/17032369

この記事へのトラックバック一覧です: ローレンツ方程式:

« Fedora 8 を使って気づいたこと | トップページ | 日経Linux 12月号 »