Last active
November 17, 2018 16:30
-
-
Save tsuzu/a8bf34f85963aafd97b7f267a88cb668 to your computer and use it in GitHub Desktop.
Cプロ第2回大レポート
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#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; | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment