第15章 H25.2.27: 数値的に積分する

積分領域を微小領域に区切り, 各領域で被積分関数$f(x)$を多項式で近似することで, 数値的に積分値を求めることができる.

• 0次式: 領域 $[x_i,x_{i+1}]$ において被積分関数$f(x)$を $f(x_i)$で近似する.

• 1次式 (台形公式): 領域 $[x_i,x_{i+1}]$ において被積分関数$f(x)$を $(x_i, f(x_i))$, $(x_{i+1}, f(x_{i+1}))$を通る1次式(直線)で近似する.

• 2次式 (シンプソンの公式): 領域 $[x_{2i},x_{2i+2}]$ において被積分関数$f(x)$を $(x_{2i}, f(x_{2i}))$, $(x_{2i+1}, f(x_{2i+1}))$, $(x_{2i+2}, f(x_{2i+2}))$ を通る2次式で近似する.

15.1 本日の課題

• 台形公式による数値積分のプログラム[15.3]を理解し, 自分で書けるようにすること.

• 数値積分の方法に台形公式とシンプソンの公式がある. この2つの方法を用いて, 各自で被積分関数と積分領域を与え, 刻み幅と誤差の関係を調べる. また, 2つの方法の差を調べる.

• fscanf を用いてファイルからの入力するプログラム[15.5.1](特にfopenの行)と, ファイルの最後まで値を入力するプログラム[15.5.3](特にfscanfの行)を理解し, 自分で書けるようにする.

最後に, 以下の自分で作成, 実行した結果をメールで送る.

• ファイルから入力するプログラム, それを実行するために必要なデータのファイル.

• 台形公式とシンプソンの公式で刻み幅と誤差の関係を調べるプログラム. 調べたときの画面への出力, SVGファイル.

svg1.h は自分で変更した部分が無い場合は, 送らないこと.

【注意】

• ここで挙げるプログラムのうち描画を行うものは先頭で #include "svg1.h"と記述されている. これをコンパイルするためには, プログラムと同じフォルダ/ディレクトリに 予め svg1.h[14.3] を作成しておく必要がある.

• gcc で数学関数を用いるプログラムをコンパイル場合は -lm オプションをつける必要がある.

15.2 0次式 (長方形) で近似

定積分を数値的に求める. そこでまず, 積分区間を分割し, 各区間の面積を長方形の面積で近似する. これは積分区間 $[x_i,x_{i+1}]$の関数$f(x)$の値を, 一定の値$f(x_i)$で近似することを意味する. すなわち各区間内の関数$f(x)$の値は0次式で近似されていることになる.

15.2.1 まずは 10 区間に分割して足し合わせてみる

定積分

$$\displaystyle \int_{\mathrm{x\_init}}^{\mathrm{x\_final}} f(x) \mathrm{d}x$$ (38)

を数値的に求めることを考える.

Figure 6: 数値的な積分.

図6に示すように積分区間 [x_init, x_final] を $n$ 分割する. $i$ 番目の区間で決まる部分の面積 $s_i$ は 近似的に $s_i \simeq f(x_i) \times h$ で表される.

具体的な例として定積分

$$\displaystyle \int_0^{1} \cos(x) \mathrm{d}x$$ (39)

を求めてみる. まず全体を 10 分割し全体の面積を求めてみる．

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

double f(double x)
{
  double fx;

  fx = cos(x);
  return fx;
}

int main(void)
{
  int n;
  double x_init, x_final, h;
  double x0, x1, x2, x3, x4, x5, x6, x7, x8, x9;
  double s;

  x_init = 0.0;
  x_final = 1.0;
  n = 10;

  h = (x_final - x_init) / n;
  x0 = x_init;
  s = 0.0;

  s = s + f(x0) * h;

  x1 = x0 + h;
  s = s + f(x1) * h;

  x2 = x1 + h;
  s = s + f(x2) * h;

  x3 = x2 + h;
  s = s + f(x3) * h;

  x4 = x3 + h;
  s = s + f(x4) * h;

  x5 = x4 + h;
  s = s + f(x5) * h;

  x6 = x5 + h;
  s = s + f(x6) * h;

  x7 = x6 + h;
  s = s + f(x7) * h;

  x8 = x7 + h;
  s = s + f(x8) * h;

  x9 = x8 + h;
  s = s + f(x9) * h;

  printf("%g, %g\n", s, sin(1.0));
  
  return 0;
}

15.2.2 for を使って簡略化

for を使って下のように簡潔に表現できる．

