next up previous contents
Next: 2.3.2 Runge-Kutta法のプログラム例 Up: 2.3 C言語のプログラム例 Previous: 2.3 C言語のプログラム例

2.3.1 Euler法のCプログラム例

euler1ex1.c

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

#include <stdio.h>

double k = 1.0;

int main(void)
{
  int i, N;
  double t, x, dt;
  double f(double, double), x0;
  double Tmax;
  // 初期値
  x0 = 1.0;
  // 最終時刻
  Tmax = 1.0;
  // 時間刻み
  printf("# N: "); scanf("%d", &N);
  dt = Tmax / N;
  // 初期値
  t = 0.0;
  x = x0;
  printf("# t x\n");
  printf("%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);
  }
  return 0;
}

double f(double x, double t)
{
  return k * x;
}
cc でコンパイル&実行
% cc euler1ex1.c
% ./a.out
# N: 100
# t x
0.000000 1.000000
中略
1.000000 2.704814
%

cglsc でコンパイル＆実行
% cglsc euler1ex1.c
% ./euler1ex1
# N: 100
(以下略)
%

解曲線を描いてみよう。この場合は gnuplot が簡単である。
解曲線を描く
% cc euler1ex1.c
% ./a.out > euler1ex1.data
100
% gnuplot
(ここで gnuplot の起動メッセージが表示されるが省略)
gnuplot> plot "euler1ex1.data" with lp
(これでグラフが表示される。以下は画像の保存。)
gnuplot> set term png
gnuplot> set output "euler1ex1.png"
gnuplot> replot
gnuplot> quit
%

図: $ \D x/\D t=x$, $ x(0)=1$ に対する Euler法の解曲線
Image euler1ex1

gnuplot については、ネットに解説があふれている (例えば桂田 [7])。 この問題について役に立ちそうなことをいくつか説明しておく。

gnuplot の細かな工夫に関する注意

問1      上の実行例で最後に出力した数値 2.704814 の誤差はいくらか。 $ N$ を変えることで誤差がどのように変わるか調べよ (両側対数目盛でプロットすることを勧める)。

問2      上の実行例では、出力をリダイレクトすることでデータファイルを作っているが、 fopen(), fclose(), fprintf() を使ってデータファイルを作れ (参考 「簡単なファイル入出力」, 解答は付録 F.1 を見よ)。

問3      Cのプログラムの中で popen() という関数を使うと、 外部プログラムを起動して通信ができる。 gnuplot を起動して解曲線を描く C プログラムを作れ。 (解答は付録 F.2 を見よ)。

問4      現象数理学科 Mac では、 cglsc コマンドでコンパイルすることによって、 GLSC というグラフィックス・ライブラリィが利用できる。 それを用いて解曲線を描く C プログラムを作れ。


複数の初期値についての解を表示してみよう。 gnuplot に数値データを読ませているとき、 e 1文字からなる行と空行 (行の最初で改行する) があると、 データの区切りになり、その後は別のプロットをすることになる。

図: $ \D x/\D t=x$, $ x(0)=1$ に対する Euler法の解曲線 ($ x(0)$を$ 0.2$刻みで変化)
Image euler1ex1multi

euler1ex1multi.c

/*
 * euler1ex1multi.c (Euler method for Malthusian model)
 *   複数の初期条件に対応する解を計算する
 */

#include <stdio.h>

double k = 1.0;

int main(void)
{
  int i, N;
  double t, x, dt;
  double f(double, double), x0;
  double Tmax;
  // 最終時刻
  Tmax = 1.0;
  // 時間刻み
  printf("# N: "); scanf("%d", &N);
  dt = Tmax / N;
  // 初期値 x(0)=x0 を 0.2 から 2.0 まで 0.2 刻みで変更する
  for (x0 = 0.2; x0 <= 2.0; x0 += 0.2) {
    // x(0)=x0 を初期条件とする。
    t = 0.0;
    x = x0;
    printf("# t x\n");
    printf("%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);
    }
    printf("e\n\n"); // e 1文字の行と空行
  }
  return 0;
}

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

next up previous contents
Next: 2.3.2 Runge-Kutta法のプログラム例 Up: 2.3 C言語のプログラム例 Previous: 2.3 C言語のプログラム例
桂田 祐史