課題4プログラム例

課題４−１

微分方程式の数値解法関数を別ファイルにする

ここで、プログラム開発の標準的なルールを説明します。

• プログラムは長くなると読みにくくなるので、プログラム単位ごとに別のファイルにし、コンパイルするときに統合します。 (100行を超えるようになれば分割を考えよう。)
• 分割した時には、関数定義はヘッダーファイルに書き、include文でそれらを読み込みます。ヘッダーファイルの拡張子は".h"です。

以下、汎用の微分方程式数値計算関数を別ファイルにする例を示します。

ヘッダーファイルはdiffEqMethod.h、本体の記述部分をdiffEqMethods.cとします。(名前の付け方は自由ですが、中身がわかるように名付けることが重要です。)

これらを使うprac4.cをコンパイルするには、コマンドラインにソースファイル名を並べます。ヘッダーファイルのインクルードはそれぞれのソースの最初の方に書き込みます。

cc -o prac4 prac4.c diffEqMethods.c -lm -lcpgplot -lpgplot -lX11

diffEqMethods.hの例

/* 微分方程式の数値計算関数の定義
 *  通常は関数ごとに引数のコメントをきちんと書く
 */
void euler(int n, double *xin, double *xout, double h, void (*f)(double *, double *, double *), double *params);
void rk4(int size_of_x, double *x0, double *x, double h, void (*f)(), double *params);

diffEqMethods.cの例

/* 微分方程式の数値解法関数 */
#include <stdlib.h>
#include "diffEqMethods.h"

/* オイラー法 */
void euler(int size_of_x, double *x0, double *x, double h, void (*f)(double *, double *, double *), double *params) {

  int i;
  double *funcVal;

  funcVal = (double *) malloc(size_of_x * sizeof(double));

  f(x0, funcVal, params);  // 関数値の計算 (入力x0に対する微分方程式の右辺の値をfuncValに受け取る
  for(i=0; i<size_of_x; i++) {
    x[i] = x0[i] + h * funcVal[i]; //オイラー法
  }

}

/* ルンゲ・クッタ法 */
void rk4(int size_of_x, double *x0, double *x, double h,  void (*f)(), double *params) {

  int i;
  double *funcVal[4];

  for(i=0;i<4;i++) {
   funcVal[i] = (double *) malloc(size_of_x * sizeof(double));
  }
  f(x0, funcVal[0], params);
  for ( i = 0 ; i < size_of_x ; i++)
    x[i] = x0[i] + h/2.0*funcVal[0][i];
  f(x, funcVal[1], params);
  for ( i = 0 ; i < size_of_x ; i++)
    x[i] = x0[i] + h/2.0* funcVal[1][i];
  f(x, funcVal[2],params);
  for ( i = 0 ; i < size_of_x ; i++)
    x[i] = x0[i] + h* funcVal[2][i];
  f(x, funcVal[3], params);

  for ( i = 0 ; i < size_of_x ; i++)
    x[i] = x0[i] + h/6.0 * (funcVal[0][i]+2.0*funcVal[1][i]+2.0*funcVal[2][i]+funcVal[3][i]);

}

課題４−1のプログラム例

上記のdiffEqMethods.cはそのまま使う．func()を別のファイルにしてしまうこともあり得る．あるいは，func()をdiffEqMethods.cに入れる構成でもよい．

/* 交通流のOVモデル(Optimal Velocity Model) */

#include "cpgplot.h"
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include "diffEqMethods.h"

#define N 80 /* 車の数 */
#define DT  0.002
#define UpdateMax 200000.0

void rk4(int size_of_x, double *x0, double *x, double h, void (*f)(),  double *params);
void func(double *x, double *funcVal, double *params); /* xは現在の値（入力）、funcValは返値（微分方程式の右辺）*/

double L=1000.0;
double c=2.0, tanhc;    /* 最適速度に含まれる係数*/