分割数 $n$ を大きくすると刻み幅 $h$ が小さくなる. このとき誤差はどのようになるか調べよ.

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

double f(double x)
{
  double fx;

  fx = cos(x);
  return fx;
}

int main(void)
{
  int i, n;
  double x_init, x_final, x, h;
  double s;

  x_init = 0.0;
  x_final = 1.0;
  n = 10;

  h = (x_final - x_init) / n;
  x = x_init;
  s = 0;

  for(i=0; i<n; i++){
    s = s + f(x) * h;
    x = x + h;
  }

  printf("numerical value: %g, analytical value: %g\n", s, sin(1.0));
  
  return 0;
}

15.3 1次式(台形公式)による積分

区間 $[x_i,x_{i+1}]$における関数$f(x)$の値を, 点$(x_i, f(x_i))$ と 点 $(x_{i+1}, f(x_{i+1}))$ を通る直線で近似すると, より精度の高い近似となる.

このとき面積$S$は,

$$\begin{eqnarray*} S &\approx& S_0 + S_1 + \cdots + S_i + \cdots + S_{n-1}\\ &=&... ...)}{2} + \sum_{i=1}^{n-1} f(x_i) + \frac{f(x_n)}{2} \right\} h\\ \end{eqnarray*}$$

ただし各区間の長さ（すなわち台形の底辺の長さ）を $x_{i+1} - x_i = h (x=0, 1, ..., n-1)$ とした.

Figure 7: 台形法による積分.

例として,

$$\displaystyle \int_0^1 \frac{4}{1+x^2} \mathrm{d} x = [\arctan x]{}_0^1 = \pi$$ (40)

を数値的に積分してみると以下のようになる. このプログラムでは, func に積分関数を設定しておき, daikei 関数に積分範囲の初めと終わり, および区間の分割数を与えて呼び出している.

台形公式による数値積分のプログラム:

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

/* Eqn (5.7) in p. 112. */

double func(double x)
{
	 double ret;

	 ret = 4.0/(1.0+x*x);
	 return ret;
}

double daikei(double x0, double xn, int n)
{
	 int i;
	 double y0 = func(x0);
	 double yn = func(xn);
	 double h = (xn - x0) / n;
	 double sum = 0.0;
	 double yi, xi;
  
	 sum = (y0 + yn)/2.0;
  
	 for(i=1; i<n; i++){
		  xi = x0 + h*i;
		  yi = func(xi);
		  sum += yi;
	 }
	 return sum * h;
}

int main(void)
{
	 double I;

	 I = daikei(0.0, 1.0, 10);
	 printf("daikei: I=%g\n", I);
  
	 return 0;
}

分割数を増やすと誤差がどのように変化するか調べてみる. この場合, 解析解(analytical_value=$\pi$)が分かっているので, 数値解(numerical_value)の誤差は $\pi$ - numerical_value である. この値は何桁も変化するので, $\log_{10} (\texttt{analytical\_value} - \texttt{numerical\_value})$ の値をプロットすることにする. この値は何桁目まで正しいかを表している. たとえば n=4000のとき誤差は$10^{-8}$の程度である. ただし, 一般にはanalytical_value - numerical_value は負の値になる場合もあるので, エラーにならないように, 実数の絶対値を返すfabs関数により, log10関数の引数を常に正になるようにしている.

また, M_PIは 環境によっては math.hの中で有効桁の範囲で円周率$\pi$に定義されている. 定義されていない場合は次のように自分で定義しておく必要がする.

#define M_PI 3.14

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


double func(double x)
{
    double ret;

    ret = 4.0/(1.0+x*x);
    return ret;
}

double daikei(double x0, double xn, int n)
{
    int i;
    double y0 = func(x0);
    double yn = func(xn);
    double h = (xn - x0) / n;
    double sum = 0.0;
    double yi, xi;
  
    sum = (y0 + yn)/2.0;
  
    for(i=1; i<n; i++){
	xi = x0 + h*i;
	yi = func(xi);
	sum += yi;
    }
    return sum * h;
}

