プログラミング演習1―課題4.1―

プログラミング演習1

― 課題4 数値解析法 ―

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

1.1 数学関数

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

sin(x) 三角関数 sin(x) を計算する。

cos(x) 三角関数 cos(x) を計算する。

tan(x) 三角関数 tan(x) を計算する。

asin(x) 逆三角関数 sin^-1(x) を計算する。(-π/2～π/2)

acos(x) 逆三角関数 cos^-1(x) を計算する。(0～π)

atan(x) 逆三角関数 tan^-1(x) を計算する。(-π/2～π/2)

atan2(y,x) 逆三角関数 tan^-1(y/x) を計算する。(-π～π)

sinh(x) 双曲線関数 sinh(x) を計算する。

cosh(x) 双曲線関数 cosh(x) を計算する。

tanh(x) 双曲線関数 tanh(x) を計算する。

exp(x) 指数関数 e^x を計算する。

log(x) 自然対数 ln(x) を計算する。(x > 0)

log10(x) 10 を底とする対数 log₁₀(x) を計算する。(x > 0)

pow(x,y) べき乗 x^y を計算する。

sqrt(x) 実数 x の平方根を求める。(x ≥ 0)

ceil(x) 実数 x の小数部を切り上げた整数を返す。

floor(x) 実数 x の小数部を切り捨てた整数を返す。

fabs(x) 絶対値|x|を求める。

ldexp(x,n) x・2ⁿ を計算する。

frexp(x,&exp) 実数 x を仮数部と指数部(2のべき乗) exp に分解する。

modf(x,&ip) 実数 x を整数部 ip と小数部に分ける。

fmod(x,y) 実数の剰余(x % y に相当)を求める。

練習1 functions.c (数学関数による様々な計算)
 1: #include 
 2: #include 
 3:
 4: int main(void)
 5: {
 6:  double x=_ ,y=_ ;
 7:
 8:     printf("x = %f, y = %f\n", x, y);
 9:     printf("X^Y = %f\n", pow(x,y));
10:     printf("SQRT_X = %f, SQRT_Y = %f\n", sqrt(x), sqrt(y));
11:     printf("EXP_X = %f, EXP_Y = %f\n", exp(x), exp(y)); 
12:     return 0;
13: }
◎x,y に値を入力して計算させてみよ。
[実行結果（例）]
－－－－－－－－－－－－－－－－－－－－－－－
x = 5.000000, y = 2.000000
X^Y = 25.000000
SQRT_X = 2.236068, SQRT_Y = 1.414214
EXP_X = 148.413159, EXP_Y = 7.389056
－－－－－－－－－－－－－－－－－－－－－－－
x = 10.000000, y = 5.000000
X^Y = 100000.000000
SQRT_X = 3.162278, SQRT_Y = 2.236068
EXP_X = 22026.465795, EXP_Y = 148.413159
－－－－－－－－－－－－－－－－－－－－－－－
x = 3.000000, y = 2.500000
X^Y = 15.588457
SQRT_X = 1.732051, SQRT_Y = 1.581139
EXP_X = 20.085537, EXP_Y = 12.182494

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

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

atof(s) 文字列 s を double型浮動小数点数に変換する。

atoi(s) 文字列 s を int型整数に変換する。

atol(s) 文字列 s を long型整数に変換する。

rand() 0～RAND_MAXの範囲の整数乱数を発生する。

srand(seed) rand() の種として seed をセットする。

abs(n) int型整数 n の絶対値 |n|を返す。

labs(n) long型整数 n の絶対値 |n|を返す。

練習2 rand10.c (乱数を発生する)
 1: #include 
 2: #include 
 3: 
 4: int main(void)
 5: {
 6:      int i, seed;
 7:      
 8:      seed=_;
 9:      srand(seed);
10: 
11:      printf("RAND_MAX = % d\n", RAND_MAX);
12:      printf("seed = %d\n", seed);
13:      for (i=0; i<10; i++) printf("%2d: %d\n", i, rand());
14:      return 0;
15: }
◎seedに値を入力して計算させてみよ。
[実行結果（例）]
－－－－－－－－－－－－－－－－－－－－－－－
RAND_MAX = 32767
seed = 1
0: 41
1: 18467
2: 6334
3: 26500
4: 19169
5: 15724
6: 11478
7: 29358
8: 26962
9: 24464
－－－－－－－－－－－－－－－－－－－－－－－
RAND_MAX = 32767
seed = 5
0: 54
1: 28693
2: 12255
3: 24449
4: 27660
5: 31430
6: 23927
7: 17649
8: 27472
9: 32640
－－－－－－－－－－－－－－－－－－－－－－－
RAND_MAX = 32767
seed = 10
0: 71
1: 16899
2: 3272
3: 13694
4: 13697
5: 18296
6: 6722
7: 3012
8: 11726
9: 1899