int main(int ARGC, char *ARGV[]) {

  double x[2*N+2], xold[2*N+2], xold4draw[2*N+2];
  /* 車の位置,座標を記憶する配列 0〜Nまでが位置，N+1〜2N+1までが速度*/
  /* xold4drawは、描画用の「以前の位置」 */
  long n = 0, m;
  double t=0, told=0;
  double params[1];

  tanhc = tanh(c);

  printf("道路の長さは?(%4.0f程度) ", c*N); scanf("%lf", &L);
  printf("感度aは? (臨界速度=%.3f)", 2./cosh(L/N-2.0)); scanf("%lf", &params[0]);

  cpgopen("/xserv");         /* Xwindow(Xサーバ)上へ描画 */
  cpgpap(12.0,0.6);    /* ウィンドウの大きさと縦横比 */
  cpgenv(0.0, (float) UpdateMax*DT, 0.0,  (float)L, 0,0);
                          /* 座標の設設定 */
  srand(323211);  /* 乱数の種 (適当な整数) */

  /* 初期位置と速度を乱数で構成 */
  for(n = 0; n < N; n++) {
    xold[N+1+n] = (double)rand()/RAND_MAX * tanh(2.0);
    xold[n] = L * (double)n/N ;  /* 初期位置は等間隔 */
  }
  xold[N] = xold[0];
  xold[2*N+1] = xold[N+1];
  memcpy(x, xold, sizeof(x));
  memcpy(xold4draw, xold, sizeof(x));

  /* update */
  while (t < UpdateMax*DT) {
	  told = t;
	  for(n=0; n< 100; n++) {
		  /* 時間発展 */
	    rk4(2*N+2, xold, x, DT, func, params);
		  memcpy(xold, x, sizeof(x));

		  for(m=0; m<N; m++) {
			  if(xold[m]> L) xold[m] -= L; // 周期境界条件
		  }

		  xold[N]=xold[0];
		  t+=DT;
	  }

	  /* 描画 */
	  cpgsci(2);
	  cpgslw(2);
	  for (n=0; n< N; n++) {
		  if(xold4draw[n] < x[n]) { // 境界を超えたときは書かない
			  cpgmove((float) told, xold4draw[n]);
			  cpgdraw((float) t, x[n]);
		  }
	  }
	  memcpy(xold4draw, x, sizeof(x)); // 描画した前の位置を覚えておく
  }
  cpgclos();
  exit(0);
}

void func(double *x, double *funcVal, double *params) {
  long i;
  double headway;
  double sensitivity = params[0];

  for(i = 0; i < N; i++) {
    headway = x[i+1] - x[i];
    if(headway < 0) headway += L;
    funcVal[i]  = x[N+1+i];
    funcVal[N+1+i] = sensitivity * (tanh( headway - 2.0) + tanh(2.0)
			  - x[N+1+i]);
  }
  /* 周期境界条件 */
  funcVal[N] = funcVal[0];
  funcVal[2*N+1] = funcVal[N+1];
}

コマンドライン引数としてパラメータを読み込む方法

  printf("道路の長さは?(%4.0f程度) ", c*N); scanf("%lf", &L);
  printf("感度aは? (臨界速度=%.3f)", 2./cosh(L/N-2.0)); scanf("%lf", &params[0]);

  if (ARGC != 3) {
    printf("コマンドライン引数としてLとaを与えてください。(典型例: 160 1.5)\n");
    exit(0);
  }
  L = strtof(ARGV[1],NULL);
  params[0] = strtof(ARGV[2], NULL);
  printf("L=%f, a=%f\n", L, params[0]);

コマンドライン引数は文字列として受け取るので、数値(flot)に変換する。変換にはatofという関数が用いられるが、エラー処理の安全性のためにはstrtofが良いという解説もあるので、ここではstrtofを使ってみた。(atof(ARGV[0])のような記述でも良い。)

課題4-3について： 数値シミュレーションとして行うこと

こちらに手順を書いた。

複数のパラメータがある場合は、それぞれどれか一つのパラメータに対する注目して、値をすこしずつかえ変化の様子を観察することが重要である。

数値データとしてまず考えられるのは「平均速度」である。

• 密度（道路長と台数）を固定して、感度パラメータaを少しずつ変えて流れの様子を観察する。また、平均速度とaの関係をプロットしてみる。
• aと台数を固定して、道路長Lを変えて密度変化による流れの違いや平均速度がどう変化するかを調べてみる。
• a, Lに対する平均速度の変化を図示する鳥瞰図をつくってみてもおもしろい。
• 交通システムとしては、ある点を単位時間に流れる台数、つまり流量と密度の関係である。これを交通工学の分野では基本図とよぶ。
• 渋滞クラスターが生ずる場合、つまり、車間距離が狭い領域と広い領域の２つに別れる場合の両者の特徴を定量的に捉えることも重要である。速度の分布（ヒストグラム）や車間距離の分布など。