C言語でニュートン法を実装してみた

ニュートン法とは

ニュートン法については、ここでは割愛します。
http://qiita.com/PlanetMeron/items/09d7eb204868e1a49f49
この方の記事がわかりやすいかも

微分

導関数の定義に従って微分します。

# define DX 0.00000001

double f(double x) {
    return - 2*pow(x,4) + 3*pow(x,3) + 4*x - 3;
}

// 微分を求める
double fd(double x) {
    return (f(x+DX) - f(x)) / DX;
}

dxは0に近い適当な値を入れておきます。
fはxに対しての返り値で、fdはその微分を求めます。

接線を求める

dfを傾きとして、

y = df(currentX) * (x - currentX) + f(currentX)

currentXは接線を求めたい位置のx座標

x切片を求める

// x切片を返す、y=ax+bの形、ax+b=0のとき
double searchInterceptX(double a, double b) {
    return -b/a;
}

y=ax+bの切片は簡単にもとまりますね！

完成したファイル

Newton.c

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

# define DX 0.00000001

// 入力した関数
double f1(double x) {
    return - 2*pow(x,4) + 3*pow(x,3) + 4*x - 3;
}
double f2(double x) {
    return -2 * sin(x) + exp(x+1) + x - 10;
}

// 微分を求める
double fd(double (*func)(double), double x) {
    return (func(x+DX) - func(x)) / DX;
}

// x切片を返す、y=ax+bの形、ax+b=0のとき
double searchInterceptX(double a, double b) {
    return -b/a;
}

// ニュートン法で次のXを求める
double nextX(double (*func)(double), double currentX) {
    double dif = fd(func, currentX);
    return searchInterceptX(dif, func(currentX) - dif*currentX);
}

// 解いて結果を返す
// 本当は誤差評価すべき
double solve(double (*func)(double), double initialX) {
    double currentX = initialX;
    for(int i=0; i<10; i++) {
        printf("%lf\n", currentX);
        currentX = nextX(f1, currentX);
    }
    return currentX;
}

int main() {
    printf("f1\n");
    printf("result = %lf\n", solve(f1, 1));

    printf("\n");

    printf("f2\n");
    printf("result = %lf\n", solve(f2, M_PI));
}

関数ポインタを使って実装しました。
solve関数の第二引数は、ニュートン法の初期x座標になっています。求めたい解の近くの値にしておいてください。

誤差評価はすべきですが、今回はループ回数のみにしました