Part.8

定数・数学関数・二分法・gnuplot

内容

  • 定数の使い方
  • 数学関数の使い方
  • 2分法
  • gnuplotの簡単な使い方

目標

2分法のプログラムを理解し，gnuplot の使い方の基本を身につける．

定数（#define 文の使い方）

#define 文を使うと定数を定義することができます．通常、#include 文の後に書きます．例えば、

#define  N  (5)

と書いた場合、後のプログラム中の N は全て (5) に置き換わります．例えば、次のように使えます．

#include <stdio.h>

#define A (5)
#define PI (3.1415926)

int main()
{
    printf("%d, %f\n", A, PI);
}

#define 文は、実は単に文字列を置き換えるだけです．つまり、printf 関数部分に A および PI がありますが、それらはコンパイル時に次のように書き換えられます．

printf("%d, %f\n", (5), (3.1415926));

#define 文を上手に使うと、プラグラムが大変見通しのよいものとなります．

（参考）#define 文を使う場合私はカッコを必ず付けますが、それには理由があります。まず、通常の定数として使う場合には、問題ありません。例えば、先の例は、

#define N 5

と書いても同じ結果になります。しかしながら、例えば式を書く場合、

#define N1 (5+2)

と

#define N2 5+2

では、使い方に寄っては結果が異なる場合があります。

#include <stdio.h>

#define N1 (5+2)
#define N2 5+2

int main()
{
  printf("%d, %d\n", N1*2, N2*2);
}

実行すると、

14, 9

と表示されます。先の #define 文は文字列の置き換えということを思い出せば、

printf("%d, %d\n", (5+2)*2, 5+2*2);

となるので違いがわかるでしょう。カッコを付けておけば、こういうミス（バグ）が減ると考えているので、必ずカッコを付けるようにしています。

数学関数（三角関数等）の使い方

C言語では三角関数等の数学関数を用いることができます．数学関数を用いる場合には、

#include <math.h> をプログラムの先頭付近に書き入れることと、コンパイル時に、ファイル名を sample0.c とすると、

cc sample0.c -o sample0 -lm

というふうに、最後に -lm を付け足す必要があります（付けなくても良い場合もある）．三角関数を用いた例を次に示します．

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

#define PI (3.1415926)
#define N (10)

int main()
{
    int i;
    float dtheta;

    dtheta = 2*PI/(N-1);

    for(i = 0; i < N; i++)
    {
        printf("%f %f¥n", i*dtheta, sin(i*dtheta));
    }
}

#define 文の使い方にも注目してください．このように定数としてNを定義しておけば、Nの定義を変更するだけですみ、他の部分の変更は必要ありません．

2分法

講義で2分法やニュートン法といった非線形方程式の根を求めるアルゴリズムを習ったことと思います．ここでは2分法によって根を求めるプログラムを実際に作成し，アルゴリズムとプログラムとの対応をみましょう．

次の問題を考えます．次の関数 f(x) を考える．

f(x) = x2 - 2

f(x) = 0 となる正の実数 x を求めよ．もちろん答えは x = √2 ですが，それを2分法を用いて近似的に求めましょう．

gnuplot 入門

二分法では求めたい根を挟み込み，徐々に幅をせばめていくアルゴリズムですが，根を挟み込む値はこちらが前もって与える必要があります．今回の例については簡単にグラフを書くことができますが，そうでないこともありますし，関数のグラフを簡単に確かめたいこともあるでしょう．そういう場合に便利なのが gnuplot と呼ばれるアプリケーションです．ここでは簡単に gnuplot を用いたグラフの描画を実践してみましょう．

ターミナルで，

gnuplot

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

plot x*x - 2

とすることにより，グラフを得ることができます．x の範囲は自動的に決められますが，こちらで設定することもできます．x2は x**2 と表記します．

plot [0:2] x**2 - 2

同時に f(x) = 0 のグラフも描くとわかりやすいので，書き加えます．

plot [0:2] x**2 - 2, 0

以下のような結果が得られ，求めたい根は 1 と 2 の間にあることがわかります．

gnuplot は，

quit

で終了できます．

たとえば，次のようなものを試してみてください．

plot sin(x)
plot x**3 - x**2 + 5