1.3 整数に関する制限

ヘッダ <limits.h> には整数型のサイズを表す定数が定義されている。

CHAR_BIT char型のビット数

CHAR_MAX char型の最大値

CHAR_MIN char型の最小値

INT_MAX int型の最大値

INT_MIN int型の最小値

LONG_MAX long型の最大値

LONG_MIN long型の最小値

SHRT_MAX short型の最大値

SHRT_MIN short型の最小値

UCHAR_MAX unsigned char型の最大値

UINT_MAX unsigned int型の最大値

ULONG_MAX unsigned long型の最大値

USHRT_MAX unsigned short型の最大値

練習3 chkint.c (標準的な整数のビットサイズと制限の確認)
 1: #include stdio.h;
 2: #include limits.h;
 3:
 4: int main(void)
 5: {
 6:     printf("int: %u bit, unsigned: %u bit\n", sizeof(int)*8, sizeof(unsigned)*8);
 7:     printf("INT_MAX = % d\n", INT_MAX);
 8:     printf("INT_MIN = % d\n", INT_MIN);
 9:     printf("UINT_MAX = %u\n", UINT_MAX);
10:     return 0;
11: }
◎上記プログラムを実行してみよ。
[実行結果]
－－－－－－－－－－－－－－－－－－－－－－－
int: 32 bit, unsigned: 32 bit
INT_MAX = 2147483647
INT_MIN = -2147483648
UINT_MAX = 4294967295
－－－－－－－－－－－－－－－－－－－－－－－

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

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

単精度

_FLT_RADIX 指数表現における基数

_FLT_ROUNDS 丸めのモード

FLT_DIG 精度の10進けた数

FLT_EPSILON 計算機イプシロン

FLT_MANT_DIG 浮動小数点数の有効数字部（_FLT_RADIXを基数とする）のけた数

FLT_MAX 表現可能な最大の有限浮動小数点数

FLT_MAX_EXP 表現可能な数の指数の最大値

FLT_MIN 最小の正規化された正の浮動小数点数

FLT_MIN_EXP 最小正規化数の指数の値

倍精度

_DBL_RADIX 指数表現における基数

_DBL_ROUNDS 丸めのモード

DBL_DIG 精度の10進けた数

DBL_EPSILON 計算機イプシロン

DBL_MANT_DIG 浮動小数点数の有効数字部（_DBL_RADIXを基数とする）のけた数

DBL_MAX 表現可能な最大の有限浮動小数点数

DBL_MAX_EXP 表現可能な数の指数部の最大値

DBL_MIN 最小の正規化された正の浮動小数点数

DBL_MIN_EXP 最小正規化数の指数部の値

練習4 chkflt.c (浮動小数点数のビットサイズと制限の確認)
 1: #include stdio.h;
 2: #include float.h;
 3:
 4: int main(void)
 5: {
 6:     printf("float: %u bit, double: %u bit\n", sizeof(float)*8, sizeof(double)*8);
 7:     printf("FLT_MAX =     % e\n", FLT_MAX);
 8:     printf("FLT_MIN =     % e\n", FLT_MIN);
 9:     printf("FLT_EPSILON = % e\n", FLT_EPSILON);
10:     printf("DBL_MAX =     % e\n", DBL_MAX);
11:     printf("DBL_MIN =     % e\n", DBL_MIN);
12:     printf("DBL_EPSILON = % e\n", DBL_EPSILON);
13:     return 0;
14: }
◎上記プログラムを実行してみよ。
[実行結果]
－－－－－－－－－－－－－－－－－－－－－－－
float: 32 bit, double: 64 bit
FLT_MAX = 3.402823e+038
FLT_MIN = 1.175494e-038
FLT_EPSILON = 1.192093e-007
DBL_MAX = 1.797693e+308
DBL_MIN = 2.225074e-308
DBL_EPSILON = 2.220446e-016

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

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

