C

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/10

NTEmacs + MinGW + MSYSでWindows上にC言語開発環境を構築してみる - その5

前回までの作業でC言語開発環境は整いましたが、実は問題があります。

その問題は次のようなコードを実行してみるとわかります。

#include <stdio.h>

int main(void)
{
  int num;
  
  printf("Please Input Num:");
  scanf("%d", &num);
  printf("Num:%d\n", num);
  
  return 0;
}

これをMSYSから実行しても次のようになってしまいます。

$ ./a.exe
100
Please Input Num:Num:100

scanf()の直前のprintf()で表示するはずの文字列が表示されません。これはMSYSで使っているrxvtの実装に問題があるためです。

この問題の対処法はrxvtを使用せず、Windows標準のコマンドプロンプトを使用するのが簡単です。具体的に「はC:\msys\1.0」内にある「msys.bat」の41行目の

if EXIST rxvt.exe goto startrxvt

をコメントアウトして

rem if EXIST rxvt.exe goto startrxvt

とすればMSYS起動時にrxvtが呼び出されなくなります。

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

2008/10/09

NTEmacs + MinGW + MSYSでWindows上にC言語開発環境を構築してみる - その4

1.動作確認

まずはデスクトップ上に作成されたMSYSのショートカットをダブルクリックしてMSYSを起動して次のようにコマンドを入力してください。

user@host ~
$ pwd
/c/Home

user@host ~
$ gcc
gcc.exe: no input files

user@host ~
$ make
make: *** No targets specified and no makefile found.  Stop.

user@host ~
$ gdb
GNU gdb 6.3
Copyright 2004 Free Software Foundation, Inc.
GDB is free software, covered by the GNU General Public License, and you are
welcome to change it and/or distribute copies of it under certain conditions.
Type "show copying" to see the conditions.
There is absolutely no warranty for GDB.  Type "show warranty" for details.
This GDB was configured as "i686-pc-mingw32".
(gdb) quit

user@host ~
$ 

gcc、make、gdbはこのようなメッセージが表示されればOKです。また、pwdコマンドの実行結果のとおりホームディレクトリは「C:\Home」になっています。MSYSのデフォルトのホームディレクトリはC:\msys\1.0\homeなのですが、環境変数HOMEが設定されているとHOMEで指定されているパスをホームディレクトリを認識します。

2.設定(お好みで)

MSYSを起動してみてわかると思いますが、MSYSのシェルは白地にグリーンとブルーのプロンプトで目がチカチカします。そこで色を変更してみます。

シェルの色を変更するには「C\msys\1.0」内にある「msys.bat」を編集します。

55行目にある
if "x%MINGW32BGCOLOR%" == "x" set MINGW32BGCOLOR=LightYellow
if "x%MINGW32FGCOLOR%" == "x" set MINGW32FGCOLOR=Navy
を次のようにコメントアウトしてこのような2行を追加
rem if "x%MINGW32BGCOLOR%" == "x" set MINGW32BGCOLOR=LightYellow
rem if "x%MINGW32FGCOLOR%" == "x" set MINGW32FGCOLOR=Navy
if "x%MINGW32BGCOLOR%" == "x" set MINGW32BGCOLOR=Black
if "x%MINGW32FGCOLOR%" == "x" set MINGW32FGCOLOR=White

また、Emacsの背景も黒にしたい場合はemacsコマンドに-rvオプションをつけて起動すればEmacsが黒基調になりますが、いつもオプションをつけるのは面倒なので'emacs -rv'を'emacs'と認識するようにエイリアスを登録しておきます。

エイリアスを登録するには「C:\msys\1.0\etc」内にある「profile」に

alias emacs='emacs -rv'

と追記しておく。(MSYSを再起動すると変更が反映される)

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

2008/10/08

NTEmacs + MinGW + MSYSでWindows上にC言語開発環境を構築してみる - その3

今度はMinGWとMSYSをインストール。

1.MinGWのインストール

まずはSourceForgeからMinGW-5.1.4.exeをダウンロードします。ダウンロードしたらインストーラを起動します。(このとき新しくフォルダを作って、その中でインストーラを起動したほうがいい)

インストール途中にインストールするコンパイラの選択画面が表示されるので欲しいインストーラを選択します。その後インストール先の指定画面が表示されますが、もしパスに空白が含まれているとトラブルの元なのでインストール先は「C:\MinGW」とします。

2.MSYSのインストール

SourceForgeからMSYS-1.0.10.exeをダウンロードしてきて、インストーラを起動します。(MinGWのインストーラと同じフォルダ内で起動する)

インストール先を聞かれたら「C:\msys\1.0」とします。

インストールが開始されてプログレスバーが100%になるとコマンドプロンプトが立ち上がり、いくつかの質問に答えることになります。

Q.1 Do you wish to continue with the post install? [yn ] ・・・ y

Q.2 Do you have MinGW installed? [yn ] ・・・ y

Q.3 Please answer the following in the from of c:/foo/bar.
       Where is your MinGW installation?
MinGWをインストールしたパスを入力します。このとき「\」ではなく「/」を使います。

このあと、もし「Makeがありません」といった内容のメッセージが表示されても問題ないので先に進んでください。

3.GDBのインストール

一応、デバッガをインストールしておきます。

SourceForgeからgdb-6.3-2.exeをダウンロードしてきてインストーラを起動します。インストール先はMinGWをインストールしたパスと同じに。「C:\MinGW」

4.環境変数の設定

PATHに、「C:\MinGW\bin;C:\msys\1.0\bin」を追加

また以下を新規作成

C_INCLUDE_PATH=C:\MinGW\include