以下に2分法のアルゴリズムを示します．なお，:= は代入を示します（C言語では = に対応）．

2分法のアルゴリズム

f(a) < 0, f(b) > 0 を満たすように変数 a, b の値を設定する
正定数 e を設定する
|a - b|/2 > e の間以下を繰り返す
{
  c := (a+b)/2
  fc := f(c)
  もし fc > 0 ならば b := c
  もし fc < 0 ならば a := c
}
c を答えとする

C言語による実装例

上記アルゴリズムをC言語のプログラムとして実装した例を示します．ここでは，fabs 関数，exit 関数と，break 文を使っています．

fabs 関数は引数として実数値をとり，戻り値としてその絶対値を返します．先に作成した myabs 関数と同様です．

exit 関数は引数として整数値をとり（与える整数値には意味はあるが，当面何でも良い），プログラムを終了させます．ここでは，A, B の値が不適切であれば，メッセージを表示して終了するようにしてあります．

break 文は，繰り返し文から強制的に抜ける場合に用います．ここでは fc == 0.0 という起こりそうに無い状況が起きた場合に使っています．

なお，このプログラムでは数学関数である fabs 関数を使っているので，コンパイル方法がこれまでと少し異なります．以下のプログラムをたとえば Part8.c として打ち込んだ場合，コンパイルは以下のようにします．

cc Part8.c -o Part8 -lm

最後の -lm が変更点で（マイナス・小文字のエル・小文字のエム），数学関数は標準関数では無いため，このような指定が必要です．不正確な説明ですが，数学関数を使うので，その数学関数集として定義されたものを自前のプログラムにひっつけるということです．たとえば sin 関数等も数学関数として使うことができますが，その場合にも同様に -lm が必要になります．

以下，プログラムソース．

/* 二分法 */
#include <stdio.h>
/* exit 関数を用いる為に必要 */
#include <stdlib.h>
/* fabs 関数を用いる為に必要 */
#include <math.h>

/* 関数 f(x) のプロトタイプ宣言 */
double f(double x);

/* f(A) < 0 となるような A を設定 */
#define A (1.0)
/* f(B) > 0 となるような B を設定 */
#define B (2.0)
/* 収束判定に用いる */
#define E (1.0e-10)

int main()
{
  /* 使用する変数の宣言 */
 double a, b, c, e, fc;
 
 /* 変数の初期化 */
 a = A;
 b = B;
 e = E;
 
 /* A, B の値の妥当性チェック */
 if(!((f(a) < 0) && (f(b) > 0)))
 {
  printf("A, B の値が不適切です．\n");
  /* プログラムを終了させる標準関数 */
  exit(0);
 }
 
 /* 二分法アルゴリズム */
 /* |a-b|/2 > e の間繰り返す */
 while(fabs(a - b)/2.0 > e)
 {
  c = (a + b)/2.0;

  fc = f(c);
  
  /* c の符号により a または b に c を代入 */
  if(fc > 0.0)
  {
   b = c;
  }
  else if(fc < 0.0)
  {
   a = c;
  }
  /* fc == 0.0 であれば，c は根であるから繰り返し文を強制的に抜ける */
  else
  {
   /* break 文によって強制的に繰り返し文を抜けることができる */
   break;
  }
 }
 
 /* c の値を結果として表示 */
 /* %10.8f とは，全部で10文字，小数点以下8桁で表示するという意味 */
 printf("答 = %10.8f\n", c);
}

/* 関数 f(x) の定義 */
double f(double x)
{
 return (x*x - 2.0);
}

なお，上記のサンプルプログラムを少々変更したプログラムを次に載せます。 sample2.c というファイル名で保存して下さい．

/* 二分法  sample2.c */
#include <stdio.h>
/* exit 関数を用いる為に必要 */
#include <stdlib.h>
/* fabs 関数を用いる為に必要 */
#include <math.h>

/* 関数 f(x) のプロトタイプ宣言 */
double f(double x);

/* f(A) < 0 となるような A を設定 */
#define A (1.0)
/* f(B) > 0 となるような B を設定 */
#define B (2.0)
/* 収束判定に用いる */
#define E (1.0e-10)

