プログラミング演習1 課題3 数値解析法

1. Cによる数値解析のための準備１

数学関数

Cには標準ライブラリの中に数学関数が用意されており、ヘッダ <math.h> をインクルードすることで使用可能となる。ヘッダ <math.h>内には下記の関数などが具体的に宣言されている。

練習1: 数学関数による様々な計算

以下のコードの(1)と(2)に値を入力して数学関数を使った計算を実行しよう。

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

int main(void)
{
    double x=__(1)__ ,y=__(2)__ ;

    printf("x = %f, y = %f\n", x, y);
    printf("x^y = %f\n", pow(x,y));
    printf("atan2(y/x) = %f\n", atan2(y, x));
    printf("sqrt(x*x+y*y) = %f\n", sqrt(pow(x,2) + pow(y,2)));
    return 0;
}

[実行結果（例）]
－－－－－－－－－－－－－－－－－－－－－－－
x = 3.000000, y = 4.000000
x^y = 81.000000
atan2(y/x) = 0.927295
sqrt(x*x+y*y) = 5.000000
－－－－－－－－－－－－－－－－－－－－－－－
x = 5.000000, y = 12.000000
x^y = 244140625.000000
atan2(y/x) = 1.176005
sqrt(x*x+y*y) = 13.000000   

数に関するユーティリティ

ヘッダ<stdlib.h>には、いくつかの有用な数値変換関数や整数乱数発生関数が宣言されている。

浮動小数点数に関する定数

ヘッダ <float.h>には浮動小数点数に関する様々な定数が定義されている。

練習2: 浮動小数点数のビットサイズと制限の確認

基本的に、数値計算には倍精度のdouble型を用いる。以下のコードを実行して、float型とdouble型の違いを確認しよう。

#include <stdio.h>
#include <float.h>

int main(void)
{
    printf("float: %u bit, double: %u bit\n", sizeof(float)*8, sizeof(double)*8);
    printf("FLT_MAX =     % e\n", FLT_MAX);
    printf("FLT_MIN =     % e\n", FLT_MIN);
    printf("DBL_MAX =     % e\n", DBL_MAX);
    printf("DBL_MIN =     % e\n", DBL_MIN);

    float f1  = 3.14159265358979323846;  // float
    double d1 = 3.14159265358979323846; // double
    printf("Value of f1 (float): %.20f\n", f1);
    printf("Value of d1 (double): %.20lf\n", d1);

    return 0;
}

[実行結果]
－－－－－－－－－－－－－－－－－－－－－－－
float: 32 bit, double: 64 bit
FLT_MAX =      3.402823e+038
FLT_MIN =      1.175494e-038
DBL_MAX =      1.797693e+308
DBL_MIN =      2.225074e-308
Value of f1 (float): 3.14159274101257324219
Value of d1 (double): 3.14159265358979311600

プログラムで間違いやすい分数式

(1) プログラムでは分数式のまま記述することができないので、割り算で代用することになる。その場合に分子・分母に括弧をつけなければならなくなるので注意すること。

たとえば、\(\displaystyle \frac{a}{b}\) は a/b でよいが、\(\displaystyle \frac{a+b}{c}\) は a+b/c ではなく、(a+b)/c のようにカッコを付ける。 同様に、\(\displaystyle \frac{a}{b+c}\) は a/b+c ではなく a/(b+c)、\(\displaystyle \frac{a}{bc}\) は a/b*c ではなく a/(b*c) となる。

(2) 整数どうしの割り算は整数除算を行うので、実数除算にするためには分子か分母を実数化する必要がある（※プログラミング入門Iでも説明）。整数値を実数化するには小数点をつければよい。

#include <stdio.h>

int main(void)
{
    double x1, x2, x3, x4;

    x1 = 1 / 2;   /* 整数除算 */
    x2 = 1./ 2;   /* 分子を実数化 */
    x3 = 1 / 2.;  /* 分母を実数化 */
    x4 = 1./ 2.;  /* 分子・分母を実数化 */

    printf("x1 = %f, x2 = %f, x3 = %f, x4 = %f\n", x1, x2, x3, x4);
    return 0;
}

[実行結果]
－－－－－－－－－－－－－－－－－－－－－－－
x1 = 0.000000, x2 = 0.500000, x3 = 0.500000, x4 = 0.500000
－－－－－－－－－－－－－－－－－－－－－－－

練習3: 型キャスト

上記の他に整数変数を実数化するための方法としては、キャスト演算子(double)などを用いた型変換がある。以下のコードを実行して、計算結果を確認してみよう。

#include <stdio.h>

int main(void)
{
    int a=1, b=2;
    double x1, x2, x3, x4;

    x1 = a / b;                  /* 整数除算 */
    x2 = (double)a / b;          /* 分子を実数化 */
    x3 = a / (double)b;          /* 分母を実数化 */
    x4 = (double)a / (double)b;  /* 分子・分母を実数化 */

    printf("a = %d, b = %d\n", a, b); 
    printf("x1 = %f, x2 = %f, x3 = %f, x4 = %f\n", x1, x2, x3, x4);
    return 0;
}

[実行結果]
－－－－－－－－－－－－－－－－－－－－－－－
a = 1, b = 2
x1 = 0.000000, x2 = 0.500000, x3 = 0.500000, x4 = 0.500000
－－－－－－－－－－－－－－－－－－－－－－－

for文による等間隔実数列の生成（復習）

関数をグラフで表すと、xとyの関係を理解する大きな手助けとなる。このため、数値解析では一定間隔に刻んだ数列を用いて計算することが多い。以下は、柴田望洋「新版 明解C言語 入門編」の第7章にある、 0 から1 まで 0.01 間隔の実数列を生成するプログラム List 7-9 を少し手直ししたものである。

/*
   0 から1 まで 0.01 間隔の実数列を生成する。
*/
#include <stdio.h>

int main(void)
{
    float x;

    for (x=0; x<=1; x+=0.01) printf(" x = %g\n", x);

    return 0;
}

[実行結果]
－－－－－－－－－－－－－－－－－－－－－－－
 x = 0
 x = 0.01
 x = 0.02
 ...
 x = 0.979999
 x = 0.989999
 x = 0.999999
－－－－－－－－－－－－－－－－－－－－－－－

このように、実数(単精度浮動小数点数)を直接、for文で増減させると、刻み幅(上記の場合は0.01)を2進数にしたときの丸め誤差(表現誤差)が累積されて次第に値を正確に刻んでくれなくなる場合が多い。このようなことを避けるためにも、for文は必ず整数で制御するべきである。

/*
   0 から1 まで 0.01 間隔の実数列を生成する。（ループを整数で制御）
*/
#include <stdio.h>

int main(void)
{
    int i;
    float x;

    for (i=0; i<=100; i++) {
        x = i / 100;
        printf(" x = %g\n", x);
    }
    return 0;
}

練習4: 等間隔実数列の生成

上記のプログラムは、前述の考え方を適用した「明解C言語」の第7章 List 7-10とほぼ同じものであるが、整数除算になっているため0しか出力されない。以下の実行結果になるように修正しよう。

[実行結果]
－－－－－－－－－－－－－－－－－－－－－－－
 x =  0
 x =  0.01
 x =  0.02
 ...
 x =  0.98
 x =  0.99
 x =  1
－－－－－－－－－－－－－－－－－－－－－－－