F..2 問3

「Cからgnuplotを呼び出す」

toi3.c

/*
 * toi3.c --- euler1ex1.c (Euler method for Malthusian model)
 */

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

double k = 1.0;

int main(void)
{
  int i, N;
  double t, x, dt;
  double f(double, double), x0;
  double Tmax;
  FILE *of;

  if ((of = fopen("euler1ex1.data", "w")) == NULL) {
    fprintf(stderr, "cannot open euler1ex1.data");
    exit(1);
  }
  // 初期値
  x0 = 1.0;
  // 最終時刻
  Tmax = 1.0;
  // 時間刻み
  printf("# N: "); scanf("%d", &N);
  fprintf(of, "# N = %d\n", N);
  dt = Tmax / N;
  // 初期値
  t = 0.0;
  x = x0;
  printf("# t x\n");
  fprintf(of, "# t x\n");
  printf("%f %f\n", t, x);
  fprintf(of, "%f %f\n", t, x);
  // Euler 法
  for (i = 0; i < N; i++) {
    x += dt * f(x, t);
    t = (i + 1) * dt;
    printf("%f %f\n", t, x);
    fprintf(of, "%f %f\n", t, x);
  }
  fclose(of);
  if ((of = popen("gnuplot", "w")) == NULL) {
    fprintf(stderr, "cannot open pipe.\n");
    exit(1);
  }
  fprintf(of, "plot \"euler1ex1.data\" with l\n");
  fflush(of);
  pclose(of);
  return 0;
}

double f(double x, double t)
{
  return k * x;
}

もちろん gnuplot は事前にインストールしておく必要がある。 コンパイル＆実行すると、N の値を尋ねてくるので、 例えば 100 と入力すると、 gnuplot が起動されてグラフが描かれるはず。

popen(), pclose() は、C++ でも使えるので (別の関数を使うべきか？)、 次でとりあえず動く。

toi3.cpp

/*
 * toi3.cpp --- euler1ex1.c (Euler method for Malthusian model)
 */

#include <iostream>
#include <fstream>
using namespace std;

double k = 1.0;

int main(void)
{
  int i, N;
  double t, x, dt;
  double f(double, double), x0;
  double Tmax;
  ofstream ofs("euler1ex1.data");
  FILE *of;
  if (!ofs) {
    cerr << "cannot open euler1ex1.data" << endl;
    exit(1);
  }
  // 初期値
  x0 = 1.0;
  // 最終時刻
  Tmax = 1.0;
  // 時間刻み
  cout << "# N: "; cin >> N;
  ofs << "# N = " << N << endl;
  dt = Tmax / N;
  // 初期値
  t = 0.0;
  x = x0;
  cout << "# t x" << endl;
  ofs  << "# t x" << endl;
  cout << fixed;
  cout << t << " " << x << endl;
  ofs  << fixed;
  ofs << t << " " << x << endl;
  // Euler 法
  for (i = 0; i < N; i++) {
    x += dt * f(x, t);
    t = (i + 1) * dt;
    cout << t << " " << x << endl;
    ofs  << t << " " << x << endl;
  }
  ofs.close();

  if ((of = popen("gnuplot", "w")) == NULL) {
    fprintf(stderr, "cannot open pipe.\n");
    exit(1);
  }
  fprintf(of, "plot \"euler1ex1.data\" with l\n");
  fflush(of);
  pclose(of);

  return 0;
}

double f(double x, double t)
{
  return k * x;
}