Cプロ第2回大レポート

074-1W183088.c

#include <stdio.h>
#include <math.h>

#define NEWTON_MAX_TIMES 20
#define EPS (1e-6)
#define DIFFERENTIAL_H (1e-4)

/* f(x) を定義する */
double Fx(double x, int pm /* 1 or -1 */) {
    return x + pm * sqrt(1 - x * x) + 1./2;
}

/* f(x)で符号がプラス */
double FxPlus(double x) { /* f(x) を定義する */
    return Fx(x, 1);
}
/* f(x)で符号がマイナス */
double FxMinus(double x) {
    return Fx(x, -1);
}

/* 一次関数を用いxからyを算出する */
double CalcY(double x) {
    return x + 1./2;
}

/* doubleのswap関数 */
void SwapDouble(double* x, double* y) {
    double tmp = *x;
    *x = *y;
    *y = tmp;
}

/* intのswap関数 */
void SwapInt(int* x, int* y) {
    int tmp = *x;
    *x = *y;
    *y = tmp;
}

/* 微分値を求める関数 */
double CalcDifferential(double (*calc)(double), double x, double h) {
    return (calc(x + h) - calc(x)) / h;
}

/* ニュートン法で近似値を求める */
double NewtonRaphson(double (*calc)(double), double init, double eps, double h, int max_times, int *counter) {
    *counter = 0;
    int i;
    double x = init, delta;

    for(i = 0; i < max_times; ++i) {
        delta = calc(x) / CalcDifferential(calc, x, h);
        ++*counter;

        x -= delta;
        if(isnan(x))
            return 0. / 0.; /* 微分値が0でdeltaおよびxの値がNaNとなる */

        if (fabs(delta) < eps)
            return x;
    }
    
    return 0. / 0.; /* max_times回で収束しなかった */ 
}

/* 二分法で近似値を求める */
double Dichotomy(double (*calc)(double), double a, double b, double eps, int* counter) {
    *counter = 0;
    double c, fc;
    
    if(calc(a) > 0 && calc(b) < 0) {
        /* f(a) < 0 && f(b) > 0となるようにswapする */
        SwapDouble(&a, &b);
    }else if(calc(a) > 0 || calc(b) < 0) { /* a<=c<=bにf(c)=0となるcが存在しない */
        return 0.;
    }
    
    while(c = (a + b) / 2, fabs(a - b)  >= 2*eps) {
		fc = calc(c);
        ++*counter;
        if(fc > 0) {
            b = c;
        }else if(fc < 0) {
            a = c;
        }else if(fc == 0) {
            break;
        }
    }
    
    return c;
}

/* フォーマットに沿って出力させる */
void Printer(char const* str, double x, double y, int counter) {
    printf("%s part: ", str);
    
    if(isnan(x)) {
        printf("not converge\n");
    }else {
        printf("(%d times) : (%.6f, %.6f)\n", counter, x, y);
    }
}

int main(void) {
    int counter[2];
    double x[2];
    
    /* ニュートン法で2つのf(x)でf(x)=0となる解を求める */
    x[0] = NewtonRaphson(FxPlus, 0, EPS, DIFFERENTIAL_H, NEWTON_MAX_TIMES, &counter[0]);
    x[1] = NewtonRaphson(FxMinus, 0, EPS, DIFFERENTIAL_H, NEWTON_MAX_TIMES, &counter[1]);
    
    /* UpperとLowerを調整する(一方がNaNの時はLowerがNaNとなるように) */
    if(x[0] < x[1] || isnan(x[0])) {
        SwapDouble(&x[0], &x[1]);
        SwapInt(&counter[0], &counter[1]);
    }
    
    /* ニュートン法の結果を出力 */
    printf("Newton-Raphson method:\n");
    Printer("Upper", x[0], CalcY(x[0]), counter[0]);
    Printer("Lower", x[1], CalcY(x[1]), counter[1]);
 
	puts("");

    /* 二分法で2つのf(x)でf(x)=0となる解を求める */
    x[0] = Dichotomy(FxPlus, -1., 1., EPS, &counter[0]);
    x[1] = Dichotomy(FxMinus, -1., 1., EPS, &counter[1]);
    
    /* UpperとLowerを調整する(一方がNaNの時はLowerがNaNとなるように) */
    if(x[0] < x[1] || isnan(x[0])) {
        SwapDouble(&x[0], &x[1]);
        SwapInt(&counter[0], &counter[1]);
    }
    
    /* 二分法の結果を出力 */
    printf("Bisection method:\n");
    Printer("Upper", x[0], CalcY(x[0]), counter[0]);
    Printer("Lower", x[1], CalcY(x[1]), counter[1]);
    
    return 0; 
}