ニュートン法について

Question

ニュートン法を用いた問題で
方程式 f(x) = x^2 − 2 として，方程式 f(x) = 0 を
実行時に与えた初期値のもとで
Newton 法によって解くプログラムを作成したいのですが思い通りの結果になりません。
Newton 反復を 1 回行う関数の戻り値がよくわかっていないためだと思います。
変更すべき点を教えていただければ幸いです。

#include <stdio.h>
#include <math.h>
#define K_MAX 100
#define DX  1.0e-10
#define EPS 1.0e-10

void input(double *x);
double f(double x);
double df(double x);
double g(double x);
int newton(double *x);

int main(void) {
	double x;
	int k;

	input(&x);
	k = newton(&x);
	printf("\n<%d> f(%.16f) = %.16f\n",k,x,f(x));

	return 0;
}

void input(double *x) {
	printf("input x: ");
	scanf("%lf",x);
}

double f(double x) {
	return (x * x - 2);
}

double df(double x) {
	/* return (2 * x); */
	return ( (f(x + DX) - f(x)) / DX );
}

double g(double x) {
	return (x * x - 2);　
}

int dgtof(int m) {
    int ans = 1;

    while (m >= 10) {
        ans++;
        m = m / 10;
    }

    return ans;
}

int newton(double *x) {
	double newx, oldx = *x;
    int dgt = dgtof(K_MAX);
	int k;

	for (k = 1; k <= K_MAX; k++) {
		newx = g(oldx);
		printf("<%*d> g(%.16f) = %.16f\n",dgt,k,oldx,newx);
		if (fabs(newx - oldx) < EPS)
			break;
		oldx = newx;
	}
	*x = newx;

	return ((k < K_MAX) ? k : -1);
}

Comments

[mather]

「思い通りの結果になりません」期待値とどういう結果になったか書いてください。

[fana]

結果がどうの以前に，まずNewton法の実装自体が存在しないように見えるのですが…

Answer 1

ニュートン法の更新式がおかしいです。

newx = g(oldx);

↓

newx = oldx - f(oldx) / df(oldx);

サンプルコードを以下に記載します。

サンプルコード

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

#define K_MAX 100
#define DX 1.0e-10
#define EPS 1.0e-10

void input(double *init_x);
double f(double x);
double df(double x);
bool newton(double init_x, double *x, int *num_iters);

int main()
{
    double init_x;  // 初期値
    double x;       // 解
    int num_iters;  // 反復回数

    printf("input x: ");
    scanf("%lf", &init_x);

    if (newton(init_x, &x, &num_iters))
        printf("number of iteration: %d, x=%.15f\n", num_iters, x);
    else
        printf("newton method failed.");
}

/**
   @brief x^2 の値を計算する。
 */
double f(double x)
{
    return (x * x - 2);
}

/**
   @brief x^2 の導関数の値を計算する。
 */
double df(double x)
{
    return (f(x + DX) - f(x)) / DX;
}

/**
   @brief newton ニュートン法を実行する。
   @param init_x 初期値
   @param x 解
   @param num_iters 反復回数
 */
bool newton(double init_x, double *x, int *num_iters)
{
    double x1 = init_x, x2;
    int i;

    for (i = 1; i <= K_MAX; i++) {
        if (fabs(df(x1)) < EPS)
            return false;  // ニュートン法が失敗する場合

        x2 = x1 - f(x1) / df(x1);
        printf("%d: x=%.15f\n", i, x2);

        if (fabs(x2 - x1) < EPS) {
            *x = x2;
            *num_iters = i;
            return true;
        }

        x1 = x2;
    }

    return false;  // 収束しなかった場合
}

input x: 3
1: x=1.833333429863758
2: x=1.462121242835444
3: x=1.414998368736482
4: x=1.414213779259034
5: x=1.414213562373021
6: x=1.414213562373095
number of iteration: 6, x=1.414213562373095

Comments

[fana]

重箱の隅的な話ですが，df(x1)の値にチェックかけた方が良いです

[tiitoi]

コメントありがとうございます。チェックを入れてみました。

[fana]

dx(x1)の値次第でx2が遥か彼方の値に吹っ飛ぶかもしれないことを危惧するなら，もうちょっと厳し目な「dx(x1)の絶対値がとても小さいなら諦める」みたいな形でチェックするのかもしれませんね．（妥当な閾値の決め方が謎ですが…用途次第？）

[tiitoi]

たしかにご指摘のとおり、0近辺で値が発散するかもしれないので、<EPS としました。

[fana]

fabs( dx(x1) ) としないとまずいです

[fana]

× dx()
○ df()
ですね…　私のコメントもおかしい．

[EUROPEAN]

ありがとうございました！無事に求めていた結果がでました＾＾

Answer 2

更新式ですが、x_(n+1) = x_n - (f(x_n))/(f'(x_n)) を手計算で求め、算出する関数を実装すればよいです。このf(x)ですが、質問者さんが記載されているものとおなじです。

詳しくは、次のサイトを参考にしてください。

1. ニュートン法とは何か？？ニュートン法で解く方程式の近似解
2. 【微分法】接線の方程式の求め方