プログラミング演習1


― 課題4 数値解析 ―

宿題4-1の解答

要求された機能を満たす2分法のプログラムを示す。

(1) 精度、初期区間の入力機能追加
(1-1) 初期区間 [x1, x2] の値と 精度 eps を入力できるプログラム (bisec1_1.c)。 
/*
    bisec1_1.c: 初期区間、精度入力バージョン
 */
#include <stdio.h>

double f(double x)
//  f(x) = x^2 - 2
{
    return x * x - 2.0;
}

int main(void)
{
    int n=0;
    double x1=1, x2=2, x3, err;
    double eps = 1.0e-6;
    char s[128];

    fprintf(stderr, "# x1, x2 = "); fgets(s, 128, stdin);
    sscanf(s, "%lf%lf", &x1, &x2);
    fprintf(stderr, "# eps = "); fgets(s, 128, stdin);
    sscanf(s, "%lf", &eps);
    printf("# n, x, err\n");
    do {
        n++;
        x3 = (x1 + x2) / 2.0;
        if (f(x3) > 0.0) x2 = x3;
        else             x1 = x3;
        err = x2 - x1;
        printf("%4d, % .15e, % .15e\n", n, x3, err);
    } while (err >= eps);
    printf("\n# sqrt(2) = % .15e\n", (x1 + x2) / 2.0);
    return 0;
}
(1-2) 負の領域の解も求められるプログラム (bisec1_2.c)。 
/*
    bisec1_2.c: 負領域可能バージョン
 */
#include <stdio.h>

double f(double x)
//  f(x) = x^2 - 2
{
    return x * x - 2.0;
}

int main(void)
{
    int n=0;
    double x1=1, x2=2, x3, err;
    double eps = 1.0e-6;
    char s[128];

    fprintf(stderr, "# x1, x2 = "); fgets(s, 128, stdin);
    sscanf(s, "%lf%lf", &x1, &x2);
    fprintf(stderr, "# eps = "); fgets(s, 128, stdin);
    sscanf(s, "%lf", &eps);
    printf("# n, x, err\n");
    do {
        n++;
        x3 = (x1 + x2) / 2.0;
        if (f(x1)*f(x3) < 0.0) x2 = x3;
        else                   x1 = x3;
        err = x2 - x1;
        printf("%4d, % .15e, % .15e\n", n, x3, err);
    } while (err >= eps);
    printf("\n# sqrt(2) = % .15e\n", (x1 + x2) / 2.0);
    return 0;
}
(2) エラー処理の追加
(2-1) 反復が止まらない場合にエラーメッセージ “*** 反復回数が100回を超えました。***”を表示して 有限回(N = 100)で止まるプログラム (bisec2_1.c)。 
/*
    bisec2_1.c: 反復回数制限バージョン
 */
#include <stdio.h>
#define N 100

double f(double x)
//  f(x) = x^2 - 2
{
    return x * x - 2.0;
}

int main(void)
{
    int n=0;
    double x1=1, x2=2, x3, err;
    double eps = 1.0e-6;
    char s[128];

    fprintf(stderr, "# x1, x2 = "); fgets(s, 128, stdin);
    sscanf(s, "%lf%lf", &x1, &x2);
    fprintf(stderr, "# eps = "); fgets(s, 128, stdin);
    sscanf(s, "%lf", &eps);
    printf("# n, x, err\n");
    do {
        n++;
        if (n > N) {
            fprintf(stderr, "*** 反復回数が100回を超えました。***");
            break; // →[付録1]
        }
        x3 = (x1 + x2) / 2.0;
        if (f(x1)*f(x3) < 0.0) x2 = x3;
        else                   x1 = x3;
        err = x2 - x1;
        printf("%4d, % .15e, % .15e\n", n, x3, err);
    } while (err >= eps);
    printf("\n# sqrt(2) = % .15e\n", (x1 + x2) / 2.0);
    return 0;
}
(2-2) 入力値が x1 > x2 でも解が求められるプログラム (bisec2_2.c)。 
/*
    bisec2_2.c: 初期区間逆入力可能バージョン
 */
#include <stdio.h>
#include <math.h>
#define N 100

double f(double x)
//  f(x) = x^2 - 2
{
    return x * x - 2.0;
}

int main(void)
{
    int n=0;
    double x1=1, x2=2, x3, err;
    double eps = 1.0e-6;
    char s[128];

    fprintf(stderr, "# x1, x2 = "); fgets(s, 128, stdin);
    sscanf(s, "%lf%lf", &x1, &x2);
    fprintf(stderr, "# eps = "); fgets(s, 128, stdin);
    sscanf(s, "%lf", &eps);
    printf("# n, x, err\n");
    do {
        n++;
        if (n > N) {
            fprintf(stderr, "*** 反復回数が100回を超えました。***");
            break;
        }
        x3 = (x1 + x2) / 2.0;
        if (f(x1)*f(x3) < 0.0) x2 = x3;
        else                   x1 = x3;
        err = fabs(x2 - x1);
        printf("%4d, % .15e, % .15e\n", n, x3, err);
    } while (err >= eps);
    printf("\n# sqrt(2) = % .15e\n", (x1 + x2) / 2.0);
    return 0;
}
(2-3) 初期区間の数値が異常な場合 (例えばf (x3)とf (x1), f (x2)の符号がすべて同じ場合) は入力エラーメッセージ“*** 初期区間の数値が異常です。***” を表示して終了するプログラム (bisec2_3.c)。 
/*
    bisec2_3.c: 異常入力検出バージョン
 */
