Numerical Computing

2008/12/23

BiCGSTAB法

双共役勾配(BiCG)法の残差を減少させ安定化させたのがBiCGSTAB法。BiCGのCGとの違いは非対称行列の反復計算に使用できること。(詳しく解説しようとすると、とても難しいので割愛)

とりあえず、サンプルプログラムを作ってみた。(用いた方程式は以前CG法のプログラムで使ったものと同様)

/*==================================================*/
// Bi-CGSTAB法テストプログラム
/*==================================================*/

#include <cstdio>
#include <cstdlib>
#include <cmath>

using namespace std;

const int    N      = 10;
const int    ITRMAX = 100;
const double EPS    = 1.0e-6;

typedef double Vector[N];      // ベクトル
typedef double Matrix[N][N];   // 行列

// ベクトルに行列を作用 y = Ax
void matrix_x_vector(Vector y, const Matrix a, const Vector x)
{
  for(int i = 0; i < N; i++) {
    double mxv = 0.0;
    for(int j = 0; j < N; j++) {
      mxv += a[i][j]*x[j];
    }
    y[i] = mxv;
  }
}

// 内積を計算
double inner_product(const Vector x, const Vector y)
{
  double inner_p = 0.0;
  for(int i = 0; i < N; i++) {
    inner_p += x[i]*y[i];
  }
  return inner_p;
}

// ベクトルノルムを計算
// ベクトルノルム := sgm(0〜N-1)|x[i]|
double vector_norm(const Vector x)
{
  double norm = 0;
  for(int i = 0; i < N; i++){
    norm += fabs(x[i]);
  }
  return norm;
}

// BiCGSTAB法
/*
  [Algorithm] Ax = b を BiCGSTAB法で解く
  1) Compute r0 = b - Ax (Choose r0* arbitary, r0* = r0, (r0, r0*) != 0)
             p = r = r0
     (c1 = (r0*, r0) = (r0, r0)とする)
  2) Loop...
     c2    = (r0*, Ap)
     alpha = c1 / c2
     e     = r - alpha*Ap
     c3    = (e, Ae) / (Ae, Ae)

     x_new = x + alpha*p + c3*e
     r_new = e - c3*Ae
     if(EPS > ||r||/||b||) break

     c1    = (r_new, r0*) = (r_new, r0)
     beta  = c1 / (c2*c3)
     p_new = r_new + beta*(p - c3*Ap)
*/
void bicgstab(const Matrix a, Vector x, const Vector b)
{
  static Vector r0;  // 初期残差ベクトル
  static Vector r;   // 残差ベクトル
  static Vector p;   // 探索方向ベクトル

  static Vector ax;
  matrix_x_vector(ax, a, x);
  // r0,r,pを計算 r0 := b - Ax
  // p = r = r0
  for(int i = 0; i < N; i++) {
    r0[i] = b[i] - ax[i];
    p[i]  = r[i] = r0[i];
  }
  // (r0, r0*) != 0,   r0* = r0
  double c1 = inner_product(r0, r0);
  if(c1 == 0.0) {
    fprintf(stderr, "(r0, r0*) == 0!!\n");
    exit(1);
  }

  for(int iter = 1; iter < ITRMAX; iter++) {
    static Vector ap;
    matrix_x_vector(ap, a, p);
    double c2    = inner_product(r0, ap);
    double alpha = c1 / c2;

    static Vector e;
    for(int i = 0; i < N; i++) {
      e[i] = r[i] - alpha*ap[i];
    }
    static Vector ae;
    matrix_x_vector(ae, a, e);
    double e_dot_ae  = inner_product(e, ae);
    double ae_dot_ae = inner_product(ae, ae);
    double c3        = e_dot_ae / ae_dot_ae;

    for(int i = 0; i < N; i++) {
      x[i] += alpha*p[i] + c3*e[i];
      r[i]  = e[i] - c3*ae[i];
    }
    
    double err = vector_norm(r)/vector_norm(b);
    printf("LOOP : %d\t Error : %g\n", iter, err);
    if(EPS > err) break;

    c1          = inner_product(r, r0);
    double beta = c1 / (c2*c3);
    
    for(int i = 0; i < N; i++) {
      p[i] = r[i] + beta*(p[i] - c3*ap[i]);
    } 
  }
}

