9.1.2.1 例: 複数のグラフを plot する

外部ファイルと異なり、'-' は読み直しが出来ないので、 データは一つ一つ出力する (対応する 'e' もグラフの個数分必要になる) というのがミソです。

3sukumi2.c

/*
 * 3sukumi2.c --- 3すくみ方程式の数値シミュレーション
 *   http://phi.med.gunma-u.ac.jp/oldlec/ecology_p14.html の 3sukumi.R を
 *   C 言語に焼き直し
 */

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

#define N (200)

void runif(double *v, int n, double a, double b)
{
  int i;
  for (i = 0; i < n; i++) v[i] = a + (b - a) * drand48();
}

int main(void)
{
  double N1[N+1],N2[N+1],N3[N],yr[N];
  double cc[3] = {0.01, 0.01, 0.01};
  double a[] = {0.2, 0.4, 0.1, 0.1, 0.2, 0.4, 0.4, 0.1, 0.2};
  double r[] = {0.5, 0.6, 0.7};
  double rand1[N], rand2[N], rand3[N];
  int i;
  FILE *gp;
  double hw;

  N1[0] = 1; N2[0] = 1; N3[0] = 1;
  /* yr = {1,2,...,N} */
  for (i = 0; i < N; i++)
    yr[i] = i + 1;
  /* drand48() の seed設定 --- 適当 */
  srand48(1L);
  /* [-0.1,0.1) の一様乱数 N 個 */
  hw = 0.1;
  runif(rand1, N, -hw, hw); runif(rand2, N, -hw, hw); runif(rand3, N, -hw, hw);
  /* 微分方程式を Euler 法で解く */
  for (i = 0; i < N-1; i++) {
    N1[i+1] = cc[0] + exp(r[0]+rand1[i]-a[0]*N1[i]-a[1]*N2[i]-a[2]*N3[i])*N1[i];
    N2[i+1] = cc[1] + exp(r[1]+rand2[i]-a[3]*N1[i]-a[4]*N2[i]-a[5]*N3[i])*N2[i];
    N3[i+1] = cc[2] + exp(r[2]+rand3[i]-a[6]*N1[i]-a[7]*N2[i]-a[8]*N3[i])*N3[i];
  }
  /* gnuplot を呼び出して描画 */
  gp = popen("gnuplot -persist", "w");
  fprintf(gp, "plot '-' with lines, '-' with lines, '-' with lines\n");
  for (i = 0; i < N; i++) fprintf(gp, "%f %f\n", yr[i], N1[i]); fprintf(gp, "e\n");
  for (i = 0; i < N; i++) fprintf(gp, "%f %f\n", yr[i], N2[i]); fprintf(gp, "e\n");
  for (i = 0; i < N; i++) fprintf(gp, "%f %f\n", yr[i], N3[i]); fprintf(gp, "e\n");
  fprintf(gp, "quit\n");
  pclose(gp);
  return 0;
}

図 19: 3sukumi2.c の実行結果

(2021/4/28追記) 上の 3sukumi.c は元ネタが自分で書いたものでなくて、 どうも良く分からないところがあるので、別のサンプル・プログラムを用意した。

lotka-volterra-rk2.c

/*
 * lotka-volterra-rk2.c (Runge-Kutta method for Lotka-Volterra model)
 * ターミナルで
 *   cc lotka-volterra-rk2.c
 *   ./a.out > lotka-volterra-rk2.data
 *   500
 *   gnuplot
 * gnuplot> に対して
 *  plot "lotka-volterra-rk2.data" using 2:3 with l
 */

#include <stdio.h>

double a,b,c,d;

int main(void)
{
  int i, N;
  double t, x, y, dt, k1x, k1y, k2x, k2y, k3x, k3y, k4x, k4y;
  double fx(double, double, double), fy(double, double, double), x0, y0;
  double Tmax;
  // パラメーター決定
  a = 2; b = 1; c = 3; d = 1;
  // 最終時刻
  Tmax = 10.0;
  // 時間刻み
  printf("# N: "); scanf("%d", &N);
  dt = Tmax / N;
  // 初期値
  x0 = 3.0;
  for (y0 = 0.5; y0 <= 5.0; y0 += 0.5) {
    // 初期値
    t = 0.0;
    x = x0;
    y = y0;
    printf("# t x y\n");
    printf("%f %f %f\n", t, x, y);
    // Euler 法
    for (i = 0; i < N; i++) {
      k1x = dt * fx(x, y, t);
      k1y = dt * fy(x, y, t);
      k2x = dt * fx(x + k1x / 2, y + k1y / 2, t + dt / 2);
      k2y = dt * fy(x + k1x / 2, y + k1y / 2, t + dt / 2);
      k3x = dt * fx(x + k2x / 2, y + k2y / 2, t + dt / 2);
      k3y = dt * fy(x + k2x / 2, y + k2y / 2, t + dt / 2);
      k4x = dt * fx(x + k3x, y + k3y, t + dt);
      k4y = dt * fy(x + k3x, y + k3y, t + dt);
      x += (k1x + 2 * k2x + 2 * k3x + k4x) / 6;
      y += (k1y + 2 * k2y + 2 * k3y + k4y) / 6;
      t = (i + 1) * dt;
      printf("%f %f %f\n", t, x, y);
    }
    printf("e\n\n"); // 線を切る
  }
  return 0;
}

double fx(double x, double y, double t)
{
  return a * x - b * x * y;
}

double fy(double x, double y, double t)
{
  return - c * y + d * x * y;
}

コンパイル、実行、可視化

% cc lotka-volterra-rk2.c
% ./a.out > lotka-volterra-rk2.data
500
% gnuplot
gnuplot> plot "lotka-volterra-rk2.data" using 2:3 with l

図: Lotka-Volterraの方程式の解軌道