#include <stdio.h>
#include <math.h>
#define N 100

double f(double x)
//  f(x) = x^2 - 2
{
    return x * x - 2.0;
}

int main(void)
{
    int n=0;
    double x1=1, x2=2, x3, err;
    double eps = 1.0e-6;
    char s[128];

    fprintf(stderr, "# x1, x2 = "); fgets(s, 128, stdin);
    sscanf(s, "%lf%lf", &x1, &x2);
    fprintf(stderr, "# eps = "); fgets(s, 128, stdin);
    sscanf(s, "%lf", &eps);
    printf("# n, x, err\n");
    do {
        n++;
        if (n > N) {
            fprintf(stderr, "*** 反復回数が100回を超えました。***");
            break;
        }
        x3 = (x1 + x2) / 2.0;
        if      (f(x1)*f(x3) < 0.0) x2 = x3;
        else if (f(x3)*f(x2) < 0.0) x1 = x3;
        else {
            fprintf(stderr, "*** 初期区間が異常です。***\n");
            return 1;  // →[付録1]
        }
        err = fabs(x2 - x1);
        printf("%4d, % .15e, % .15e\n", n, x3, err);
    } while (err >= eps);
    printf("\n# sqrt(2) = % .15e\n", (x1 + x2) / 2.0);
    return 0;
}
(3) ソースのさらなる改良
(3-1) 二分法のアルゴリズムの部分(bisec.cの20行め〜28行めに相当する部分) を bisection() と名づけ、 main関数から独立させたプログラム (bisec3_1.c)。 
/*
    bisec3_1.c: サブルーチン独立バージョン
 */
#include <stdio.h>
#include <math.h>
#define N 100

double f(double x)
//  f(x) = x^2 - 2
{
    return x * x - 2.0;
}

double bisection(double x1, double x2, double eps)
{
    int n=0;
    double x3, err;

    printf("# n, x, err\n");
    do {
        n++;
        if (n > N) {
            fprintf(stderr, "*** 反復回数が100回を超えました。***");
            break;
        }
        x3 = (x1 + x2) / 2.0;
        if      (f(x1)*f(x3) < 0.0) x2 = x3;
        else if (f(x3)*f(x2) < 0.0) x1 = x3;
        else {
            fprintf(stderr, "# 初期区間が異常です。\n");
            exit(1);  // →[付録1]
        }
        err = fabs(x2 - x1);
        printf("%4d, % .15e, % .15e\n", n, x3, err);
    } while (err >= eps);
    return (x1 + x2) / 2.0;
}

int main(void)
{
    double x1=1, x2=2;
    double eps = 1.0e-6;
    char s[128];

    fprintf(stderr, "# x1, x2 = "); fgets(s, 128, stdin);
    sscanf(s, "%lf%lf", &x1, &x2);
    fprintf(stderr, "# eps = "); fgets(s, 128, stdin);
    sscanf(s, "%lf", &eps);
    printf("\n# sqrt(2) = % .15e\n", bisection(x1, x2, eps));
    return 0;
}
(3-2) 方程式の左辺の関数 f (x) の名前をmain関数側で自由に 指定できるように bisection の引数に関数引数を用いた定義を追加したプログラム (bisec3_2.c)。 
/*
   bisec3_2.c: サブルーチン関数引数バージョン
 */
#include <stdio.h>
#include <math.h>
#define N 100

double f(double x)
//  f(x) = x^2 - 2
{
    return x * x - 2.0;
}

double bisection(double (*f)(double), double x1, double x2, double eps)
{
    int n=0;
    double x3, err;

    printf("# n, x, err\n");
    do {
        n++;
        if (n > N) {
            fprintf(stderr, "*** 反復回数が100回を超えました。***");
            break;
        }
        x3 = (x1 + x2) / 2.0;
        if (f(x1)*f(x3) < 0.0) x2 = x3;
        else if (f(x3)*f(x2) < 0.0) x1 = x3;
        else {
            fprintf(stderr, "# 初期区間が異常です。\n");
            exit(1);
        }
        err = fabs(x2 - x1);
        printf("%4d, % .15e, % .15e\n", n, x3, err);
    } while (err >= eps);
    return (x1 + x2) / 2.0;
}

int main(void)
{
    double x1=1, x2=2;
    double eps = 1.0e-6;
    char s[128];

    fprintf(stderr, "# x1, x2 = "); fgets(s, 128, stdin);
    sscanf(s, "%lf%lf", &x1, &x2);
    fprintf(stderr, "# eps = "); fgets(s, 128, stdin);
    sscanf(s, "%lf", &eps);
    printf("\n# sqrt(2) = % .15e\n", bisection(f, x1, x2, eps));
    return 0;
}
【目次】 | 【1.】 | 【2.】 | 【3.】 | 【付録1】 【付録2】