プログラミング演習1

― レポート課題3 ―

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

【ﾚﾎﾟｰﾄ遅延とならないための方法】
★ﾚﾎﾟｰﾄは、3つの各課題において“課題3用ﾃﾝﾌﾟﾚｰﾄ文書”記載の項目(課題概要・目的、実現方法、ｱﾙｺﾞﾘｽﾞﾑ、PAD、変数の説明、ｿｰｽとｺﾒﾝﾄ、実行結果とその説明、考察等)を、”全て”満たす必要があるので、完成させるﾚﾎﾟｰﾄの枚数はかなり多くなる(目安：A4で30頁程度)。
⇒したがって、「ﾚﾎﾟｰﾄ作成は講義内で終わらせる」つもりでどんどん進めること。講義終了後、まとめてﾚﾎﾟｰﾄ作成する方法はあまりお勧めできない。なお課題3-2と3-3の内容は第3回目に扱うが、“数値解析”の講義で既に学んでいるので、着手は可能であると思われる。

---

●課題3−1　発生系列が異なる疑似乱数の生成

以下は、疑似乱数を発生させて、それらの平均値と標準偏差を求めるプログラムである。
プログラム中の�@〜�Hの穴を埋めて完成させよ。
（☆穴埋めの時に全角スペース等が残っているとエラーとなるので注意すること。）

(1) seedに“100”を設定し、srand関数を同じseedで5回呼び、1〜100までの疑似乱数を
　　10個発生させること。更に、それらの疑似乱数の平均値と標準偏差を毎回求めて示すこと。
(2) seedに“5”、“10”、“50”を設定した場合の実行結果をそれぞれ示すこと。
(3) 横軸を設定したseed（5、10、50、100）、縦軸を1回目の疑似乱数の平均値とする
　　グラフ[棒グラフ]をExcelを用いて作成して示すこと※。

 1: #include <stdlib.h>
 2: #include <stdio.h>
 3: #include <math.h>
 4: #define �@＿�A＿
 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回目] ", �B＿, count1);
15:         srand( �C＿ );                              /*ヒント1：�@，�B，�Cは同じもの*/ 
16:         for (count2 = 0; count2 < �D＿; count2++){
17:             random = rand()%�E＿;　　　　　 　　　　/*ヒント2：1〜100までの表し方*/ 　 
18:             sum    += random;
19:             ssqu   += �F＿; 　　　　　　　　　　　　/*ヒント3：ssquは乱数の二乗和*/ 
20:             printf("%3d ", random);
21:         }
22:         mean = �G＿ / �H＿;
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) [挿入]→[縦棒グラフ]→[2D縦棒]を選択し、「ウィンドウ」をあげる。
(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) 初期区間 [x₁, x₂]の値と精度eps を入力できること。
(2) 負の領域の解も求められること^（※1）。

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

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

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

※2：初期区間を入力するようにしたプログラムでは、誤差を極端に小さくした場合には、 do〜while の反復が止まらなくなる。

※3：x₁, x₂ の大きさがx₁ > x₂ の場合や交点が0かまたは2個以上あるような区間から開始した場合に 誤った結果を求めてしまう。

[オプション課題]

上記【1】【2】の機能を保ったまま、以下の２つの要求を満たすプログラムを作成し、 “f(x) = x² - 5”, “eps = 1e-10”， “x₁ = 1, x₂ = 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 ² - 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 ² + (y/2)² = 1, y ² = x ² + 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() の引数に関数引数用の定義(関数ポインタ)を追加する。

---