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 % |
gnuplot については、ネットに解説があふれている (例えば桂田 [7])。 この問題について役に立ちそうなことをいくつか説明しておく。
gnuplot の細かな工夫に関する注意
test1.gp |
plot "euler1ex1.data" with lp set term png set output "euler1ex1.png" replot |
% gnuplot test1.gp |
問1 上の実行例で最後に出力した数値 2.704814 の誤差はいくらか。 を変えることで誤差がどのように変わるか調べよ (両側対数目盛でプロットすることを勧める)。
問2 上の実行例では、出力をリダイレクトすることでデータファイルを作っているが、 fopen(), fclose(), fprintf() を使ってデータファイルを作れ (参考 「簡単なファイル入出力」, 解答は付録 F.1 を見よ)。
問3 Cのプログラムの中で popen() という関数を使うと、 外部プログラムを起動して通信ができる。 gnuplot を起動して解曲線を描く C プログラムを作れ。 (解答は付録 F.2 を見よ)。
問4 現象数理学科 Mac では、 cglsc コマンドでコンパイルすることによって、 GLSC というグラフィックス・ライブラリィが利用できる。 それを用いて解曲線を描く C プログラムを作れ。
複数の初期値についての解を表示してみよう。 gnuplot に数値データを読ませているとき、 e 1文字からなる行と空行 (行の最初で改行する) があると、 データの区切りになり、その後は別のプロットをすることになる。
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; } |
桂田 祐史