第10回(2011年12月21日)


今回は冬休み用のレポート(紙面にて提出)がありますのでご確認ください. 次回は1月11日です.


今日の演習

今日は以下の内容で演習を行います.


今日の目標

gnuplot で計算結果を可視化する.


(前回の復習) gnuplot 入門

前回に引続き、gnuplot を用いて、関数やデータの可視化を行います.


まず前回の復習として、幾つかグラフを描画します.

ターミナルで,

gnuplot

と打ち込み,改行キーを押すことで gnuplot を起動することができます.起動すると,ターミナル上のプロンプトは gnuplot のものとなり,gnuplot のコマンドを入力できる状態となります.そこで例えば,

plot [-3:3][-5:5] x-2*x**2-1.5*x**3+x**4

とすると、横軸 -3 〜 3、縦軸 -5 〜 5 の範囲で x + 2x2 - 1.5x3 + x4 のグラフが描かれます.(xn は x**n と表記します.)

他に sin(x) や exp(0.5*x) といった関数も同様に描画できます.

またグラフの横軸(x)、縦軸(y)の表示範囲は

plot [(min_x):(max_x)][(min_y):(max_y)] (関数)

として、min_x、 max_x、min_y、max_y を与えることで決ります。指定しないと自動的に与えられます。


gnuplot は,

quit

で終了できます.


(前回の復習)y = sin(x) のグラフの描画

これから、gnuplot を使って y = sin(x) のグラフを描きます.もちろん,gnuplot では,

plot sin(x)

とすれば y = sin(x) のグラフが得られるのですが,gnuplot がどのようにしてそれを描いているかわかりません.そこで少々遠回りをします.

実際には y = sin(x) のグラフは滑らかな曲線ですが、コンピュータではそのような滑らかな曲線を直接は扱えませんので折れ線で近似して描画します.つまり、変数 x は、実際には連続的に変化するものですが、とびとびの値を使います.後の例で見るように、とびとびとはいえ、十分に細かく取れば、折れ線で曲線をうまく近似することができます.

図のように [0, L] 区間を N-1 等分して、その等分した小区間の幅を dx とします.N-1 等分すると図のように N 個の区切りが現れますが、それらに 0,1, ...., i, ...., N-1 と番号をつけることにします.また、 N を分割数、 dx を分割幅と呼びます.このとき、

Xi = i*dx

の部分についてそれぞれ y = sin(Xi) を計算し、それらを折れ線で結ぶことを考えます.C言語風に書けば

X[i] = dx*i;
Y[i] = sin(X[i]);