int main(void)
{
    int n; // 分割数
    
    // プロットに関係する変数.
    SVG svg_file;
    PLOT plot;
    PLOT_VIEW view = {100.0, 50.0, 800.0, 500.0};
    PLOT_WORLD world = {0.0, -10.0, 10000.0, 0.1};
    PLOT_AXIS axis = {"n", 5, "log10(error) = log10(M_PI - I)", 5};
    ATTRIBUTE circle_attrib   = {"blue", 1.0, "white"};
	 
    // 数値計算に関係する変数.
    double analytical_value; // 解析解 (数学的に得られた値)
    double numerical_value;  // 数値解 (数値計算により得られた値)
    double error_order;      // 数値解の誤差の桁 (何桁目まで正しいかを表す指標)

    // プロット初期化
    plot_file_open(&svg_file, "drawing.svg", 855.0, 605.0);
    plot_init(&plot, &svg_file, &view, &world, &axis);

    // 分割数を変えて積分を数値的に求める.
    for(n=10; n<10000; n+=100){
	numerical_value = daikei(0.0, 1.0, n);
	analytical_value = M_PI;
	error_order = log10( fabs( numerical_value - analytical_value) );
	// printf("daikei: I=%g\n", error_order);
	plot_circle_draw(&plot, n, error_order, 2.0, &circle_attrib);
    }

    // プロット終了
    plot_file_close(&plot);
    
    return 0;
}

実行結果:

15.4 2次式 (シンプソンの公式) による積分

区間 $[x_{2i},x_{2i+2}]$における関数$f(x)$の値を, 3点 $(x_{2i}, f(x_{2i}))$, $(x_{2i+1}, f(x_{2i+1}))$, $(x_{2i+2}, f(x_{2i+2}))$ を通る2次曲線$P_i(x)$で近似すると, より精度の高い近似となる. $P_i(x)$は2次多項式でこれを積分すると

$$\displaystyle I_i = \int_{x_{2i}}^{x_{2i+2}} P_i(x) \mathrm{d}x = \frac{h}{3} f(x_{2i}) + \frac{4h}{3} f(x_{2i+1}) + \frac{h}{3} f(x_{2i+2})$$ (41)

となる. ここでは, 各領域の間隔を $x_{2i+1} - x_{2i} = x_{2i+2} - x_{2i+1} = h$ とした. 計算は「シンプソンの公式の導出」[15.4.2]に示したのでそちらを参照のこと. 積分領域$[a, b]$を$2n$等分するとき 全領域の積分 $$\int _a^b f(x) \mathrm{d}x$$は

$$\int _a^b f(x) \mathrm{d}x = \int _{x_0}^{x_{2n}} f(x) \mathrm{d}x = \sum_{i=0}^{n-1} I_i$$ (42)

