4.2.3 ソースプログラム reidai7-1-glsc.c

(情報処理IIの授業で使っていたコンピューター環境は古いので、 当時使っていたプログラムを GLSC を利用するように書き換えたもの。)

/*
 * reidai7-1-glsc.c --- reidai7-1.c を GLSC を用いるように書き換えた
 * http://nalab.mind.meiji.ac.jp/~mk/program/ode_prog/reidai7-1-glsc.c
 */

/**********************************************************************
 **********************************************************************
 *
 *   2 次元の定数係数線形常微分方程式
 *       x'(t) = a x + b y
 *       y'(t) = c x + d y
 *  に初期条件
 *       x(0)=x0
 *       y(0)=y0
 *  を課した常微分方程式の初期値問題を解いて、相図を描く。
 *
 *  このプログラムは次の4つの部分から出来ている。
 *     main()
 *       主プログラム
 *       行列の係数の入力、ウィンドウのオープン等の初期化をした後に、
 *       ユーザーにメニュー形式でコマンドを入力してもらう。
 *       実際の計算のほとんどは他の副プログラムに任せている。
 *     draworbit(x0,y0,h,tlimit)
 *       (x0,y0) を初期値とする解の軌道を描く。
 *       刻み幅を h、追跡時間を tlimit とする。
 *       近似解の計算には Runge-Kutta 法を用いる。
 *     fx(x,y)
 *       微分方程式の右辺の x 成分
 *     fy(x,y)
 *       微分方程式の右辺の y 成分
 **********************************************************************/

#include <stdio.h>
#include <stdlib.h> // system()
#include <math.h>
#include <glsc.h>

// 係数行列 A の成分
double a, b, c, d;
// ウィンドウに表示する範囲
double xleft, xright, ybottom, ytop;

void draworbit(double, double, double, double);

void g_dump(char *fname, Display *display, Window wid)
{
  char command[256];
  sprintf(command, "import -window %lu %s", wid, fname);
  system(command);
}

int main(void)
{
  // 初期値
  double x0, y0;
  // 時間の刻み幅、追跡時間
  double h, tlimit, newh, newT;
  // メニューに対して入力されるコマンドの番号
  int cmd;
  // マウスのボタンの状態
  int but;
  //
  double win_width, win_height, w_margin, h_margin;
  //
  Display *display;
  Window wid; // Window ID
  // ウィンドウに表示する範囲の設定
  xleft = -1.0;
  xright = 1.0;
  ybottom = -1.0;
  ytop = 1.0;
  // 時間刻み幅、追跡時間（とりあえず設定）
  h = 0.01;
  tlimit = 10.0;
  // 行列の成分を入力
  printf(" a,b,c,d=");
  scanf("%lf%lf%lf%lf", &a, &b, &c, &d);
  // ウィンドウを開く
  win_width = 200.0; win_height = 200.0; w_margin = 10.0; h_margin = 10.0;
  g_init("GRAPH", win_width  + 2 * w_margin, win_height + 2 * h_margin);
  g_device(G_BOTH);
  g_def_scale(0, xleft, xright, ybottom, ytop,
              w_margin, h_margin, win_width, win_height);
  g_sel_scale(0);
  // x 軸、y 軸を描く
  g_move(xleft, 0.0); g_plot(xright, 0.0);
  g_move(0.0, ybottom); g_plot(0.0, ytop);
  // メイン・ループの入口
  do {
    // メニューを表示して、何をするか、番号で選択してもらう
    printf(" したいことを番号で選んで下さい。\n");
    printf("  -1:メニュー終了, 0:初期値のキーボード入力, 1:初期値のマウス入力\n");
    printf("   2: 刻み幅,追跡時間変更(現在 h=%7.4f, T=%7.4f)\n", h,tlimit);
    scanf("%d", &cmd);
    // 番号 cmd に応じて、指示された仕事をする
    if (cmd == 0) {
      // 初期値の入力
      printf(" 初期値 x0,y0=");
      scanf("%lf%lf", &x0, &y0);
      draworbit(x0,y0,h,tlimit);
    }
    else if (cmd == 1) {
      do {
        printf("マウスの左ボタンで初期値を指定して下さい（右ボタンで中止）。\n");
        g_mouse_sence(&x0, &y0, &but);
        printf("x0=%g, y0=%g, but=%d\n", x0, y0, but);
        if (but == 1) {
          printf("左ボタン, (x0,y0)=%f %f\n", x0,y0);
          draworbit(x0,y0,h,tlimit);
        }
        else if (but == 2) {
          printf("中ボタン, (x0,y0)=%f %f\n", x0,y0);
          draworbit(x0,y0,-h,tlimit);
        }
      }
      while (but != 3);
      printf("右ボタンがクリックされました。マウスによる初期値の入力を打ち切ります。\n");
    }
    else if (cmd == 2) {
      // 時間刻み幅、追跡時間の変更
      printf("時間刻み幅 h (%g), 追跡時間 T (%g): ", h, tlimit);
      scanf("%lf%lf", &newh, &newT);
      if (newh != 0 && newT != 0) {
        h = newh;
        tlimit = newT;
        printf("新しい時間刻み幅 h = %g, 新しい追跡時間 T = %g\n", h, tlimit);
      }
      else {
        printf("h=%g, T=%g は不適当です。\n", newh, newT);
      }
    }
  }
  while (cmd != -1);

  g_dump("reidai7-1.png", g_get_display(), g_get_window());

  printf("GLSCウィンドウを左ボタンでクリックして下さい\n");
  g_sleep(-1.0);
}
    