int main()
{
  // 連立方程式 Ax = b
  Matrix a = { {5.0, 2.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	       {2.0, 5.0, 2.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	       {0.0, 2.0, 5.0, 2.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	       {0.0, 0.0, 2.0, 5.0, 2.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	       {0.0, 0.0, 0.0, 2.0, 5.0, 2.0, 0.0, 0.0, 0.0, 0.0},
	       {0.0, 0.0, 0.0, 0.0, 2.0, 5.0, 2.0, 0.0, 0.0, 0.0},
	       {0.0, 0.0, 0.0, 0.0, 0.0, 2.0, 5.0, 2.0, 0.0, 0.0},
	       {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 2.0, 5.0, 2.0, 0.0},
	       {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 2.0, 5.0, 2.0},
	       {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 2.0, 5.0} };
  Vector b = {3.0, 1.0, 4.0, 0.0, 5.0, -1.0, 6.0, -2.0, 7.0, -15.0};
  // 初期値は適当
  Vector x = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};

  // CG法でAx=bを解く
  bicgstab(a, x, b);
  
  printf("###Calc End.###\n");
  for(int i = 0; i < N; i++) {
    printf("x[%d] = %2g\n", i, x[i]);
  }
  
  return 0;
}

実行結果は次のようになる。

LOOP : 1         Error : 0.416594
LOOP : 2         Error : 0.126537
LOOP : 3         Error : 0.0364201
LOOP : 4         Error : 0.0123132
LOOP : 5         Error : 0.00372894
LOOP : 6         Error : 0.0011802
LOOP : 7         Error : 0.00023157
LOOP : 8         Error : 4.34888e-05
LOOP : 9         Error : 4.5223e-06
LOOP : 10        Error : 7.16945e-20
###Calc End.###
x[0] =  1
x[1] = -1
x[2] =  2
x[3] = -2
x[4] =  3
x[5] = -3
x[6] =  4
x[7] = -4
x[8] =  5
x[9] = -5

解は得られているのでBiCGSTABのコードとしては多分あってるはず。

| | コメント (0) | トラックバック (0)

2008/12/09

CG法

まずはCG法(Conjugate Gradient Method)で連立方程式を解く練習。

CG方のアルゴリズムはWikipediaに任せることにして、CG法のプログラムは以下のようになった。

/*==================================================*/
// CG法テストプログラム
/*==================================================*/

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

#define N    10
#define TMAX 100
#define EPS  (1.0e-6)

typedef double Vector[N];       // ベクトル
typedef double Matrix[N][N];    // 行列

// ベクトル初期化
void init_vector(Vector x)
{
  int i;
  for(i = 0; i < N; i++){
    x[i] = 0;
  }
}

// ベクトルに行列を作用 y = Ax
void vector_x_matrix(Vector y, Matrix a, Vector x)
{
  int    i, j;
  double vxm;
  for(i = 0; i < N; i++){
    vxm = 0;
    for(j = 0; j < N; j++){
      vxm += a[i][j]*x[j];
    }
    y[i] = vxm;
  }
}

// 内積を計算
double dot_product(Vector x, Vector y)
{
  int    i;
  double dot_p = 0;
  for(i = 0; i < N; i++){
    dot_p += x[i]*y[i];
  }
  return dot_p;
}

// ベクトルノルムを計算
// ベクトルノルム := sgm(0〜N-1)|x[i]|
double vector_norm(Vector x)
{
  int    i;
  double norm = 0;
  for(i = 0; i < N; i++){
    norm += fabs(x[i]);
  }
  return norm;
}

// CG法
void cg_method(Matrix a, Vector x, Vector b)
{
  static Vector p;   // 探索方向ベクトル
  static Vector r;   // 残差ベクトル
  static Vector ax;
  static Vector ap;
  int           i, iter;

  // Axを計算
  vector_x_matrix(ax, a, x);

  // pとrを計算 p = r := b - Ax
  for(i = 0; i < N; i++){
    p[i] = b[i] - ax[i];
    r[i] = p[i];
  }

  // 反復計算
  for(iter = 1; iter < TMAX; iter++){
    double alpha, beta, err = 0;

    // alphaを計算
    vector_x_matrix(ap, a, p);
    alpha = +dot_product(p, r)/dot_product(p, ap);
   
    for(i = 0; i < N; i++){
      x[i] += +alpha*p[i];
      r[i] += -alpha*ap[i];
    }
   
    err = vector_norm(r);   // 誤差を計算
    printf("LOOP : %d\t Error : %g\n", iter, err);
    if(EPS > err) break;

    // EPS < err ならbetaとpを計算してループを継続
    beta = -dot_product(r, ap)/dot_product(p, ap);
    for(i = 0; i < N; i++){
      p[i] = r[i] + beta*p[i];
    }
  }
}

int main(void)
{
  // 連立方程式 Ax = b
  // 行列Aは正定値対象行列
  Matrix a = { {5.0, 2.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
       {2.0, 5.0, 2.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
       {0.0, 2.0, 5.0, 2.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
       {0.0, 0.0, 2.0, 5.0, 2.0, 0.0, 0.0, 0.0, 0.0, 0.0},
       {0.0, 0.0, 0.0, 2.0, 5.0, 2.0, 0.0, 0.0, 0.0, 0.0},
       {0.0, 0.0, 0.0, 0.0, 2.0, 5.0, 2.0, 0.0, 0.0, 0.0},
       {0.0, 0.0, 0.0, 0.0, 0.0, 2.0, 5.0, 2.0, 0.0, 0.0},
       {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 2.0, 5.0, 2.0, 0.0},
       {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 2.0, 5.0, 2.0},
       {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 2.0, 5.0} };
  Vector b = {3.0, 1.0, 4.0, 0.0, 5.0, -1.0, 6.0, -2.0, 7.0, -15.0};
  // 初期値は適当
  Vector x = {1.0, 1.0, 1.0, 1.0, 1.0,  1.0, 1.0,  1.0, 1.0,   1.0};
  int    i;

  // CG法でAx=bを解く
  cg_method(a, x, b);
 
  printf("###Calc End.###\n");
  for(i = 0; i < N; i++){
    printf("x[%d] = %2g\n", i, x[i]);
  }
 
  return 0;
}

これを実効すると以下のようになる。

$ ./a.out
LOOP : 1         Error : 33.625
LOOP : 2         Error : 19.5236
LOOP : 3         Error : 8.64938
LOOP : 4         Error : 3.7894
LOOP : 5         Error : 1.49629
LOOP : 6         Error : 0.757903
LOOP : 7         Error : 0.348423
LOOP : 8         Error : 0.141565
LOOP : 9         Error : 0.0454147
LOOP : 10        Error : 1.17553e-16
###Calc End.###
x[0] =  1
x[1] = -1
x[2] =  2
x[3] = -2
x[4] =  3
x[5] = -3
x[6] =  4
x[7] = -4
x[8] =  5
x[9] = -5

どうやらちゃんと解けているようである。

| | コメント (0) | トラックバック (0)

2008/10/05

2次元キャビティ流れ

2次元の正方形キャビティ内の流れを「流れ関数-渦度法」で変形した支配方程式を差分を使って解く数値計算プログラムを作っている。まだ未完成でいろいろ修正点もあるけど一応それらしい解が得られている。

Psi2_2 

この図は流れ関数の等高線をプロットしたもの。

大雑把に言えば、このプログラムは移流拡散方程式(Navieor-Stoks方程式のrotをとったもの)とポアソン方程式(渦度と流れ関数の関係式)を解いているだけだけど、(キャビティ流れの)理論的な部分を完全に理解していないので他人に説明できるようなレベルじゃない。(この状態だと、"このようにすると何故か動いてくれる"といった説明しかできない)

あと未完成だけどソースコード:「cavity.c」

/*
  [2008]
   2次元 Cavity Flow Simulation
*/

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

/* Cavity Parameters */
#define L  1.0    // 1辺の長さ
#define N  100    // メッシュの数
#define DX (L/N)  // メッシュ間隔(x)
#define DY (L/N)  // メッシュ間隔(y)

/* Fluid Parameters */
#define U  1.0       // 流速
#define RE 100       // レイノルズ数
#define NU (U*L/RE)  // 動粘性率

/* CFL Condition */
#define CFL     0.3                                    // 安全係数alpha
#define DT_ADV  (1.0/(U/DX))                           // 移流時間
#define DT_DIFF (1.0/(2*NU*(1.0/DX/DX + 1.0/DY/DY)))   // 拡散時間
#define DT      (CFL*(1.0/(1.0/DT_ADV + 1.0/DT_DIFF))) // 時間刻み

/* Calc Step */
#define TMAX 10000   // 計算回数

/* Poisson : Gauss-Sidel Parameters */
#define LOOP_MAX (N*N/4)   // ガウスザイデルの計算回数上限
#define EPS      (1.0e-4)  // 精度


/* Define New Data Type */
typedef double Array[N+1][N+1];


/* Initialize Array Type */
void initArray(Array var)
{
  int i, j;
  for(i = 0; i <= N; i++){
    for(j = 0; j <= N; j++)
      var[i][j] = 0.0;
  }
}


/* Omega's Boundary Conditions */
void setBoundary_omega(Array omega, Array psi)
{
  int ix, iy;
  for(ix = 0; ix <= N; ix++){
    omega[ix][0] = -2.0*(psi[ix][0+1]     )/DY/DY;   // bottom
    omega[ix][N] = -2.0*(psi[ix][N-1]+U*DY)/DY/DY;   // top
  }
  for(iy = 0; iy <= N; iy++){
    omega[0][iy] = -2.0*psi[0+1][iy]/DX/DX;          // left
    omega[N][iy] = -2.0*psi[N-1][iy]/DX/DX;          // right
  }
}

/* Psi's Boundary Conditions */
void setBoundary_psi(Array psi)
{
  int ix, iy;
  for(ix = 0; ix <= N; ix++){
    psi[ix][0] = 0.0;   // bottom
    psi[ix][N] = 0.0;   // top
  }
  for(iy = 0; iy <= N; iy++){
    psi[0][iy] = 0.0;   // left
    psi[N][iy] = 0.0;   // right
  }
}
 

/* Operator */
double advection(Array omega, Array psi, int x, int y)
{
  return ( (psi[x+1][y]-psi[x-1][y])*(omega[x][y+1]-omega[x][y-1])
	   -(psi[x][y+1]-psi[x][y-1])*(omega[x+1][y]-omega[x-1][y])
	   )/4.0/DX/DY;
}

double laplacian(Array omega, int x, int y)
{
  return ( (omega[x+1][y]-2.0*omega[x][y]+omega[x-1][y])/DX/DX
	   +(omega[x][y+1]-2.0*omega[x][y]+omega[x][y-1])/DY/DY );
}

/* Solve Omega */
void solveOmega(Array omega, Array psi)
{
  int ix, iy;
  static Array omega_new;

  for(ix = 1; ix < N; ix++){
    for(iy = 1; iy < N; iy++){
      omega_new[ix][iy] = omega[ix][iy]
	+DT*(advection(omega, psi, ix, iy) + NU*laplacian(omega, ix, iy));
    }
  }

  for(ix = 1; ix < N; ix++){
    for(iy = 1; iy < N; iy++){
      omega[ix][iy] = omega_new[ix][iy];
    }
  }
  
}


/* Solve Psi using Gauss-Sidel Method */
void solvePsi(Array psi, Array omega)
{
  int    ix, iy, iter;
  double psiMax = 1.0e-7;
  
  for(iter = 1; iter <= LOOP_MAX; iter++){
    double errMax = 0.0;
    for(ix = 1; ix < N; ix++){
      for(iy = 1; iy < N; iy++){
	double err     = 0.0;
	double psi_old = psi[ix][iy];
	double psi_new = (DX*DX*DY*DY*omega[ix][iy]
			  +DY*DY*psi[ix+1][iy  ]
			  +DY*DY*psi[ix-1][iy  ]
			  +DX*DX*psi[ix  ][iy+1]
			  +DX*DX*psi[ix  ][iy-1])/(2.0*(DX*DX+DY*DY));
	psi[ix][iy] = psi_new;
	if(psiMax < fabs(psi_new))
	  psiMax = psi_new;
	err = (fabs(psi_new - psi_old))/psiMax;
	if(errMax < err)
	  errMax = err;
      }
    }
    if(errMax < EPS) break;
  }
}


/* Output Data */
void output(FILE *fp, Array psi)
{
  int ix, iy;
  for(ix = 0; ix <= N; ix++){
    for(iy = 0; iy <= N; iy++){
      fprintf(fp, "%g %g %g\n", ix*DX, iy*DY, psi[ix][iy]);
    }
    fprintf(fp, "\n");
  }
}

/* Plot Data */
void plot_psi(FILE *gp, Array psi)
{
  int ix, iy;
  fprintf(gp, "set contour\n");
  fprintf(gp, "unset surface\n");
  fprintf(gp, "splot '-' w l\n");
  for(ix = 0; ix <= N; ix++){
    for(iy = 0; iy <= N; iy++){
      fprintf(gp, "%g %g %g\n", ix*DX, iy*DY, psi[ix][iy]);
    }
    fprintf(gp, "\n");
  }
  fprintf(gp, "end\n");
  fflush(gp);
}

/* MAIN */
int main(void)
{
  static Array psi;     // Stream Function
  static Array omega;   // Vorticity
  int  it;
  FILE *fp = fopen("output.dat", "w");
  FILE *gp = popen("gnuplot -persist", "w");

  /* ===Initialize=== */
  initArray(psi);
  initArray(omega);

  /* ===Set Boundary Condition=== */
  setBoundary_psi(psi);
  setBoundary_omega(omega, psi);

  /* ===Solve Omega & Psi=== */
  for(it = 1; it <= TMAX; it++){
    solveOmega(omega, psi);
    solvePsi(psi, omega);
    setBoundary_psi(psi);
    setBoundary_omega(omega, psi);
    printf("Iterator:%d\r", it);
  }
  printf("\n");

  /* ===Output Data=== */
  output(fp, psi);
  fclose(fp);

  /* ===Plot Data=== */
  plot_psi(gp, psi);
  pclose(gp);
  
  return 0;
}

| | コメント (0) | トラックバック (0)

2008/09/06

Box-Muller法で発生させた乱数で正規分布を描く

[Under Construction]

Box-Muller法で正規分布に従う乱数を発生させて、その乱数で実際に正規分布をヒストグラムで描画するプログラムを作る。

これは、平均0、分散1の標準正規分布。

Normaldsitri_2

プログラムは、「distribution.c」

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

#define NT     10000000             // 発生させる乱数の数
#define BOX    500                  // ヒストグラムの解像度
#define UNIT   0.02                 // 刻み幅
#define WIDTH  10                   // WIDTH := BOX*UNIT
#define LEDGE -0.5*WIDTH            // 最小の座標(Left Edge)
#define VALUE  (double)BOX/NT/WIDTH // ヒストグラムの単位量

void   setSeed() {srand((unsigned)time(NULL));}
double genRand() {return rand()/(1.0+RAND_MAX);}

double gaussian(double myu, double sgm)
{
  static int    frag = 0;
  static double save = 0.0;
  
  if(frag == 0){
    double u_0 = genRand();
    double u_1 = genRand();
    double v_0 = myu + sgm*sqrt(-2*log(1-u_0))*cos(2*M_PI*u_1);
    double v_1 = myu + sgm*sqrt(-2*log(1-u_0))*sin(2*M_PI*u_1);
    save       = v_1;
    frag       = 1;
    return v_0;
  }
  else{
    frag = 0;
    return save;
  }  
}

void initArray(double array[], int size)
{
  int it;
  for(it = 0; it < size; it++)
    array[it] = 0.0;
}

void histogram(double hist[], int size, double nRand)
{
  double ibox      = (nRand-LEDGE)/UNIT;
  hist[(int)ibox] += VALUE;
}

double calcMean(const double nRand[], int size)
{
  double sum = 0.0;
  int    it;
  for(it = 0; it < size; it++)
    sum += nRand[it];
  return sum/size;
}

double calcVariance(const double nRand[], int size, double mean)
{
  double sum = 0.0;
  int    it;
  for(it = 0; it < size; it++)
    sum += nRand[it]*nRand[it] - mean*mean;
  return sum/size;
}

void plotter(FILE* gp, const double hist[], int size)
{
  int it;
  fprintf(gp, "set terminal x11\n");
  fprintf(gp, "set title 'Normal Distribution'\n");
  fprintf(gp, "unset key\n");
  fprintf(gp, "plot '-' w l\n");
  for(it = 0; it < size; it++){
    fprintf(gp, "%g %g\n", LEDGE+UNIT*it, hist[it]);
  }
  fprintf(gp, "end\n");
  fflush(gp);
}

int main(int argc, char *argv[])
{
  static double hist[BOX]; // ヒストグラム
  static double nRand[NT]; // 乱数
  int    it;
  double myu    = 0.0; // Box-Muller法で発生させる分布の平均
  double sgm    = 1.0; // Box-Muller法で発生させる分布の標準偏差
  double mean   = 0.0; // 実際に計算する分布の平均
  double vari   = 0.0; // 実際に計算する分布の分散
  FILE   *gp    = popen("gnuplot -persist", "w");
  
  if(2 <= argc && argc <= 3){
    myu = atof(argv[1]);
    sgm = atof(argv[2]);
    sgm = sqrt(sgm);
  }
  else if(argc >= 4){
    fprintf(stderr, "Error! Arguments too many!!\n");
    exit(EXIT_FAILURE);
  }
  setSeed();
  initArray(hist, BOX);
  initArray(nRand, NT);

  for(it = 0; it < NT; it++){
    nRand[it] = gaussian(myu, sgm);
    histogram(hist, BOX, nRand[it]);
  }

  mean = calcMean(nRand, NT);
  vari = calcVariance(nRand, NT, mean);

  plotter(gp, hist, BOX);
  printf("========================================\n");
  printf("Mean = %g, Variance = %g\n", mean, vari);
  printf("========================================\n");
  pclose(gp);

  return 0;
}

引数に平均と分散の値を入力して実行すると、それに従った分布が描画される。(Wikipedia:正規分布)

ヒストグラムを作る関数は工夫しないと実行速度が極端に遅くなる。

| | コメント (0) | トラックバック (0)

2008/07/14

グラフの表示を遅延

Gnuplotのグラフアニメーションをじっくり見るために、時間制限付きキー入力待ちをselect()関数を使って実装しようとしていたけど、Gnuplotのpauseコマンドを使えば表示を遅延させられることを発見した。(なんで始めにGnuplotの機能を調べないんだ・・・)

使い方は簡単。

> pause time

ex)  > pause 0.1

コレだけ。timeのところに遅延させたい時間を記述する。(恐らく単位は秒)

プログラムに組み込むときも、一行で済むから簡単。

実際どんな動きになるかは、下のサンプルプログラムを実行してみるとわかる。

/*
  $ gcc -Wall -O fou1.c -lm
  でコンパイル
*/
#include <stdio.h>
#include <math.h>

/* Parameter */
#define PERIOD    2*M_PI    // 区間
#define EDGE_L   -PERIOD/2  // 区間の左の端点
#define NP        1000      // プロットする点の数
#define DP        PERIOD/NP // 点と点の距離      
#define N_MAX     101       // 重ね合わせる回数

/*
  周期関数
   f(t) = -1 (-pi < t <= 0 )
        =  1 (  0 < t <= pi)
  フーリエ級数展開
   f(t) ~ a_0 + coef * sigma(n=1 -> inf){a_n*cos(var_a) + b_n*sin(var_b)}
   f(t) ~ 1/2 + 2/pi * sigma(n=1 -> inf){sin(2n-1)t/(2n-1)}
  座標データ
   f(t) = pt[NP]
     t  = EDGE_L + DP * ip (0 <= ip < NP)
*/

void func_t(double pt[], int n_max)
{
  int ip, in;
  
  for(ip = 0; ip < NP; ip++){
    double sigma = 0;
    double t     = EDGE_L + DP * ip;
    double a_0   = 1.0/2.0;
    double coef  = 2.0/M_PI;
    for(in = 1; in <= n_max; in++){
      double a_n   = 0;
      double var_a = 0;
      double b_n   = 1.0/(2.0*in-1.0);
      double var_b = (2.0*in-1.0)*t;
      
      sigma += a_n * cos(var_a) + b_n * sin(var_b);
    }
    pt[ip] = a_0 + coef * sigma;
  }
}

void gnuplot_run(FILE *gp, const double pt[], int n_max)
{
  int ip;
  fprintf(gp, "set terminal x11\n"); // gnuplot4.2.2ではwxtが不安定
  fprintf(gp, "unset key\n");
  fprintf(gp, "set title \"n=%d\"\n", n_max);
  fprintf(gp, "plot '-' w l\n");
  for(ip = 0; ip < NP; ip++){
    double t = EDGE_L + DP * ip;
    fprintf(gp, "%g %g\n", t, pt[ip]);
  }
  fprintf(gp, "end\n");
  fprintf(gp, "pause %g\n", 0.05); // pauseコマンドで指定秒数だけ遅延
  fflush(gp);
}
	      
int main(void)
{
  static double pt[NP];
  FILE *gp = popen("gnuplot -persist\n", "w");
  
  int in;
  for(in = 1; in <= N_MAX; in++){
    func_t(pt, in);
    gnuplot_run(gp, pt, in);
  }
  pclose(gp);

  return 0;
}

このプログラムはフーリエ級数が元の周期関数に近づく様子を見るために作ったもの。

pauseコマンドの部分をコメントアウトするとスピードが速すぎてじっくり見れない(マクロN_MAXを11とかにするとそりゃぁもう)けど、pauseをいれておけば様子がじっくり見れる。

他の周期関数が見たい人は自分でフーリエ級数展開しませう。

| | コメント (0) | トラックバック (0)

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

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

| | コメント (0) | トラックバック (0)