CPLUS_INCLUDE_PATH=C:\MinGW\include

LIBRARY_PATH=C:\MinGW\lib

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

2008/10/07

NTEmacs + MinGW + MSYSでWindows上にC言語開発環境を構築してみる - その2

まずはNTEmacsのインストール。

1.NTEmacsのバイナリをダウンロードする。

NTEmacsはIMEを利用した日本語入力に難があるようなので、「NTEmacs JP Project」からIME patchが当てられたNTEmacsのバイナリをダウンロードします。(現在の最新版は22.2)

Link:NTEmacs JP Project - SourceForge.JP

2.インストール

ダウンロードしたバイナリを起動すると解凍する場所を聞いてくるので、「C:\」を指定します。(ハードディスクのルートならC:\でなくてもよい)

解凍が終わるとC:\emacsが作成されます。NTEmacsの実行ファイルは C:\emacs\22.2\bin内のrunemacs.exeです。

3.起動する前に

NTEmacsを起動する前にいくつか設定を済ませておきます。

まずはC:\に「Home」というフォルダを作成しておく。(今後、このフォルダがホームディレクトリになる)

次は環境変数PATHに、NTEmacsの実行ファイルが置かれているパス「C:\emacs\22.2\bin」を追加する(パスとパスはセミコロン「;」で区切る)。これでNTEmacsをどこからでも呼び出すことができるようになる。

さらに環境変数HOMEを作って、値には「C:\Home」を設定する。

4.NTEmacsを起動

とりあえずrunemacs.exeのショートカットをデスクトップに作成してNTEmacsを起動してみます。

Ntemacs_2 

設定ファイル「.emacs」はC:\Home内に作成します。.emacsを新規作成して

(setq inhibit-startup-message t)

と書いて保存し、NTEmacsを再起動すると先ほどの起動画面が表示されなくなることが確認できるはずです。

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

2008/10/06

NTEmacs + MinGW + MSYSでWindows上にC言語開発環境を構築してみる - その1

Windows上にC言語開発環境を構築する方法はいろいろあるが、比較的よく見かけるのは以下のようなものだろう。

  • Borland C++ CompilerもしくはVisual C++ Express Editionをインストールしてメモ帳等のエディタでコーディング
  • CygwinをインストールしてWindows上にUNIXライクな環境を構築する
  • Eclipse CDT + MinGW
  • coLinuxもしくはandLinuxを使う
  • VMware等の仮想PCソフトを使ってLinuxをエミュレート(Windows上ではないけど)

私はcoLinux/andLinux以外は全て試してみて現在はVMwareに落ち着いている。

しかし、この方法はLinuxディストリビューションをまるごと使いたいという目的のほうが強いので、単純にC言語開発環境を構築するという観点で見ると無駄が大きい。(しかもある程度のパフォーマンスを得るにはCPUやメモリ容量などの面で制約があるので一昔前のPCだと厳しい)

これはLinuxカーネルをWindows上で動作させるcoLinux/andLinuxにも言えることだと思う。

そこで今回は次のような方針で開発環境を構築する。

  • できるだけシンプルに
  • C言語入門レベルに統合開発環境はいらない
  • あまりマウスはつかいたくない
  • なるべくUNIXライク
  • Emacs使いたい

「Cygwinで全て解決」な感じもするが、CygwinもC言語開発環境としては十分大げさであることと、いつもCygwinではつまらないので、NTEmacs + MinGW + MSYSで構築してみる。(VistaでCygwin/Xが不安定であることも理由の一つ)

| | コメント (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/28

2次元ランダムウォーク

2次元ランダムウォークをプログラムにして可視化してみた。(ブラウン運動に関係しているので)

Link:2次元ランダムウォーク - Wikipedia

まず点を座標(x, y)=(0, 0)に置きます。次の瞬間この点がいる場所は確率1/4で以下のいずれかの場所にいることになる。

  • (x+1, y)
  • (x-1, y)
  • (x , y+1)
  • (x , y-1)

以上を任意のステップ数繰り返すだけ。例えばこのような軌跡を描く。

Randomwalk

プログラムは「randwalk.c」 。(そんなにいい乱数じゃないけどまぁいいや)

/*
  2008/07/28
   2次元ランダムウォーク
*/

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

#define N 100000

/* 乱数生成 */
void   setSeed() {srand(time(NULL));}
double genRand() {return rand()/(1.0+RAND_MAX);}

/* 座標 */
typedef struct {
  double x;
  double y;
} coordinate;

/* Plot Data */
void runGnuplot(FILE* gp)
{
  fprintf(gp, "set terminal x11\n");
  fprintf(gp, "unset key\n");
  fprintf(gp, "set title 'Random Walk'\n");
  fprintf(gp, "plot 'randwalk.data' w l\n");
}

int main(void)
{
  coordinate point = {0.0, 0.0};
  FILE*      fp    = fopen("randwalk.data", "w");
  FILE*      gp    = popen("gnuplot -persist", "w");
  int        i;

  setSeed();  // Set Random Seed
  
  for(i = 0; i < N; i++){
    double prob = genRand();
    if(prob < 0.25)
      point.x += 1.0;
    else if(prob >= 0.25 && prob < 0.5)
      point.x -= 1.0;
    else if(prob >= 0.5  && prob < 0.75)
      point.y += 1.0;
    else
      point.y -= 1.0;
    
    fprintf(fp, "%g, %g\n", point.x, point.y);
  }
  fclose(fp);

  runGnuplot(gp);
  pclose(gp);
  
  return 0;
}

| | コメント (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)

より以前の記事一覧