// 指示された初期値に対する解軌道を描く
void draworbit(double x0, double y0, double h, double tlimit)
{
  int in;
  double x, y, fx(double, double), fy(double, double);
  double k1x, k1y, k2x, k2y, k3x, k3y, k4x, k4y, t;
  // 時刻を 0 にセットする
  t = 0.0;
  // 初期値のセット
  x = x0;
  y = y0;
  // 初期点を描く
  if ((in = (xleft <= x && x <= xright && ybottom <= y && y <= ytop)) != 0)
    g_move(x,y);
  // ループの入口
  do {
    // Runge-Kutta 法による計算
    // k1 の計算
    k1x = h * fx(x, y);
    k1y = h * fy(x, y);
    // k2 の計算
    k2x = h * fx(x + k1x / 2.0, y + k1y / 2.0);
    k2y = h * fy(x + k1x / 2.0, y + k1y / 2.0);
    // k3 の計算
    k3x = h * fx(x + k2x / 2.0, y + k2y / 2.0);
    k3y = h * fy(x + k2x / 2.0, y + k2y / 2.0);
    // k4 の計算
    k4x = h * fx(x + k3x, y + k3y);
    k4y = h * fy(x + k3x, y + k3y);
    // (Xn+1,Yn+1) の計算
    x += (k1x + 2.0 * k2x + 2.0 * k3x + k4x) / 6.0;
    y += (k1y + 2.0 * k2y + 2.0 * k3y + k4y) / 6.0;
     // 解軌道を延ばす
    if (in) {
      if ((in = (xleft <= x && x <= xright && ybottom <= y && y <= ytop)) != 0)
        g_plot(x,y);
    }
    else {
      if ((in = (xleft <= x && x <= xright && ybottom <= y && y <= ytop)) != 0)
        g_move(x,y);
    }
    // 時刻を 1 ステップ分進める
    t += h;
  }
  while (fabs(t) <= fabs(tlimit)); // まだ範囲内かどうかチェック
}

// 微分方程式の右辺のベクトル値関数 f の x 成分
double fx(double x, double y)
{
  return a * x + b * y;
}

// 微分方程式の右辺のベクトル値関数 f の y 成分
double fy(double x, double y)
{
  return c * x + d * y;
}