int main()
{
  /* 使用する変数の宣言 */
 double a, b, c, e, fc;
 int i;
 
 /* 変数の初期化 */
 i = 0;
 a = A;
 b = B;
 e = E;
 
 /* A, B の値の妥当性チェック */
 if(!((f(a) < 0) && (f(b) > 0)))
 {
  printf("A, B の値が不適切です．\n");
  /* プログラムを終了させる標準関数 */
  exit(0);
 }
 
 /* 二分法アルゴリズム */
 /* |a-b|/2 > e の間繰り返す */
 while(fabs(a - b)/2.0 > e)
 {
  c = (a + b)/2.0;
  
  i++;
  /* c の値を結果として表示 */
  /* %10.8f とは，全部で10文字，小数点以下8桁で表示するという意味 */
  printf("%d %10.8f\n", i, c);

  fc = f(c);
  
  /* c の符号により a または b に c を代入 */
  if(fc > 0.0)
  {
   b = c;
  }
  else if(fc < 0.0)
  {
   a = c;
  }
  /* fc == 0.0 であれば，c は根であるから繰り返し文を強制的に抜ける */
  else
  {
   /* break 文によって強制的に繰り返し文を抜けることができる */
   break;
  }
 }
}

/* 関数 f(x) の定義 */
double f(double x)
{
 return (x*x - 2.0);
}

このプログラムは，繰り返しごとに繰り返し回数とそのときの c の値を表示するようにしてあります．このプログラムを実行すると，数値の列が表示されます．それを gnuplot でグラフ化してみましょう．

ターミナルで，

cc sample2.c -o sample2 -lm
./sample2

上記のようにすることで，数値の列がみられ，繰り返しごとに√2の値に近づく様子がわかります．この数値の列をグラフとしてみるには次のようにします．

./sample2 > data

このようにすることで，画面には表示されず，data というファイル名のファイルに書き込まれます．たとえばエディターで data というファイルを開いてみればわかります．さて，gnuplot で data 内の数値を可視化しましょう．

ターミナルで，

gnuplot

として gnuplot を起動し，

plot "data" with linespoints, sqrt(2)

とすると，次のようなグラフが得られます．

gnuplotの使い方２

2分法において，c の値が繰り返し毎にどのように変化するかをみました．他の a, b はどのように変化するでしょうか？すべてをまとめてグラフとしてみることはできないでしょうか？そのために，sample2.c を次のように少々変更します．変更部分は太字部分です．sample2.cを変更しsample3.cとして保存してください．

/* 二分法 sample3.c */
#include <stdio.h>
/* exit 関数を用いる為に必要 */
#include <stdlib.h>
/* fabs 関数を用いる為に必要 */
#include <math.h>

/* 関数 f(x) のプロトタイプ宣言 */
double f(double x);

/* f(A) < 0 となるような A を設定 */
#define A (1.0)
/* f(B) > 0 となるような B を設定 */
#define B (2.0)
/* 収束判定に用いる */
#define E (1.0e-10)

int main()
{
 /* 使用する変数の宣言 */
 double a, b, c, e, fc;
 int i;
 /* 変数の初期化 */
 i = 0;
 a = A;
 b = B;
 e = E;

 /* A, B の値の妥当性チェック */
 if(!((f(a) < 0) && (f(b) > 0)))
 {
  printf("A, B の値が不適切です．\n");
  /* プログラムを終了させる標準関数 */
  exit(0);
 }

 /* 二分法アルゴリズム */
 /* |a-b|/2 > e の間繰り返す */
 while(fabs(a - b)/2.0 > e)
 {
  c = (a + b)/2.0;

  i++;
  /* a,b,c の値を結果として表示 */
  /* %10.8f とは，全部で10文字，小数点以下8桁で表示するという意味 */
  printf("%d %10.8f %10.8f %10.8f\n", i, a, b, c);

  fc = f(c);

  /* c の符号により a または b に c を代入 */
  if(fc > 0.0)
  {
   b = c;
  }
  else if(fc < 0.0)
  {
   a = c;
  }
  /* fc == 0.0 であれば，c は根であるから繰り返し文を強制的に抜ける */
  else
  {
   /* break 文によって強制的に繰り返し文を抜けることができる */
   break;
  }
 }
}

/* 関数 f(x) の定義 */
double f(double x)
{
 return (x*x - 2.0);
}

ターミナルで，

