C++でニュートン法を美しく書く その２（再帰処理で書く）

前回は
C++でニュートン法を美しく書く　その１
で，ニュートン法を書いてみました．このときはwhileループを使って書きましたが，同じことを再帰処理を使って書いてみたいと思います．

#include <iostream>

auto solve_by_newton_method = []( auto func , double x ,double error , double delta ) {

    double x1 = x;
    double x2 = x1 + delta;

    double y1 = func(x1);
    double y2 = func(x2);

    double diff_x = (y2 - y1) / (x2 - x1);
    double next_x = x1 - y1 / diff_x;

    if (fabs(y1 - 0) < error) {
        return x1;
    }
    double answer = solve_by_newton_method(func, next_x, error , delta);

    return answer;
};

int main()
{
    double x;
    double y;
    double xini = 5.0;//初期値

    auto func = [](double x) {
        return x * x - 5.0;
    };

    x = solve_by_newton_method( func ,xini , 0.00001 , 0.00001);
    y = func(x);

    std::cout << "x = " << x << std::endl;
    std::cout << "y = " << y << std::endl;

    return 0;
}

はい，書けました．
次に，すこしかゆいところに手を届かせてみます．ニュートン法をはじめとする数値計算は，最終的に収束してくれれば良いのですが，収束せずに何度もループを回ってしまうことが良くあります．
そういう時に，無限ループに陥ってしまうとよろしくありませんので，ループ回数を制限することがあります．それで，ループ回数の上限に達してしまった場合は，これまでのニュートン法のステップの中で一番解に近いものを出力するようにしてみます．

一応やりたいことは，

• ループの上限を設定できる
• ループ上限に達した場合は一番解に近いものを出力

です．

#include <iostream>

auto solve_by_newton_method = [](auto func, double x, double error, double delta , int counter , int max_loop ) {
    double x1 = x;
    double y1 = func(x1);
    double diff_x = ( func(x1+delta)  - y1) / ((x1 + delta) - x1);
    double next_x = x1 - y1 / diff_x;

    if ( fabs(y1 - 0.0) < error) {
        //収束判定に収まった場合は，再帰処理の底になる．
        return x1;
    } else if ( counter > max_loop) {
        //ループ上限に達した場合は，再帰処理の底になる．
        std::cout << "loop over" << std::endl;
        return x1;
    }else {
        //どれにも当てはまらなかったら再帰処理を続ける
        double answer = solve_by_newton_method(func, next_x, error, delta, counter + 1, max_loop);

        //収束しなかったanswerが戻ってくる場合もあるので，この階層の再帰処理の結果と比較して，解に近い方を返す．
        if (fabs(y1 - 0.0) < fabs(func(answer) - 0.0)) {
            return x1;
        }
        else {
            return answer;
        }
    }
};

int main()
{
    double x;
    double y;
    double xini = 5.0;//初期値

    //解きたい関数の定義
    auto func = [](double x) {
        return x * x + 0.1;
    };

    x = solve_by_newton_method(func, xini, 0.0001, 0.00001, 1, 1000);
    y = func(x);

    std::cout << "x = " << x << std::endl;
    std::cout << "y = " << y << std::endl;

    return 0;
}

こんな感じになります．最大1000回の再帰呼び出しにしました．
ちなみに解く方程式は，

y=x^2+0.1

に変更しました．この方程式は実数解を持たないので，ニュートン法は解を見つけようとして何度も何度もループを繰り返します．