$$\simeq \displaystyle \left( \frac{h}{3} f_0 + \frac{4h}{3} f_{1} + \frac... ...left( \frac{h}{3} f_2 + \frac{4h}{3} f_{3} + \frac{h}{3} f_{4} \right) + \cdots$$ (43)

$$\displaystyle + \left( \frac{h}{3} f_{2i-2} + \frac{4h}{3} f_{2i-... ...ac{h}{3} f_{2i} + \frac{4h}{3} f_{2i+1} + \frac{h}{3} f_{2i+2} \right) + \cdots$$ (44)

$$\displaystyle + \left( \frac{h}{3} f_{2n-4} + \frac{4h}{3} f_{2n-... ...left( \frac{h}{3} f_{2n-2} + \frac{4h}{3} f_{2n-1} + \frac{h}{3} f_{2n} \right)$$ (45)

$$= \displaystyle \frac{h}{3} f_0 + \frac{4h}{3} f_{1} + \displaystyl... ...displaystyle \frac{2h}{3} f_{2n-2} + \frac{4h}{3} f_{2n-1} + \frac{h}{3} f_{2n}$$ (46)

$$= \displaystyle \frac{h}{3}f_{0} + \frac{4h}{3}\sum_{i=1}^n f_{2i-1} + \frac{2h}{3}\sum_{i=1}^{n-1} f_{2i} + \frac{h}{3}f_{2n}$$ (47)

となる.

例として,

$$\displaystyle \int_0^1 \frac{4}{1+x^2} \mathrm{d} x = [\arctan x]{}_0^1 = \pi$$ (48)

をシンプソンの公式を用いて数値的に積分してみると以下のようになる.

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

double func(double x)
{
  double ret;
  ret = 4.0/(1.0+x*x);
  return ret;
}

double Simpson(double x_0, double x_2n, int n)
{
	 double h = (x_2n - x_0) / (2.0*n);
	 double sum = 0.0;
	 int j;
	 double x_2j, x_2j1, y_2j, y_2j1;

	 sum  = h/3.0* func(x_0);

	 for(j=1; j<=n; j++){
		  x_2j1 = x_0 + (2*j-1)*h;
		  y_2j1 = func(x_2j1);
		  sum += h/3.0*4.0*y_2j1;
	 }

	 for(j=1; j<n; j++){
		  x_2j = x_0 + 2*j*h;
		  y_2j = func(x_2j);
		  sum += h/3.0*2.0*y_2j;
	 }

	 sum += h/3.0* func(x_2n);

	 return sum;
}

int main(void)
{
	 double I;
  
	 I = Simpson(0.0, 1.0, 5);
	 printf("Simpson: I=%g\n", I);
  
	 return 0;
}

15.4.1 台形公式とシンプソンの公式の比較

台形公式($I_\mathrm{D}$)とシンプソンの公式($I_\mathrm{S}$)を比較してみる.

$$I_\mathrm{D} = \displaystyle \frac{h}{2} f_0 + \left( h f_1 + h f_2 + \cdots + h f_i + \cdots h f_{2n-1} \right) + \frac{h}{2}f_{2n}$$ (49)

$$I_\mathrm{S} = \displaystyle \frac{h}{3} f_0 + \left( \frac{4h}{3} f_{1} + \disp... ...tyle \frac{2h}{3} f_{2n-2} + \frac{4h}{3} f_{2n-1} \right) + \frac{h}{3} f_{2n}$$ (50)

台形公式では, 一定の重みで加えているのに対し, シンプソンの公式では, 隣接する組を $$\frac{2}{3}:\frac{4}{3}$$の重みを付けて足し合わせている点が異なる.

以下のプログラムでは, 台形法とシンプソン法により積分 $$\int_0^1 \frac{4}{1+x^2}\mathrm{d}x$$の値を求め, 解析解との差を 対数をとってから表示している.

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

double func(double x)
{
	 return 4.0/(1.0+x*x);
}

double daikei(double x0, double xn, int n)
{
	 int i;
	 double h = (xn - x0) / n;
	 double sum = 0.0;
	 double xi;
  
	 sum = (func(x0) + func(xn)) / 2.0;
  
	 for(i=1; i<n; i++){
		  xi = x0 + h*i;
		  sum += func(xi);
	 }
	 return sum * h;
}

double Simpson(double x_0, double x_2n, int n)
{
	 double h = (x_2n - x_0) / (2.0*n);
	 double sum = 0.0;
	 int i;

	 sum  = func(x_0);
	 for(i=1; i<=n; i++)
		  sum += 4.0 * func( x_0 + (2*i-1) * h );
	 for(i=1; i<n; i++)
		  sum += 2.0* func(x_0 + 2*i * h);
	 sum += func(x_2n);

	 return sum * h / 3.0;
}

int main(void)
{
	 int n; // 分割数/2
	 
	 // プロットに関係する変数.
	 SVG svg_file;
	 PLOT plot;
	 PLOT_VIEW view = {100.0, 50.0, 800.0, 500.0};
	 PLOT_WORLD world = {0.0, -20.0, 120.0, 0.1};
	 PLOT_AXIS axis = {"n", 5, "log10 |error| = log10 | M_PI - I |", 5};
	 ATTRIBUTE daikei_attrib   = {"blue", 1.0, "white"};
	 ATTRIBUTE simpson_attrib   = {"red", 1.0, "white"};

	 // 数値計算に関係する変数.
	 double I_D, I_S;

	 // プロット初期化
	 plot_file_open(&svg_file, "drawing.svg", 855.0, 605.0);
	 plot_init(&plot, &svg_file, &view, &world, &axis);

	 // 分割数を変えて積分を数値的に求める.
	 for(n=10; n<50; n+=2){
		  I_D = daikei(0.0, 1.0, 2*n);
		  I_S = Simpson(0.0, 1.0, n);
		  plot_circle_draw(&plot, n, log10( fabs(M_PI - I_D) ), 2.0, &daikei_attrib);
		  plot_circle_draw(&plot, n, log10( fabs(M_PI - I_S) ), 2.0, &simpson_attrib);
	 }

	 // プロット終了
	 plot_file_close(&plot);
	 
	 return 0;
}

実行例:

15.4.2 シンプソンの公式の導出

シンプソンの公式は以下のようにして得られる.

関数$f(x)$を $[x_0, x_{2n}]$の区間で数値的に積分することを考える. この区間を$2n$等分に分割する. このうちの隣接した2つの領域 $[x_i, x_{i+2}]$における積分に着目する. このとき, 関数$f(x)$は, $(x_i, f(x_i)$, $x_{i+1}, f(x_{i+1})$, $(x_{i+2}, f(x_{i+2}))$ の3点を通る. シンプソン法では, これら3点を通る2次式で関数を近似し, この2次式を積分する ことで, この区間における積分値とする.

この方法に従って実際に積分を求めてみる. ここでは, 各区間の間隔を$h$とし, $x_{i} = p$, $x_{i+1} = m$, $x_{i+1} = q$と置くことにする. 3点 $(p, f(p))$, $(m, f(m))$, $(q, f(p))$を通る2次曲線を$P(x)$とすると, $P(x)$は

$$P(x) = \displaystyle \frac{(x-m)(x-q)}{(p-m)(p-q)} f(p) +\display... ...x-p)(x-q)}{ (m-p)(m-q) } f(m) +\displaystyle \frac{(x-p)(x-m)}{(q-p)(q-m)} f(q)$$ (51)

と一意に決まる. 実際3点を通る2次曲線は一意に決まるが, この$P(x)$の式は2次式であるし, $(x, y)$として 3点$(p, f(p))$, $(m, f(m))$, $(q, f(q))$をそれぞれ代入してみると, 3点とも$y=P(x)$が成立し, 曲線上にあることが分かる.

$m=(p+q)/2$であることに注意して, $P(x)$を$[p, q]$の領域で積分すると次のようになる.

$$\int_p^q P(x) \mathrm{d}x$$ (52)

$$= \displaystyle \int_p^q \frac{(x-m)(x-q)}{(p-m)(p-q)} f(p) + \frac{(x-p)(x-q)}{ (m-p)(m-q) } f(m) + \frac{(x-p)(x-m)}{(q-p)(q-m)} f(q) \mathrm{d}x$$ (53)

$$= \displaystyle \frac{(q-p)(q+2p-3m)}{ 6 (p-m)} f(p) - \frac{(q-p)^2}{(p-m)(q-m)} f(m) + \frac{(q-p)(2 q+p-3m)}{6 (q-m)} f(q)$$ (54)

$$= \displaystyle \frac{q-p}{6}f(p) + \frac{2(q-p)}{3} f(m) + \frac{q-p}{6} f(q)$$ (55)

$$= \displaystyle \frac{h}{3} f(p) + \frac{4h}{3} f(m) + \frac{h}{3} f(q) \hspace{2cm} \because q-p =2h$$ (56)

$$= \displaystyle \frac{h}{3} f(x_i) + \frac{4h}{3} f(x_{i+1}) + \frac{h}{3} f(x_{i+2})$$ (57)

この面積は$2h$分の領域の積分を与えていることに注意すること.

さらに, 次の図のように全積分領域$[x_0, x_{2n}]$を$2n$等分し積分を考える. ここで各区間の幅は$h$である. 積分 $S = \displaystyle \int_a^b f(x) dx$ は, $f(x_i) = f_i$と表記することにすると,

$$S = \int_{x_{0}}^{x_{1}} f(x) \mathrm{d}x + \int_{x_{1}}^{x_{2}} f(x)... ..._{4}}^{x_{5}} f(x) \mathrm{d}x + \int_{x_{5}}^{x_{6}} f(x) \mathrm{d}x + \cdots$$ (58)

$$+ \int_{x_{2i}}^{x_{2i+1}} f(x) \mathrm{d}x + \int_{x_{2i+1}}^{x_... ...{2i+5}} f(x) \mathrm{d}x + \int_{x_{2i+5}}^{x_{2i+6}} f(x) \mathrm{d}x + \cdots$$ (59)

$$+ \int_{x_{2n-6}}^{x_{2n-5}} f(x) \mathrm{d}x + \int_{x_{2n-5}}^{... ...{2n-2}}^{x_{2n-1}} f(x) \mathrm{d}x + \int_{x_{2n-1}}^{x_{2n}} f(x) \mathrm{d}x$$ (60)

$$= \displaystyle \left( \frac{h}{3} f_0 + \frac{4h}{3} f_{1} + \frac... ...ft( \frac{h}{4} f_{2} + \frac{4h}{3} f_{5} + \frac{h}{3} f_{6} \right) + \cdots$$ (61)

$$+ \displaystyle \left( \frac{h}{3} f_{2i-2} + \frac{4h}{3} f_{2i-... ...{h}{3} f_{2i+2} + \frac{4h}{3} f_{2i+3} + \frac{h}{3} f_{2i+4} \right) + \cdots$$ (62)

$$+ \displaystyle \left( \frac{h}{3} f_{2n-6} + \frac{4h}{3} f_{2n-... ...left( \frac{h}{3} f_{2n-2} + \frac{4h}{3} f_{2n-1} + \frac{h}{3} f_{2n} \right)$$ (63)

$$= \displaystyle \frac{h}{3} f_0 + \frac{4h}{3} f_{1} + \displaystyl... ...c{4h}{3} f_{3} + \displaystyle \frac{2h}{4} f_{2} + \frac{4h}{3} f_{5} + \cdots$$ (64)

$$+ \frac{4h}{3} f_{2i-1} + \displaystyle \frac{2h}{3} f_{2i} + \fr... ...f_{2i+1} + \displaystyle \frac{2h}{3} f_{2i+2} + \frac{4h}{3} f_{2i+3} + \cdots$$ (65)

$$+ \displaystyle \frac{4h}{3} f_{2n-5} + \displaystyle \frac{2h}{3... ...displaystyle \frac{2h}{3} f_{2n-2} + \frac{4h}{3} f_{2n-1} + \frac{h}{3} f_{2n}$$ (66)

$$= \displaystyle \frac{h}{3} f_0 + \frac{4h}{3} f_{1} + \displaystyl... ...displaystyle \frac{2h}{3} f_{2n-2} + \frac{4h}{3} f_{2n-1} + \frac{h}{3} f_{2n}$$ (67)

$$= \displaystyle \frac{h}{3}f_{0} + \frac{4h}{3}\sum_{i=1}^n f_{2i-1} + \frac{2h}{3}\sum_{i=1}^{n-1} f_{2i} + \frac{h}{3}f_{2n}$$ (68)

15.5 ファイル操作: fscanf

ここまででは, fprintf 関数によりファイルに出力することを行ってきた. ここでは, fscanf関数によりファイルから入力する.

fscanf関数 は以下の点で scanf 関数と異なるがそれ以外は同様に使用できる.

• fopen 関数により予めファイルをオープンしてファイルポインタ FILE *fp を得る.
• fscanf を呼び出すときは常に出力先のファイルとして fp を第一引数に指定する. 第2引数以降はscanfと同じ.

• fclose 関数により, ファイルへの書き込みが終了したらファイルを閉じる.

fscanf 関数

#include <stdio.h>
int fscanf(FILE *fp, const char *format, ...);
int  scanf(          const char *format, ...);

第2引数として, fscanf関数に予めオーブンしたファイルのファイルポインタを渡す. 第2引数として, fscanf 関数に与える書式 format は scanf 関数と共通. 第3引数以降には, 変数のアドレスを渡すことに注意.

以下に例を挙げる.

char c; int i; double x; char str[100]; double v[100]; // のとき

fscanf(fp, "%c",  &c);
fscanf(fp, "%d",  &i);
fscanf(fp, "%lf", &x); // double 型は "%lf". 間違え易いので注意.
fscanf(fp, "%s",  str);  // 文字列専用の書式として %s がある. str は char str[100] の先頭要素のアドレス.
fscanf(fp, "%lf", &v[2]); // 各要素に & を付けると, その要素のアドレス.

fscanf(fp, "%c%d%lf",  &c, &i, &x); // 複数の変数に入力

15.5.1 ファイルから入力する

次のプログラムは, scanf の代わりに fscanf を用いて, キーボードから入力するのと同じようにファイルから入力する. このプログラムを実行する前に, このプログラムと同じフォルダ/ディレクトリに入力用データを書き込んだファイル plot.cnf を作成しておく必要がある. 入力の順序と入力の変数の型が一致していないと正しく動作しないので注意すること.

plot.cnfの内容の例:

red
5.0
4.2

plot.cnf の内容を変数に入力し, その値を画面に出力するプログラム:

#include <stdio.h>

int main(void)
{
	 FILE *fp;
	 char color[100];
	 double x, y;
	 
	 if ( (fp=fopen("plot.cnf", "r")) == NULL ) {
		  printf("Can not open plot.cnf\n");
	 } else {
		  fscanf(fp, "%s", color);
		  fscanf(fp, "%lf", &x);
		  fscanf(fp, "%lf", &y);

		  printf("color=%s\n", color);
		  printf("x=%f\n", x);
		  printf("y=%f\n", y);

		  fclose(fp);
	 }
	 
	 return 0;
}

1. fopen関数によりファイルを開く.
   • 第1引数にファイル名, 第2引数にモードを指定する. ファイルから入力を行う場合は, "r" (read) を指定する.
   • ファイルを開くことができたときは, fopen 関数はファイルを識別する FILE 構造体へのポインタ (ファイルポインタ) を返す. 以降では, このポインタによりファイルを指定して, ファイルからの入力を行う. そのために FILE *fp; を用意して, これに代入している.
   • C言語では, 代入式は代入した値を持つので, fp=fopen("plot.cnf", "r") は fopen 関数が返した値を持つ. このプログラムでは, この値が NULL 等しいときと等しくないときでif文で分岐している.

   • ファイルを開くことができなかったとき, fopen 関数は NULL を返す. NULL は <stdio.h>の中で定義されているアドレスである (#define NULL 0). このプログラムでは, ファイルを開くことができなかったときは, メッセージを出力して終了する.
2. fopen関数が NULL 以外のアドレスを返したとき, ファイルの操作が可能となる. ここでは, fscanf 関数によりファイルから変数に入力を行っている.
3. 読み込みが終了したら, fclose 関数によりファイルを閉じる.

15.5.2 ファイルから入力し描画する.

次のプログラムは, Configuration ファイル plot.cnf に保存してある 色(文字列), X座標(double), Y座標(double)を各変数に読み込み, SVGファイルに○を描画するプログラムである. 一度コンパイルしたのちは, plot.cnf の内容を書き換えると 再度コンパイルしないで実行するだけで, プロットされる点が変わることに注意すること.

gcc でコンパイルする場合は, svg1.h の中で数学関数を用いているので -lm オプションを付けておく.

$ emacs plot.cnf
$ emacs plot.c
$ gcc -Wall plot.c -lm
$ ./a.out
$ firefix drawing.svg
...                   # この後は コンパイルし直さなくても,
$ emacs plot.cnf      # plot.cnf を書き換え
$ ./a.out             # プログラムを実行すると
$ firefox drawing.svg # plot.cnf に応じて描画結果が変わっている.

【構造体の復習】

ATTRIBUTE circle_attrib;と宣言されている構造体のメンバーへの代入は, .を用いて circle_attrib.stroke = color;のように行う. ATTRIBUTE *circle_attrib;と構造体へのポインタとして宣言されているときは, 構造体のメンバーへの代入は, ->を用いてcircle_attrib->stroke = color;のように行う. (文法上は, ポインタとして宣言されている場合は (*circle_attrib).stroke = color; も正しい. ただし, この書き方は煩雑になることが多い思う.)

ATTRIBUTE circle_attrib;とATTRIBUTE *circle_attrib;は次のように使い分けるとよい. main関数などで, 実体を ATTRIBUTE circle_attrib;のようにして用意する. どこかで必ずこの実体の確保が必要である. 関数内で使用するときは, それへのポインタ(アドレス)を渡す. この場合呼び出される関数の仮引数は, ATTRIBUTE *circle_attrib となり, 関数内では ->を用いてメンバーを参照する.

【See also】

• 「ポインタは値を代入してから用いる」[10.4]
• 「大きなオブジェクトを関数の実引数とするとき」

【実行の手順】

1. 下のプログラム plot.c を作成する.
2. svg1.h を作成する.
3. plot.c をコンパイルし実行ファイルを作成する.
4. 実行ファイルと同じフォルダ/ディレクトリに上の plot.cnfを作成しておく.
5. 実行ファイルを実行する. これにより drawing.svg のような描画ファイル が作成される.
6. Firefox で描画ファイルを表示し, 描画結果を確認する.
7. plot.cnf の内容を変更したのち, 再度実行ファイルを実行する. 先ほどとは描画結果が異なることを確認する.

plot.c の内容:

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

int main(void)
{
	 FILE *fp;
	 char color[100];
	 double x, y;
	 
	 // プロットに関係する変数.
	 SVG svg_file;
	 PLOT plot;
	 PLOT_VIEW view = {100.0, 50.0, 800.0, 500.0};
	 PLOT_WORLD world = {-10.0, -10.0, 10.0, 10.0};
	 PLOT_AXIS axis = {"x", 5, "y", 5};
	 ATTRIBUTE circle_attrib   = {"blue", 1.0, "white"};

	 // プロット初期化
	 plot_file_open(&svg_file, "drawing.svg", 855.0, 605.0);
	 plot_init(&plot, &svg_file, &view, &world, &axis);

	 if ( (fp=fopen("plot.cnf", "r")) == NULL ) {
		  printf("Can not open plot.cnf\n");
	 } else {
		  fscanf(fp, "%s", color);
		  fscanf(fp, "%lf", &x);
		  fscanf(fp, "%lf", &y);
		  circle_attrib.stroke = color;
		  plot_circle_draw(&plot, x, y, 5.0, &circle_attrib);
		  fclose(fp);
	 }

	 // プロット終了
	 plot_file_close(&plot);
	 
	 return 0;
}

15.5.3 ファイル上のデータを入力してグラフにする

次のプログラムは, data.dat ファイルに記録されているデータを入力して, その値を画面に出力するプログラムである. fscanf 関数はファイルの最後まで入力が終了し次に入力する値が無いときに EOF を返す. EOF (End of File) は <stdio.h> の中で, 例えば #define EOF (-1)のように, 定義されている. 逆に, ファイル中から正しく入力が行われたときは, 入力された変数の数を返す.

このプログラムでは, fscanf 関数がEOFを返さない間, ファイルから入力を続ける. while の条件式の中で入力した xと y の値は, 条件が成立したとき実行するブロック内 ({ ... }) で使用できる. このプログラムでは, その値を画面に出力している.

このプログラムを実行する前に, 予め下に挙げる data.dat のようなデータファイルを作成しておくこと.

#include <stdio.h>

int main(void)
{
	 FILE *fp;
	 double x, y;
	 
	 if ( (fp=fopen("data.dat", "r")) == NULL ) {
		  printf("Can not open plot.cnf\n");
	 } else {
		  while(fscanf(fp, "%lf%lf", &x, &y) != EOF){
			   printf("x=%f, y=%f\n", x, y);
		  }
		  fclose(fp);
	 }
	 
	 return 0;
}

data.dat の内容の例:

0.7  5.0
3.8  6.4
5.8  5.7
8.4  5.2
9.3  3.0
13.1 3.8
13.2 6.6

ファイルから入力した値を画面に出力する代わりに, SVGファイルに描画するようにしたプログラムを示す.

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

int main(void)
{
	 FILE *fp;
	 double x, y;
	 
	 // プロットに関係する変数.
	 SVG svg_file;
	 PLOT plot;
	 PLOT_VIEW view = {100.0, 50.0, 800.0, 500.0};
	 PLOT_WORLD world = {-5.0, -5.0, 20.0, 15.0};
	 PLOT_AXIS axis = {"x", 5, "y", 5};
	 ATTRIBUTE circle_attrib   = {"blue", 1.0, "white"};

	 // プロット初期化
	 plot_file_open(&svg_file, "drawing.svg", 855.0, 605.0);
	 plot_init(&plot, &svg_file, &view, &world, &axis);

	 if ( (fp=fopen("data.dat", "r")) == NULL ) {
		  printf("Can not open data.dat\n");
	 } else {
		  while(fscanf(fp, "%lf%lf", &x, &y) != EOF){
			   plot_circle_draw(&plot, x, y, 2.0, &circle_attrib);
		  }
		  fclose(fp);
	 }

	 // プロット終了
	 plot_file_close(&plot);
	 
	 return 0;
}

描画結果が出力されている SVGファイル drawing.svg は, inkscape, Illustrator などで編集できるので試してみるとよい. inkscape で編集する場合は, 予め $ inkscape &のように起動しておき, File メニューから Open する. もしくは, ファイルブラウザ (nautilus) を起動して, 編集したいファイルを inkscape にドラッグする. (本来 $ inkscape drawing.svg & のようにファイル名を指定して起動できますが, 全学計算機システムでは設定がうまくいっていないようです.)

15.5.4 CSVファイルの読み込み

CSV ファイルは, コンマで各フィールドを区切ったデータファイルである. たとえば 次のような内容を持つ. これをfscanf で入力するときは, "%lf,%lf" のように書式に ,(コンマ)を入れおく.

data.csv の内容の例:

0.7,  5.0
3.8,  6.4
5.8,  5.7
8.4,  5.2
9.3,  3.0
13.1, 3.8
13.2, 6.6

#include <stdio.h>

int main(void)
{
	 FILE *fp;
	 double x, y;
	 
	 if ( (fp=fopen("data.csv", "r")) == NULL ) {
		  printf("Can not open plot.cnf\n");
	 } else {
		  while(fscanf(fp, "%lf, %lf", &x, &y) != EOF){
			   printf("x=%f, y=%f\n", x, y);
		  }
		  fclose(fp);
	 }
	 
	 return 0;
}

CSVファイルには, 次に挙げるようにテキストのフィールドは "で囲んだ上で, コンマで区切るものもある. 上のプログラムでは入力できないので修正する必要がある.

"x", "y"
0.7,  5.0
3.8,  6.4
5.8,  5.7
8.4,  5.2
9.3,  3.0
13.1, 3.8
13.2, 6.6

要望があれば続きを書きます...