をそれぞれ求め、(X[0], Y[0]), (X[1],Y[1]), ...., (X[N-1], Y[N-1]) を結ぶ折れ線を描くことになります.これをプログラムにすると次のようになります.(前回のファイル 1214-1.c で N = 50 としたものです. 復習:#define 文で定数が定義され、また数学関数を使うので#include が加わっています.)


- - - - - - - - - - - - - - - -

#include <stdio.h>
/* 次の文は数学関数を使う場合に必要 */
#include <math.h>

/*** 定数の定義 ***/
#define N (50)
#define PI (3.1415926)
#define L (2*PI)

main()
{
  int i;
  float X[N], Y[N], dx;

  /* 分割幅 dx を求める */
  dx = L/(N - 1);

  /* X[i] を求める */
  for(i = 0; i < N; i++)
  {
    X[i] = i*dx;
  }

  /* Y[i] を求める */
  for(i = 0; i < N; i++)
  {
    Y[i] = sin(X[i]);
  }

  /* 結果を数値として表示 */
  for(i = 0; i < N; i++)
  {
    printf("%f %f\n", X[i], Y[i]);
  }
}

- - - - - - - - - - - - - - - -

データファイルの作成

上記プログラムを使って、近似的な sin(x) のグラフを見てみましょう.まず,プログラムをコンパイルします.プログラム内で数学関数(ここでは sin(x) )を使っているので、最後に -lm を付けてください.

cc 1214-1.c -o 1214-1 -lm

続いて実行してみましょう.

./1214-1

とすると、数値の列が表示されます.これら数値をgnuplotでグラフ化しましょう.この数値の列からなる、データファイル作成するには、次のようにします.

./1214-1 > N=50.data

上記のようにすることで、数値列は画面に表示されず、N=50.data というファイル名のファイルが作成され、そこに書き込まれます.

(ターミナル上で、more N=50.data などとすると、確認できます. ちなみに more はファイルの中身を見るコマンドで、more (ファイル名) というふうに使います.)
なお今回は、N=50なのでこのような名前にしましたが、ファイル名は基本的に自由です.

gnuplot を起動し,

plot [0:2*pi] "N=50.data" with linespoints, sin(x)

緑色で書かれたグラフが描きたいグラフ (y = sin(x) )です.赤色が描こうとしている折れ線です.(復習:カンマで区切ることで、複数の関数やファイルを同時に表示できます.)N=50 程度とれば(赤線)、見た目はほぼ y=sin(x) のグラフに見えます.(前回みたように、 N を小くすると、近似が悪くなって、形がかけはなれていきます。)


(前回の課題)円の描写

1214-1.cを以下のように少し改造すると、円を描くプログラムができます. 赤字の部分が改造個所です. ここでは Y[i], Z[i] を2次元座標とし、X[i] を媒介変数としています(N = 50). サンプルを参考に円や楕円を描いてみてください. sin や cos の係数を変えると楕円が描けます.

- - - - - - - - - - - - - - - -

#include <stdio.h>
/* 次の文は数学関数を使う場合に必要 */
#include <math.h>

/*** 定数の定義 ***/
#define N (50)
#define PI (3.1415926)
#define L (2*PI)

main()
{
  int i;
  float X[N], Y[N], Z[N], dx;

  /* 分割幅 dx を求める */
  dx = L/(N - 1);

  /* X[i] を求める */
  for(i = 0; i < N; i++)
  {
    X[i] = i*dx;
  }

  /* Y[i], Z[i] を求める */
  for(i = 0; i < N; i++)
  {
    Y[i] = sin(X[i]);
    Z[i] = cos(X[i]);
  }

  /* 結果を数値として表示 */
  for(i = 0; i < N; i++)
  {
    printf("%f %f %f\n", X[i], Y[i], Z[i]);
  }
}

- - - - - - - - - - - - - - - -

このプログラムを実行 して データファイルを作成し、gnuplot 内で

gnuplot> plot [-2:2][-2:2]'(データファイル名)' using 2:3 with lines
とすると、半径1の円が描かれます。

ちなみに

gnuplot> set size square
gnuplot> plot [-2:2][-2:2]'(データファイル名)' using 2:3 with lines

とすると、縦軸、横軸の縮尺が同じグラフが描けます.(set size squre で gnuplot の画面が正方形になります.)


先のプログラム例では、 説明の為,少々冗長なプログラムとなっていました.実際には以下のプログラムでも同様の結果が得られます.各人確認してください.

- - - - - - - - - - - - - - - -

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

/* 定数の定義 */
#define N (50)
#define PI (3.1415926)
#define L (2*PI)

main()
{
  int i;
  float dx;

  /* dx を求める */
  dx = L/(N - 1);

  /* 結果を画面に表示 */
  for(i = 0; i < N; i++)
  {
    printf("%f %f %f\n", i*dx, sin(i*dx), cos(i*dx));
  }
}

- - - - - - - - - - - - - - - -

課題1

円を描くプログラムを少し改造して( N を変化させて)、正三角形、正方形、正五角形、正六角形を描いてください.(N を適度に変化させて)


空白行のあるデータファイルの作成

前回見てきたように、データファイル中に空白行があると、その前後のデータを別々に表示するので、複数のグラフや図形を描写できます. Cプログラムでは、printf("\n")(何も書かずに改行のみする)を使うと、空白行のあるデータを作成できます。

例として、次のような空白行のあるデータを作るプログラムのサンプルを示しておきます.

0 0
0 1
0 2

1 0
1 1
1 2

2 0
2 1
2 2

(サンプルプログラム)

- - - - - - - - - - - - - - - -

#include <stdio.h>

main()
{

  int i,j;

  /** 繰り返し ( i ) **/
  for(i=0; i<3; i++)
    {

    /** 繰り返し ( j ) **/
     for(j=0; j<3; j++)
       {

       /** i, j ( = 0, 1, 2 ) を表示 **/
       printf("%d %d\n", i, j);

     }

  /** j についての繰り返し処理後に何も書かずに改行(空白行をつくる) **/
  printf("\n");

  }

}

- - - - - - - - - - - - - - - -

これをコンパイルして実行すると、上記のように出力されます. データファイルを作成して gnuplot で描かせると( plot [-1:3][-1:3]'(データファイル名)' with lines )、以下のように3本の直線が現れます。

以下の課題2では複数の円を描いてもらいますが、この空白行の生成法を応用してください。


課題2

以下の様な、6個の独立した(線等で繋がっていない)直径1(半径 0.5)の円が、正六角形状に並んでいる絵を描くプログラムを作成し、gnuplotで描画してください。円一つにつき50個程度の点(N=50)を使い、例えばデータファイルを en6-data とし

plot 'en6-data' w l

としたときに、以下と同様の絵が描かれること。(円同士は重なっていても離れていても構いません. プログラム中では、半径は 1/2 ではなく 0.5、 もしくは 1.0/2.0 として下さい.)

(ヒント)課題1で行ったように六角形の頂点を決めて、その各頂点を中心に円を一つづつ(繰り返し文を使って)描いていく。つまり課題1では六角形の各点を一つづつ順々(繰り返し文を使って)に出力(printf())していたのだけど, その代わりに一つ円を(繰り返し文を使って)描くようにすればよい(2重繰り返しになりますね.). ただし円を一つ描く度に改行コマンド printf("\n") を使って、データファイルに空白行をつくる(上のサンプルを参考に.).


数列の収束の様子をみる.その1

課題3

次の数列 Xn = ( 1 + 1/n )n ( 1 + 1/n の n 乗) を、n : 1 〜 200 程度まで求めるプログラムを作成し、n と Xn の関係を gnuplot で可視化せよ. ( n → ∞ で Xn → e (ネイピア数)ですね。収束の様子を確認してください。)

課題4

次の等式 π/4 = 1 - 1/3 + 1/5 - 1/7 + ... + (-1)n-1/(2n-1) + ... (ライプニッツ-グレゴリーの等式)を利用して、π の近似値を求めるプログラムを作成せよ。右辺の項の数とπの近似値 との関係を gnuplot で可視化せよ. ( n : 1 〜 1000 程度までで、おおよその傾向を見てください。)


数列の収束の様子をみる.その2(漸化式)

課題5

次の漸化式で表される数列 Xn を求めるプログラムを作成し、n (0 ~ 100)と Xn の関係を gnuplot で可視化せよ。

Xn+1 = 1 + 1/Xn

X0 = 1

( n → ∞ で Xn はいわゆる黄金比に収束します。収束の様子を確認してください。)

課題6

次の漸化式で表される数列 (0 < Xn < 1) を求めるプログラムを作成し、n と Xn の関係を gnuplot で可視化せよ.

Xn+1 = A * Xn * (1.0 - Xn)

A < 1 なら Xn → ∞ → 0 となるのは、すぐ分ると思います. では A > 1 なら...

A = 0.5, 2.0, 3.5, 3.7, 4.0 等(A : 0 〜 4)それぞれの場合について、初期値 X0 を 0.3, 0.6, 0.61 等(0 〜 1)と変化させ、実行、描画して見比べてください. X0, X1, X2 あたりの様子から X9, X49 あたりの値が予測可能なら予測せよ.(これはロジスティックマップとよばれるものです. 興味がありましたら何者なのか調べてみてください. )


課題7

3次元グラフを用いて2変数関数を表示してみましょう.

z=f(x,y), f(x,y) = sin(x)cos(y) を表示させるには gnuplot を起動して

splot [0:2*pi][0:2*pi][-2:2] sin(x)*cos(y) using 1:2:3 with linespoints

とすれば表示できます.

またこのサンプルは前回の y=sin(x) を描くプログラムを改造して、z = sin(x)cos(y) をgnuplot で書くためのデータを出力するものです. 改造していろんな関数を

splot [0:2*pi][0:2*pi][-2:2] "ファイル名" using 1:2:3 with linespoints

もしくは

splot [0:2*pi][0:2*pi][-2:2] "ファイル名" using 1:2:3 with pm3d

等として表示してみてください.

レポート問題(冬休み)

数値計算の参考書等で「オイラー法」について調べ、A4用紙1〜2枚程度にまとめよ。

レポート提出期限:2012年1月11日(水)当演習開始時に回収 

注意: 

レポートの始めに学籍番号と名前を書く事(表紙はいりません.).
レポートは「手書き」で書く事。丁寧に書いてください。
参考にした文献の出典(著者、書名、出版社、出版年)も明記すること。 


今日学んだ事

  • gnuplot とCプログラムでグラフ、図形を描画しました。

戻る