(2)整数どうしの割り算は整数除算を行うので、実数除算にするためには分子か分母を実数化する必要がある。
整数値を実数化するには小数点をつければよい。
 1: #include 
 2: 
 3: int main(void)
 4: {
 5:     double x1, x2, x3, x4;
 6: 
 7:     x1 = 1 / 2;   /* 整数除算 */
 8:     x2 = 1./ 2;   /* 分子を実数化 */
 9:     x3 = 1 / 2.;  /* 分母を実数化 */
10:     x4 = 1./ 2.;  /* 分子・分母を実数化 */
11:
12:     printf("x1 = %f, x2 = %f, x3 = %f, x4 = %f\n", x1, x2, x3, x4);
13:     return 0;
14: }
[実行結果]
－－－－－－－－－－－－－－－－－－－－－－－
x1 = 0.000000, x2 = 0.500000, x3 = 0.500000, x4 = 0.500000
－－－－－－－－－－－－－－－－－－－－－－－
他に、整数変数を実数化するための方法としては、キャスト演算子　(double)　
を用いた型変換がある。
練習5 fraction.c (プログラムで間違いやすい分数式)
 1: #include 
 2: 
 3: int main(void)
 4: {
 5:    int a=1, b=2;
 6:    double x1, x2, x3, x4;
 7: 
 8:     x1 = a / b;                  /* 整数除算 */
 9:     x2 = (double)a / b;          /* 分子を実数化 */
10:     x3 = a / (double)b;          /* 分母を実数化 */
11:     x4 = (double)a / (double)b;  /* 分子・分母を実数化 */
12: 
13:     printf("a = %d, b = %d\n", a, b); 
14:     printf("x1 = %f, x2 = %f, x3 = %f, x4 = %f\n", x1, x2, x3, x4);
15:     return 0;
16: }
◎上記プログラムを実行してみよ。
[実行結果]
－－－－－－－－－－－－－－－－－－－－－－－
a = 1, b = 2
x1 = 0.000000, x2 = 0.500000, x3 = 0.500000, x4 = 0.500000
－－－－－－－－－－－－－－－－－－－－－－－

1.6 for文による等間隔実数列の生成

数値解析では一定間隔に刻んだ実数列を用いて計算することが多い。以下は、柴田望洋「新版明解C言語入門編」の第7章にある、 0 から1 まで 0.01 間隔の実数列を生成するプログラム List 7-9 を少し手直ししたもの(forflt.c)である。
 1: /*
 2:     0 から1 まで 0.01 間隔の実数列を生成する。
 3:  */
 4: #include <stdio.h>
 5:
 6: int main(void)
 7: {
 8:     float x;
 9:
10:     for (x=0; x<=1; x+=0.01) printf(" x = % g\n", x);
11:
12:     return 0;
13: }
[実行結果]
－－－－－－－－－－－－－－－－－－－－－－－
 x = 0
 x = 0.01
 x = 0.02
 ...
 x = 0.979999
 x = 0.989999
 x = 0.999999
－－－－－－－－－－－－－－－－－－－－－－－
このように、実数(単精度浮動小数点数)を直接、for文で増減させると、刻み幅(上記の場合は0.01)を2進数にしたときの丸め誤差(表現誤差) が累積されて次第に値を正確に刻んでくれなくなる場合が多い。このようなことを避けるためにも、for文は必ず整数で制御するべきである。以下のプログラム(forflt2.c)は、この考え方を適用した「明解C言語」の第7章 List 7-10と同じものである。
 1: /*
 2:     0 から1 まで 0.01 間隔の実数列を生成する。（整数で制御）
 3:  */
 4: #include <stdio.h>
 5:
 6: int main(void)
 7: {
 8:     int i;
 9:     float x;
10:
11:     for (i=0; i<=100; i++) {
12:         x = i / 100.0;
13:         printf(" x = % g\n", x);
14:     }
15:     return 0;
16: }
また、実数(浮動小数点数)宣言を double にすることで表現誤差が非常に小さくなる。 float を用いていて結果がおかしいときには double にするだけで直ることがある。
練習6 (等間隔実数列の生成)
(1)上記プログラム（forflt2.c）を実行してみよ。

[実行結果]
－－－－－－－－－－－－－－－－－－－－－－－
 x =  0
 x =  0.01
 x =  0.02
 ...
 x =  0.98
 x =  0.99
 x =  1
－－－－－－－－－－－－－－－－－－－－－－－
(2)上記プログラムの12行めの‘x = i / 100.0' 分母を‘100’，‘100.’にしてみよ。

[実行結果：‘100’の場合]
－－－－－－－－－－－－－－－－－－－－－－－
 x =  0
 x =  0
 x =  0
 ...
 x =  0
 x =  0
 x =  1
