きんじかい

近似解を求める。

2分法
逆線形補間法

ニュートン・ラプソン法

一般的に、計算速度はニュートン・ラプソン法、逆線形補間法、2分法の順に長くなります

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

double NewtonRaphson(double x, double eps, int max, int *n);
double Bisection(double x1, double x2, double eps, int max, int *n);
double InverseLiner(double x1, double x2, double eps, int max, int *n);
double func(double x);
double dfunc(double x);

int main(void)
{
    double eps = 0.00001;
    double x1 = -1.0, x2 = 1.0, ret;
    int n, max = 1000;

    /*
    printf("x1 = ");
    scanf("%lf", &x1);
    printf("x2 = ");
    scanf("%lf", &x2);
    printf("ε = ");
    scanf("%lf", &eps);
    printf("最大反復回数 = ");
    scanf("%d", &max);
    */

    ret = Bisection(x1, x2, eps, max, &n);
    printf("--Bisection method--\n");
    printf("反復回数 : %d\n", n);
    printf("x = %lf\n\n", ret);

    ret = InverseLiner(x1, x2, eps, max, &n);
    printf("--Inverse Liner Interpolation method--\n");
    printf("反復回数 : %d\n", n);
    printf("x = %lf\n\n", ret);

    ret = NewtonRaphson(x1, eps, max, &n);
    printf("--Newton Raphson method--\n");
    printf("反復回数 : %d\n", n);
    printf("x = %lf\n\n", ret);

    return 0;
}

double NewtonRaphson(double x, double eps, int max, int *n)
{
    int i;
    double next_x;

    for (i = 0;i < max;i++) {
        next_x = x - func(x) / dfunc(x);
        if (fabs(x - next_x) < eps) {
            *n = i + 1;
            return x;
        }
        x = next_x;
    }
    *n = -1;
    return x;
}

double Bisection(double x1, double x2, double eps, int max, int *n)
{
    double x;
    int i;

    for (i = 0;i < max;i++) {
        x = (x1 + x2) / 2.0;
        if (func(x)*func(x1) < 0)
            x2 = x;
        else
            x1 = x;
        if (fabs(x1-x2) < eps) {
            *n = i + 1;
            return x;
        }
    }
    *n = -1;
    return x;
}

double InverseLiner(double x1, double x2, double eps, int max, int *n)
{
    double x;
    int i;

    for (i = 0;i < max;i++) {
        x = (x1*func(x2) - x2*func(x1)) / (func(x2) - func(x1));
        if (func(x)*func(x1) < 0)
            x2 = x;
        else
            x1 = x;
        if (fabs(x1-x2) < eps) {
            *n = i + 1;
            return x;
        }
    }
    *n = -1;
    return x;
}

double func(double x)
{
    return x - cos(x);
}

double dfunc(double x)
{
    return 1 + sin(x);
}