プログラミング演習 第 6 回関数 II

目次

1 目標 1 2 sin 関数の作成 1 3 台形公式による定積分の計算 (変数のキャスト) 3

3.1 台形公式と誤差 . . . . . . . . . . . . . . . . . . . . . . . . . . 6 4 問題 8

1 目標

• 関数 sin をつくる.

• 台形公式による定積分の計算

2 sin 関数の作成

sin 関数は C 言語の数学ライブラリーに定義されておりますが, 勉強のために自分の sin 関数 (= mysin) を作ってみましょう.

x = 0 におけるマクローリン展開は,

sin x = x − \frac{x^{3}}{3!} + \frac{x^{5}}{5!} − \frac{x^{7}}{7!} + ⋯

です. よって, mysin(double) のアルゴリズムは, 再帰性を使って書くと,

double mysin(double x)
{
double xx, kai, diff, sum;
int p;
p = 3;
sum = x;                          /* 第一項までの和 */
xx = x * x;
kai = 1;                          /* 分母の初期化 */
while(1)                          /* 無限ループ */
{
  x = -x*xx;                      /* 分子の計算 */
  kai *= (p - 1)*p;
  p += p + 2;
  diff = xx/kai                  /* sin 関数の p 項 */
  sum += diff;                   /* 第 p 項までの和 */
  if(fabs(diff) < eps)           /* 求めた第 p 項が */
   {                             /* 許容誤差 eps 以下で */
     break;                      /* ループを抜け出る */
     }
  }
    return sum;
}

作った mysin() 関数とライブラリーに登録されている sin() 関数を比較するプログラムを書きましょう. 但し, 範囲は 0 ≤ x ≤ 3.0 です.

ファイル名は, 6-1.c です.

#include<stdio.h>
#include <math.h>
#define eps 10e-11
double mysin(double);
main()
{
   double x, y;
   printf("\tx\t\tsin x\t\tmysin x\n");
   for(x=0.0; x < 3.1; x += 0.1)
   {
     printf("\t%3.1f\t%12.10f\t%12.10f\n", x, sin(x), mysin(x));
     }
}
double mysin(double x)
{
    double xx, kai, diff, sum;
    int p;
    p = 3;
    sum = x;                          /* 第一項までの和 */
    xx = x * x;
    kai = 1;                          /* 分母の初期化 */
    while(1)                          /* 無限ループ */
    {
      x = -x*xx;                      /* 分子の計算 */
      kai *= (p - 1)*p;
      p = p + 2;
      diff = x/kai;                  /* sin 関数の p 項 */
      sum += diff;                   /* 第 p 項までの和 */
      if(fabs(diff) < eps)           /* 求めた第 p 項が */

       {                             /* 許容誤差 eps 以下で */
         break;                      /* ループを抜け出る */
         }
      }
        return sum;
}

main 内の printf 関数で使われている%12.10f の意味は、実数を 12 桁で表示し, そのうち小数点以下 10 桁とする, ということです.

while 文内の if 文にある許容誤差 eps は main() 関数の前で定義されています. 特定の数や言葉などを使う場合は

     define (定義名) (定義したいもの)

であらかじめ定義しておくと, エラーを防げると同時に, 間違いを訂正しやすいです.^{1}

上の場合, 小数点以下 10 桁を表示しますから, 許容誤差として 10^{−11} でよろしいわけです.

cc 6-1.c -lm でコンパイルした後./a.out で実行さすると, データが多いためウインドウないで全て見れないかも知れません. そのような場合, リダイレクションをつかってデータをファイルに保存すればよろしいです. 例えば,

     %./a.out > sindata

とすれば, データは, sindata に保存されていますので, less コマンドでファイルの中身を見ればよろしいです.

3 台形公式による定積分の計算 (変数のキャスト)

計算機による数値計算の 1 例として, 台形公式を用いた定積分の計算を取り上げます. ここでは, ∫_{0}^{1} \frac{4 dx}{1 + x^{2}} = π を取り上げましょう.

以下は「高等学校数学Ｃ（数研出版）」から抜粋です．

連続な関数 y = f(x) は，区間 [a, b] で常に f(x) ≥ 0 であるとする．図 1 のように，この区間を n 等分して，分点を小さい方から順に

a = x_{0} < x_{1} < ⋯ < x_{n} = b

とし，次のようにおく．

h = \frac{b − a}{n}, y_{i} = f(x_{i})

小区間 [x_{i−1}, x_{i}] において，区間の幅 h を高さとし， y_{i−1} と y_{i} をそれぞれ，

1 間違いを訂正することを「デバックをする」といいます

図 1:

上低，下低とする台形の面積を S_{i} とすると

S_{i} = \frac{h}{2}(y_{0} + y_{1}), S_{2} = \frac{h}{2}(y_{1} + y_{2}), ⋯, S_{n} = \frac{h}{2}(y_{n−1} + y_{n})

これらの台形の面積の総和を T_{n} とすると

T_{n} = S_{1} + S_{2} + ⋯ + S_{n}

は， n が十分に大きいとき，定積分 ∫_{a}^{b}f(x) dx の近似値を与える.

したがって，定積分の近似式を与える次の台形公式が得られる．

∫_{a}^{b} f(x) dx ≒\frac{b − a}{2n}\{y_{0} + 2(y_{1} + y_{2} + ⋯ + y_{n−1}) + y_{n}\}

この台形公式は，常に f(x) ≥ 0 でない場合でも成り立つ．

それでは例を考えましょう. プログラムの中で, daikei と言う関数が台形公式を与え, それの引数は順に, 積分区間の右端, 左端, 分割の個数, と関数です.

math.h を include して, main() の中の daikei という関数の引数の fun の部分を例えば sin に変えると, このままで ∫_{0}^{1} sin xdx の近似値計算ができます (この場合当然数学関数もリンクしなければなりません). 他の引数に付いても同様ですし, 関数 fun の記述を変える事により, いろいろな定積分の近似計算ができます.

プログラムで注意して欲しいのは, 下から 9 行目の fun の引数の計算 a + interval*i の部分です. ここでは a, interval は double 型, i は int 型で宣言されています. このような異なる型の変数の間の計算は, 変数の型変換が自動的に行われます. 今の場合 i を double 型に変換して計算が実行されます.

型変換の規則は, char < int < double で与えられます. すなわち, 計算時にはこれらの内, より右側の型に値の型を合わせて計算されます.

/* ファイル名 6-2.c */
#include <stdio.h>
double daikei(double a, double b, int n);
double fun(double x);
main()
{
     printf("Pi is about %1.20lf.\n", daikei(0.0, 1.0, 10000));
}
double daikei(double a, double b, int n)
{
     double sum=fun(a);
     double interval=(b-a)/n;   /* The length of the small interval */
     int i;
     for (i=1; i < n ; i++)      /* add to (n-1)-th value */
               sum = sum + 2.0*fun(a+interval*i);
     return interval*(sum + fun(b))/2.0;
}

図 2:

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

上のプログラムは, π を近似的に求める方法の一つです.

3.1 台形公式と誤差

定積分に関する次の部分積分法の公式

∫_{α}^{β} h(x)g′(x) dx = [h(x)g(x)]_{α}^{β} − ∫_{α}^{β} h′(x)g(x) dx  (1)

を用いると，台形公式による数値積分の誤差について調べることができます．

定積分 ∫_{α}^{β} f(x) dx は，図 2 の斜線の面積で近似するとその誤差は，次のようになります．

e = \frac{β − α}{2}\{f(α) + f(β)\} − ∫_{α}^{β} f(x) dx  (2)

部分積分法の公式を用いて，まず次の等式が成り立つことを示そう．

e = \frac{1}{2} ∫_{α}^{β}(x − α)(β − x)f′′(x) dx  (3)

公式 (1) において， h(x) として (x − α)(β − x)， g(x) として f′(x) をとると， (3) の右辺は次のようになります．

\frac{1}{2}[(x − α)(β − x)f′(x)]_{α}^{β} − ∫_{α}^{β}(−2x + α + β)f′(x) dx
= \frac{1}{2} ∫_{α}^{β}(2x − α − β)f′(x) dx

再び公式 (1) において h(x) として 2x − α − β， g(x) として f(x) をとると，上の式は次のようになります.

\frac{1}{2} [(2x − α − β)f(x)]_{α}^{β} − ∫_{α}^{β} f(x) dx

これは (2) の右辺に等しいから， (3) の等式が成り立ちます．

前節において，定積分 ∫_{x_{i−1}}^{x_{i}} f(x) dx を S_{i} で近似したときの誤差は

E_{i} = \frac{x_{i} − x_{i−1}}{2} f(x_{i−1}) + f(x_{i}) − ∫_{x_{i−1}}^{x_{i}} f(x) dx

これは，公式の (2) で α = x_{i−1}, β = x_{i} とおいた式であるから (3) により

E_{i} = \frac{1}{2} ∫_{x_{i−1}}^{x_{i}}(x − x_{i−1})(x_{i} − x)f′′(x) dx

ここで

\frac{1}{2} ∫_{x_{i−1}}^{x_{i}}(x − x_{i−1})(x_{i} − x) dx = \frac{h^{3}}{12} ただし， h = x_{i} − x_{i−1}

であるから，小区間 [x_{i−1}, x_{i}] において， f′′(x) の最大値を M_{i}，最小値を m_{i} とすると，次の公式が成り立ちます．

\frac{h^{3}}{12} m_{i} ≤ E_{i} ≤ \frac{h^{3}}{12} M_{i}

関数 f′′(x) が連続であるとき，分割の個数 n を大きくとって，小区間の幅 h を小さくすると

m_{i} ≒ f′′(x_{i}) ≒ M_{i} したがって E_{i} ≒\frac{h^{3}}{12} f′′(x_{i})

となる．よって，定積分 ∫_{a}^{b} f(x) dx の台形公式による近似値の誤差は

E = T_{n} − ∫_{a}^{b} f(x) dx = ∑_{i=1}^{n} E_{i} ≒\frac{h^{2}}{12} ∑_{i=1}^{n} f′′(x_{i})h

ここで，関数 f′′(x) の区分求積法による数値積分を考えると

∑_{i=1}^{n} f′′(x_{i})h ≒ ∫_{a}^{b} f′′(x) dx = [f′(x)]_{a}^{b} = f′(b) − f′(a)

ゆえに

E ≒\frac{h^{2}}{12} \{f′(b) − f′(a)\}

このことから，この誤差は，小区間の幅の平方 h^{2} にほぼ比例し， |f′(b) − f′(a)| が小さければ，誤差も小さいことがわかります．

4 問題

1. 関数mycos を定義し, 0.0 ≤ x ≤ 3.0 の範囲の値を求めるプログラムを書き, 出力せよ. 但し, 許容誤差は 10^{−10} の範囲で.

2. (レポート問題)

Euler 関数 ϕ(n) は, 2 以上の整数 n に対して 1 から n までの数で, n と互いに素な素数の個数として定義される. 例えば, ϕ(3) = 2, ϕ(4) = 2 となる. オイラー関数を C 言語の関数と定義して, n を入力したときその値を出力するプログラムを書け.

3. (レポート問題)

台形公式を用いて log 2 = ∫_{0}^{1} \frac{1}{1+x}dx を計算するプログラムを書き, 区間の分割数を 10, 100, 1000 としたときの値を比較せよ.

締切は 5 月 30 日とする.