プログラミング演習1


― レポート課題3 ―

以下の課題を解き、レポートとして提出せよ。

【ﾚﾎﾟｰﾄ遅延に御注意を！】ﾚﾎﾟｰﾄは、3つの課題(3-1, 3-2, 3-3)において“課題3用ﾃﾝﾌﾟﾚｰﾄ文書”記載の項目(課題概要・目的、実現方法、ｱﾙｺﾞﾘｽﾞﾑ、PAD、変数の説明、ｿｰｽとｺﾒﾝﾄ、実行結果とその説明、考察等)を、”全て”満たす必要があります！したがって、完成ﾚﾎﾟｰﾄの枚数は非常に多くなります(目安：A4で20頁～)。
⇒その結果、「ﾚﾎﾟｰﾄ作成は講義内で終わらせる」つもりでどんどん進めないと終わりません！講義終了後、まとめてﾚﾎﾟｰﾄ作成に着手する方法はあまりお勧めしません。なお課題3-2と3-3の内容は第3回目に扱いますが、“数値解析”で既に学習しているので、着手可能と思われます。頑張って下さい！
★なお、各課題のプログラムのヒントは、講義内で説明したプログラムテクニックの随所にみられますので、課題を読んでから講義ページを見返すと良いでしょう
●課題3－1　発生系列が異なる疑似乱数の生成
以下は、疑似乱数を発生させて、それらの平均値と標準偏差を求めるプログラムである。 プログラム中の①～⑨の穴を埋めて完成させよ。 （☆穴埋めの時に全角スペース等が残っているとエラーとなるので注意すること。）
(1) seedにある値*を設定し、srand関数を同じseedで5回呼び、1～100までの疑似乱数を10個発生させるプログラムを作成すること。また、疑似乱数の平均値と標準偏差を毎回求めて示すこと。
*seedには“5”、“10”、“50”、“100”を設定し、それぞれの場合について実行結果を示すこと。
(2) 横軸を設定したseed（5、10、50、100）、縦軸を1回目の疑似乱数の平均値とする 　　グラフ[棒グラフ]をExcelを用いて作成して示すこと※。 
 1: #include <stdlib.h>
 2: #include <stdio.h>
 3: #include <math.h>
 4: #define ①＿②＿
 5: 
 6: int main(void)
 7: {
 8:    int count1, count2, random;
 9:    double sum, ssqu, mean;
10: 
11:     for (count1 = 1; count1 <= 5; count1++){
12:         sum = ssqu = mean = 0.0;
13:         random = 0; 
14:         printf("seed=%2d  [%2d回目] ", ③＿, count1);
15:         srand( ④＿ );                              /*ヒント1：①，③，④は同じもの*/ 
16:         for (count2 = 0; count2 < ⑤＿; count2++){
17:             random = rand()%⑥＿;　　　　　 　　　　/*ヒント2：1～100までの表し方*/ 　 
18:             sum    += random;
19:             ssqu   += ⑦＿; 　　　　　　　　　　　　/*ヒント3：ssquは乱数の二乗和*/ 
20:             printf("%3d ", random);
21:         }
22:         mean = ⑧＿ / ⑨＿;
23:         printf("Mean:%4.2lf, SD:%4.2lf\n", mean, sqrt(fabs(ssqu/count2 - pow(mean,2))));
24:    }
25:    return 0;
26: }
【※作図のポイント】
(i) seed, Mean, SDの対応表を作成する。
(ii) [挿入]→[2-D縦棒]→[集合縦棒]を選択し、「ウィンドウ」をあげる。
(iii)そのウィンドウ上で右クリック→[データの選択]→データソースの選択ウィンドウの[追加]を選択する。
(iV) 系列名：疑似乱数、系列値→Meanの4個のデータを選択[OK]する。
(V) データソースの選択ウィンドウ[編集]を選択し、seedの4つの設定値(5,10,50,100)を選択する。
(Vi) 図の右上の“＋”の[誤差範囲]→[その他の誤差範囲オプション]を選択する。
　　 誤差範囲の書式設定の[ユーザ設定]のボタンをクリックし、[値の指定]を選択する。
　　 正の誤差の値、負の誤差の値共に、SDの4個のデータを選択[OK]する。
(Vii)図の右上の“＋”の軸ラベルを選択し、縦軸、横軸ラベルを記載するなど、図としての完成度をあげる。
 <補足資料>
[オプション課題]
(1)上記のプログラムを改良し、time関数（ヘッダをインクルード）を用いて seedに“現在の時刻”を設定し、srand関数を同じseedで10回呼び、 1～1000までの疑似乱数を5個発生させるプログラムを作成せよ。 なお、実行結果は3種類以上示すこと。
(2)毎回 [○回目の中で] 異なる乱数を発生させるために、プログラムの改良すべき点を説明し、 実装せよ。
●課題3－2　二分法
以下の要求を満たす二分法のプログラムを作成せよ。 なお、ソースリストや実行結果などは【1】【2】毎に示せばよい。

【1】二分法のプログラムに以下の２つの仕様を追加したプログラムを作成し、 “初期区間 [-2, -1]”, “eps = 1e-10”とした場合の実行結果を示すこと。

(1) 初期区間 [x1, x2]の値と精度eps を入力できること。
(2) 負の領域の解も求められること（※1）。

