第 9 章 ファイルとの入出力 45
9.4 提出課題
第9章 ファイルとの入出力 53
第9章 ファイルとの入出力 54 なお，入力ファイルの例は以下のとおり．一行目にデータ数と時間増分(計測時間幅), そのあとに温度のデータを並べて記述すること．
20 2.0 18.0 22.0 24.0 38.0 31.0 32.0 . . . 28.0 25.0 22.0 20.0
• 提出先: [email protected]
• 提出期限: 12月22日(月) 17:00
• ファイル名: 08XXXXXXX-No9-HW.c (XXXXは学生番号)
55
第 10 章 応用プログラミング
今週は，C言語によるプログラミングの総まとめとして，数値計算を中心とした応用プ ログラミングについて学ぶ．具体的なプログラミング事例として，数値積分と方程式の求 解を取り上げ，例題演習を通して実践的なプログラミング技法を習得する．
10.1 数値積分
10.1.1 もっとも原始的な数値積分 ( 長方形近似 )
数値積分とは,積分区間の被積分関数 f(x)と x 軸に挟まれた面積を細分化した棒グラ フで近似し,それらの棒グラフの面積の総和(Σ) をとることである．以下のような手順に より，積分を数値的に求めるための概念について学ぼう．
a b
k=1 k=2k=3
k=N
x y
Fig.1 長方形近似による数値積分
上の図に示されているように，関数f(x)を区間[a, b]について積分する場合について考 える．まず，区間[a, b]をN等分し，それぞれの区間の中点の座標xkについて考える．
xk =a+ 1
2∆x(2k−1), (k = 1,2, ..., N) (10.1) ただし，
∆x= b−a
N (10.2)
第10章 応用プログラミング 56 である．また，上記の座標xkにおける関数値は，
yk =f(xk) (10.3)
であり，各区間の長柱の面積は，幅×高さで求められるから，
Sk=f(xk)∆x (10.4)
となる．結果的に，区間[a, b]における面積の総和は以下のように求められる．
S =
∑N k=1
Sk =
∑N k=1
f(xk)∆x (10.5)
以上の手順に基づいて，プログラムを作成してみよう．
例題10–1 区分的な面積を長方形により近似して，関数y= 3x2 の面積を区間[1,2]
において求めよ．分割数Nはキーボードより入力できるようにせよ．
解答例10–1
1: #include <stdio.h>
2: #include <stdlib.h>
3: #include <math.h>
4:
5: int main(void){
6: int k, N;
7: double a, b; /*積分区間*/
8: double x; /*積分座標*/
9: double dx; /*きざみ幅*/
10: double s; /*面積 */
11:
12: printf("分割数は．．．\n");
13: scanf("%d", &N);
14: a = 1.0;
15: b = 2.0;
16: dx = ( b - a ) / (double)N;
17:
18: for( k = 1, s = 0.0 ; k <= N ; k++ ){
19: x = a + 0.5*dx*(2.0*(double)k - 1.0);
20: s = s + 3.0 * x * x * dx ; 21: }
22: printf("面積は %f\n", s);
23: return(0);
24: }
※ Nやkは整数であるから，doubleに型変換されている点に注意すること！
第10章 応用プログラミング 57
10.1.2 台形公式
a b
k=1 k=2
k=N
x y
k=3
Fig.2 台形公式による数値積分
先程の例では，各区間を長方形により近似した．しかしながら，長方形の近似では各区 間における長方形の面積と実際の面積の差＝誤差がやや大きく，分割数を比較的大きくと らないと正解値が得られない．そこで次に，各区間を台形により近似して面積を求める方 法を考えてみよう．先の場合と同様に，区間[a, b]をN等分すると，それぞれの区間(k番 目の区間)における左側の座標xLk と右側の座標xRk は以下のようになる．
xLk =a+ (k−1)∆x, (k= 1,2, ..., N) (10.6) xRk =a+k∆x, (k = 1,2, ..., N) (10.7) また，上記の座標における関数値は，
ykL=f(xLk), ykR=f(xRk) (10.8) であり，各区間の面積は，台形の公式より上低(yLk)＋下低(ykR)×高さ(∆x)÷2で求め られるから，
Sk= 1
2{f(xLk) +f(xRk)}∆x (10.9) と求められる．結果的に，区間[a, b]における面積の総和は以下のようになる．
S =
∑N k=1
Sk =
∑N k=1
1
2{f(xLk) +f(xRk)}∆x (10.10) さらに，総和に関する計算を簡略化すると，結果的に台形公式による区間[a, b]におけ る面積の総和は以下のように計算できることがわかる．
S =
[1 2
{
f(a) +f(b)}+
N∑−1 k=1
f(xk)
]
∆x, xk =k∆x (10.11)
第10章 応用プログラミング 58 例題10–2 区分的な面積を台形により近似して，関数y = 3x2 の面積を区間[1,2]に
おいて求めよ．なお，分割数Nはキーボードより入力できるようにせよ．
解答例10–1
1: #include <stdio.h>
2: #include <stdlib.h>
3: #include <math.h>
4: int main(void){
5: int k, N;
6: double a, b; /*積分区間*/
7: double xl, xr; /*積分座標*/
8: double yl, yr; /*関数値 */
9: double dx; /*きざみ幅*/
10: double s; /*面積 */
11: double f( double ); /*関数f(x)のプロトタイプ宣言*/
12: printf("分割数は．．．\n");
13: scanf("%d", &N);
14: a = 1.0;
15: b = 2.0;
16: dx = ( b - a ) / (double)N;
17: for( k = 1, s = 0.0 ; k <= N ; k++ ){
18: xl = a + dx*((double)k - 1.0);
19: xr = a + dx*(double)k;
20: yl = f(xl);
21: yr = f(xr);
22: s = s + 0.5 * ( yl + yr ) * dx;
23: }
24: printf("面積は %f\n", s);
25: return(0);
26: }
27: /* 関数 f(x) = 3*x*x */
28: double f ( double x ){
29: return ( 3.0 * x * x );
30: }
上記のプログラムは，式(10)に従い，各区間毎に台形の面積を求めてから総和を計算 する場合の例である．前ページの式(11)に示されているように，‘台形公式’を適用して計 算を簡略化する場合は，上記プログラムの17行目から23行目を，以下のように書き換え るとよい．
s = 0.5 * ( f(a) + f(b) ) * dx;
for( k = 1 ; k < N ; k++ ){
xl = a + (double)k * dx;
s = s + f(xl) * dx ;
}
第10章 応用プログラミング 59