－－－－－－－－－－－－－－－－－－－－－－－

[実行結果：‘100.’の場合]
－－－－－－－－－－－－－－－－－－－－－－－
 x =  0
 x =  0.01
 x =  0.02
 ...
 x =  0.98
 x =  0.99
 x =  1
－－－－－－－－－－－－－－－－－－－－－－－
(3)1つめのプログラム(forflt.c)の8行目の‘float’を‘double’にしてみよ。

[実行結果]
－－－－－－－－－－－－－－－－－－－－－－－
 x =  0
 x =  0.01
 x =  0.02
 ...
 x =  0.97
 x =  0.98
 x =  0.99
－－－－－－－－－－－－－－－－－－－－－－－
【目次】 | 【2.】| 【3.】

sin(x)	三角関数 sin(x) を計算する。
cos(x)	三角関数 cos(x) を計算する。
tan(x)	三角関数 tan(x) を計算する。
asin(x)	逆三角関数 sin^-1(x) を計算する。(-π/2～π/2)
acos(x)	逆三角関数 cos^-1(x) を計算する。(0～π)
atan(x)	逆三角関数 tan^-1(x) を計算する。(-π/2～π/2)
atan2(y,x)	逆三角関数 tan^-1(y/x) を計算する。(-π～π)
sinh(x)	双曲線関数 sinh(x) を計算する。
cosh(x)	双曲線関数 cosh(x) を計算する。
tanh(x)	双曲線関数 tanh(x) を計算する。
exp(x)	指数関数 e^x を計算する。
log(x)	自然対数 ln(x) を計算する。(x > 0)
log10(x)	10 を底とする対数 log₁₀(x) を計算する。(x > 0)
pow(x,y)	べき乗 x^y を計算する。
sqrt(x)	実数 x の平方根を求める。(x ≥ 0)
ceil(x)	実数 x の小数部を切り上げた整数を返す。
floor(x)	実数 x の小数部を切り捨てた整数を返す。
fabs(x)	絶対値\|x\|を求める。
ldexp(x,n)	x・2ⁿ を計算する。
frexp(x,&exp)	実数 x を仮数部と指数部(2のべき乗) exp に分解する。
modf(x,&ip)	実数 x を整数部 ip と小数部に分ける。
fmod(x,y)	実数の剰余(x % y に相当)を求める。

atof(s)	文字列 s を double型浮動小数点数に変換する。
atoi(s)	文字列 s を int型整数に変換する。
atol(s)	文字列 s を long型整数に変換する。
rand()	0～RAND_MAXの範囲の整数乱数を発生する。
srand(seed)	rand() の種として seed をセットする。
abs(n)	int型整数 n の絶対値 \|n\|を返す。
labs(n)	long型整数 n の絶対値 \|n\|を返す。

CHAR_BIT	char型のビット数
CHAR_MAX	char型の最大値
CHAR_MIN	char型の最小値
INT_MAX	int型の最大値
INT_MIN	int型の最小値
LONG_MAX	long型の最大値
LONG_MIN	long型の最小値
SHRT_MAX	short型の最大値
SHRT_MIN	short型の最小値

UCHAR_MAX	unsigned char型の最大値
UINT_MAX	unsigned int型の最大値
ULONG_MAX	unsigned long型の最大値
USHRT_MAX	unsigned short型の最大値

単精度
_FLT_RADIX	指数表現における基数
_FLT_ROUNDS	丸めのモード
FLT_DIG	精度の10進けた数
FLT_EPSILON	計算機イプシロン
FLT_MANT_DIG	浮動小数点数の有効数字部（_FLT_RADIXを基数とする）のけた数
FLT_MAX	表現可能な最大の有限浮動小数点数
FLT_MAX_EXP	表現可能な数の指数の最大値
FLT_MIN	最小の正規化された正の浮動小数点数
FLT_MIN_EXP	最小正規化数の指数の値

倍精度
_DBL_RADIX	指数表現における基数
_DBL_ROUNDS	丸めのモード
DBL_DIG	精度の10進けた数
DBL_EPSILON	計算機イプシロン
DBL_MANT_DIG	浮動小数点数の有効数字部（_DBL_RADIXを基数とする）のけた数
DBL_MAX	表現可能な最大の有限浮動小数点数
DBL_MAX_EXP	表現可能な数の指数部の最大値
DBL_MIN	最小の正規化された正の浮動小数点数
DBL_MIN_EXP	最小正規化数の指数部の値