※1：x < 0 では関数が右下がりになるので、区間を選択する条件が変わる。 if 文を用いて条件を切り替える方法もあるが、たとえば f (x1)*f (x3) の符号で判断すれば、 余計な場合分けをしなくて済む。

【2】更に以下の３つの仕様を追加してプログラムを改良し、“eps = 1e-16” とした場合、“[x1, x2] = [2, 6], [-6, -2], [-2, 6], [-6, 2]” とした場合、“x1 = 2, x2 = 1”とした場合の各々について実行結果を示し、 それをExcelによりグラフ描画して示すこと。なお、図のx軸はn、y軸はxとerrとすること。

(1) 反復が止まらない場合はエラーメッセージ “*** 反復回数が100回を超えました。***” を表示して有限回(N = 100)で止まること（※2）。
(2) 入力値が x1 > x2でも解が求められること（※3）。
(3) 初期区間の数値が異常な場合(例えばf (x3)とf (x1), f (x2)の符号がすべて同じ場合) は入力エラーメッセージ“*** 初期区間の数値が異常です。***”を表示して終了すること（※3）。

※2：初期区間を入力するようにしたプログラムでは、誤差を極端に小さくした場合には、 do～while の反復が止まらなくなる。
※3：x1, x2 の大きさがx1 > x2 の場合や交点が0かまたは2個以上あるような区間から開始した場合に 誤った結果を求めてしまう。

[オプション課題]
上記【1】【2】の機能を保ったまま、以下の２つの要求を満たすプログラムを作成し、 “f(x) = x2 - 5”, “eps = 1e-10”， “x1 = 1, x2 = 3” とした場合の 実行結果を示すこと。

(1)二分法のアルゴリズムの部分(bisec.cの20行め～28行めに相当する部分)を bisection() と名づけ、 main関数から独立させること。 
 1: double bisection(double x1, double x2, double eps)
 2: {
 3:     …………………………
 4:     二分法のアルゴリズム
 5:     …………………………
 6:     return (x1 + x2) / 2.0;
 7: }

(2)方程式の左辺の関数 f (x) の名前をmain関数側で自由に指定できるように 
bisection の引数に関数引数用の定義(関数ポインタ)を追加すること。

 1: double bisection(double (*f)(double), double x1, double x2, double eps)
 2: {
 3:     …………………………
 4:     二分法のアルゴリズム
 5:     …………………………
 6:     return (x1 + x2) / 2.0;
 7: }

●課題3－3　ニュートン法
以下の要求を満たすニュートン法のプログラムを作成せよ。 ニュートン法のプログラムに以下の４つの改良を加えて、a = 3 および a = 5 の場合 についての実行結果を示し、それをExcelによりグラフ描画して示すこと。 なお、図のx軸はn、y軸はxとerrとすること。

(1) 式(2.1)を用い、方程式の左辺関数f (x) (= x 2 - a ) および その導関数 f '(x) を独立させること。ただし、入力変数 a はグローバル変数としてよい。 導関数 f '(x) は、プログラムでは例えばdf (x) と名づけておくこととする。

(2)ニュートン法のアルゴリズムの部分(newton.c の 20行目～28行目)を newton() と名づけ、 main関数から独立させること。 
 1: double newton(double x, double eps)
 2: {
 3:     …………………………………
 4:     ニュートン法のアルゴリズム
 5:     …………………………………
 6:     return x;
 7: }

(3) 方程式の左辺の関数 f (x) とその導関数の名前をmain関数側で自由に指定できるように、
newton() の引数に関数引数用の定義(関数ポインタ)を追加すること。

 1: double newton(double x, double (*f)(double), double (*df)(double), double eps)
 2: {
 3:     …………………………………
 4:     ニュートン法のアルゴリズム
 5:     …………………………………
 6:     return x;
 7: }

(4) a の数値をコマンド行から入力できるようにすること。
[オプション課題]
2変数のニュートン法のプログラム newton2.c を以下の３つの要求を満たすように改良し、 x 2 + (y/2)2 = 1, y 2 = x 2 + 1 の場合について第1象限の解を求めること。

(1)関数f1(), f2()を配列 f[1], f[2]、関数J11()～J22()を2次元配列 J[1][1]～J[2][2] に割り当てることにして、方程式の左辺関数f (x)とそのヤコビ行列 J (x)を各々1つの サブルーチン(func(), Jcalc())にまとめること。ただし、関数の戻り値で渡していた値は 引数で受け渡すことにし、手続きの戻り値は終了コード(整数)とする。 
 1: int func(double x, double y, double f[])
 2: {
 3:     ……………………
 4:     f[1] = …………;
 5:     f[2] = …………;
 6:     return 0;
 7: }
また、Cでの配列の定義は0から始まることに注意し、一つ多めに宣言しておく。 配列の数(ND)をヘッダのインクルードの下で以下のように定数マクロで宣言しておくとよい。 
   #define ND 3
2次元配列を引数でうまく渡せなかった人は次善の策として1次元配列にして渡すか、 大域変数で引き渡すこと。

(2) 変数x, yを配列x[1], x[2]、変数ux, uyをu[1], u[2]に割り当ててサブルーチン等の 引数の数を減らし、n 変数に移行しやすくする。

(3) 方程式の左辺の関数 f (x) とその導関数の名前をmain関数側で自由に指定できるように、 newton2() の引数に関数引数用の定義(関数ポインタ)を追加する。