cc sample3.c -o sample3 -lm
./sample3

上記のようにすることで，数値の列がみられます．sample2 では2列のデータでしたがsample3では4列のデータが得られます．プログラムを見ればわかりますが，1列目は繰り返し数，2列目はaの値，3列目はbの値，4列目はcの値となっています．先ほどと同様に，

./sample3 > data3

このようにすることで，画面には表示されず，data3 というファイル名のファイルに書き込まれます．たとえばエディターで data3 というファイルを開いてみればわかります．さて，gnuplot で data3 内の数値を可視化しましょう．

ターミナルで，

gnuplot

として gnuplot を起動し，

plot "data3" with linespoints, sqrt(2)

とすると，（期待に反して）先とほぼ同様の結果が得られます．gnuplotは何も指示しなければ，横軸に1列目を用い，縦軸として2列目を用います．同時に a, b, c の時間変化をみたい場合には次のようにする必要があります．gnuplotに続けて以下を打ち込んでみてください．

plot "data3" using 1:2 with linespoints, "data3" using 1:3 with linespoints, "data3" using 1:4 with linespoints

using によってどの列を横軸:縦軸に使うかを指示できます．上のグラフによって，c（青線)はa（赤線）とb（緑線）の中点となっていることがわかります．グラフの見栄えをもう少しよくしましょう．with linespoints は w lp と省略できます．

set title "Graph"
set xlabel "i"
set ylabel "a,b,c"
plot "data3" using 1:2 w lp t "a", "data3" using 1:3 w lp t "b", "data3" using 1:4 w lp t "c"

さて，上記のようなグラフを画像ファイルとして保存し，ホームページに貼り付けたりするようなことがあるかもしれません．そういう場合には，これまでに続いて次のようにします．

set term png
set output "graph.png"
replot
quit

これで graph.png という画像ファイルができているはずです．グラフの png ファイルを作りたい場合の手順をまとめると，gnuplot を起動し，次のようにします．

set term png
set output "sin_graph.png"
set title "Graph of sin(x)"
set xlabel "x"
set ylabel "sin(x)"
plot sin(x)
quit

title や plot の部分は適宜変更します．png 部分（2カ所）を jpg に変更すれば jpeg 形式で保存されます．詳しくは，gnuplot を起動して，

set term

と打ち込むと出力ファイル形式の一覧が得られますから参考にしてください．

課題１

2分法において，点(a,b) が時間とともにどのように変化するかを gnuplot を用いてみてみよう．次のコマンドにてグラフを描いてみてください．2番目の例は，横軸および縦軸の範囲を指定する方法を示すためのものです．点(a,b)が点(√2,√2)に近づく様子が見られるはずである．

plot "data3" using 2:3 w lp
plot [1.4:1.42] [1.35:1.45] "data3" using 2:3 w lp

課題２

明日から毎日体重を量り，ノートに記録開始日からの経過日数および体重を書いておく．その記録を冬休み明けにテキストエディターに打ち込み，gnuplot を用いて体重の変化をグラフとしてみてみよう（もちろん提出の必要は無い）．（エクセル等を使えばよりスマートにできますが，ここでは gnuplot の使い方の練習なのでgnuplot を使ってみてください）

課題３

2分法ではなくニュートン法を用いたプログラムを作成し，2分法に比べてニュートン法の方が収束がはやいことを確かめてください．（実用的にはニュートン法を使います．ニュートン法は必ず調べましょう．）

グラフを印刷したい場合もあるかもしれません．その場合には，ps ファイルに保存すると便利です．その場合の手順は次のようになります．(後々論文に図を貼り込むこともあるかもしれないので，epsという形式で保存することにしましょう)

set term postscript eps
set output "sin_graph.eps"
set title "Graph of sin(x)"
set xlabel "x"
set ylabel "sin(x)"
plot sin(x)
quit

これによって，"sin_graph.eps"というファイルができあがります．そのファイルをファイルブラウザでダブルクリックすると適当なpsファイルビューアが起動するでしょうからそれから印刷してください．(ps ファイルとは PostScriptファイルの略です．PostScriptとはプリンタの制御言語で，プリンタに絵を描かせる命令の集まりです．)

Part.8で学んだ事

  • 2分法
  • gnuplot を使いました