Help us understand the problem. What is going on with this article?

C++でニュートン法を美しく書く その1

More than 1 year has passed since last update.

誰もが一度は通る道かもしれませんが,数値計算の例として,ニュートン法を書いてみたいと思います.なぜc++かというと,私がc++でいろいろと数値計算をしているからです.数値計算というとやっぱりFortranと言う気もしますが,c++はそこそこ実行速度も速いし,いろいろライブラリも整ってるし,オブジェクト指向はそこそこ便利と言うこともあって,私はC++は数値計算には結構向いていると思います.あともちろん,自分が所属している組織のルールで言語指定っていうのもあるしね.

それで,一応例題としては,

y=x^2-5

で,y=0となるときのxを求めてみます.つまり$\sqrt{5}$(2.236・・・)を求めるってことですね.
($-\sqrt{5}$については取りあえずおいておきます.)

早速ですが,ニュートン法の第1バージョン

main.cpp
#include <iostream>

void main()
{
    double x;
    double y;

    x = 5.0;//初期値

    double y1;
    double y2;
    double x1;
    double x2;
    double diff_x;
    double error = 0.00001;
    double delta = 0.00001;

    double next_x;

    x1 = x;

    while (1) {

        x2 = x1 + delta;

        //ここがy=x^2-5
        y1 = x1*x1 - 5.0;
        y2 = x2*x2 - 5.0;

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

        if (fabs(y1 - 0.0) < error) {
            break;
        }
        x1 = next_x;
    }

    x = x1;
    y = y1;

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

}

とりあえずで良いならこれで何とか解は求まる.
でも,いろいろ改善したいことはあって,

  • 式を直打ちしているので,ニュートン法の部分が使い回せない.しかも同じ式を二回も打っていて無駄だしバグになりやすい.
  • 今回はたまたま解に収束してくれたから良いけど,何回もループを回して解にたどりつかなかったらエラーを返して欲しい.

とか,挙げだしたら切りが無い.
これを少しずつクールにして行こう.
まずは,解きたい関数を直打ちするのは良くない.直そう.同じ式を何度も打つのはバグの温床になる.

main.cpp
#include <iostream>

int main()
{
    double x;
    double y;

    x = 5.0;

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

    double y1;
    double y2;
    double x1;
    double x2;
    double diff_x;
    double error = 0.00001;
    double delta = 0.00001;
    double next_x;

    x1 = x;

    while (1) {

        x2 = x1 + delta;

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

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

        if (fabs( y1 - 0 ) < error) {
            break;
        }
        x1 = next_x;
    }

    x = x1;
    y = y1;

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

}

よし,これで関数部分に$y=x^2-5$を書き出したから,汎用性は(少し)出てきた.(ちなみにラムダ式を使っているのは私の勉強も兼ねている.いろいろと新しいものを取り入れてみたいのである)
でもこのままだと,while文とかを直打ちしないといけないから,まだまだ汎用性は無い.
ということで,while文のあたりをまるまる関数にしてしまおう.

main.cpp
#include <iostream>

double solve_by_newton_method(double x ) {

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

    double y1;
    double y2;
    double x1;
    double x2;
    double diff_x;
    double error = 0.00001;
    double delta = 0.00001;
    double next_x;

    x1 = x;

    while (1) {

        x2 = x1 + delta;

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

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

        if (fabs(y1 - 0) < error) {
            break;
        }
        x1 = next_x;
    }

    return x1;

}

void main()
{
    double x;
    double y;
    double xini;

    xini = 5.0;//初期値

    x = solve_by_newton_method(xini);
    y = x * x - 5.0;

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

    return 0;
}

うーっん.main文はすっきりしたけど,今度はsolve_by_newton_method関数の中に$y=x^2-5$が直打ちされてしまった.こりゃだめだ.しかもmain文の最後で直打ちが発生してしまっている.
ってことで,なんとかして$y=x^2-5$を外に持ってこなくてはならない.もういちどfuncの定義をmain文に持ってきたい.
と言うことで次のバージョンは以下の通り.

main.cpp
#include <iostream>

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

    double y1;
    double y2;
    double x1;
    double x2;
    double diff_x;
    double error = 0.00001;
    double delta = 0.00001;
    double next_x;

    x1 = x;

    while (1) {

        x2 = x1 + delta;

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

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

        if (fabs(y1 - 0) < error) {
            break;
        }
        x1 = next_x;
    }

    return x1;
};

int main()
{
    double x;
    double y;
    double xini;

    xini = 5.0;//初期値

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

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

    std::cout << "x = " << x << std::endl;
    std::cout << "y = " << y << std::endl;
    return 0;
}

関数だと上手くかけなかったので,solve_by_newton_methodはラムダ式にしてあります.

コメントで教えてもらったんですが,
c++
template<typename Func>
double solve_by_newton_method(Func&& func, double x, double error, double delta){

これでも大丈夫でした.

そしてmain文にラムダ式で$y=x^2-5$を書いてしまい,ラムダ式をそのままsolve_by_newton_methodに渡すことで,ニュートン法の解きたい関数をそのまま渡せるのです.

だから,ラムダ式の形で解きたい方程式を定義すれば,そのままsolve_by_newton_methodへ渡せば解ける!あとは,かゆいところに手を届かせたいので,ニュートン法のパラメータとなるdeltaとerrorも外部から入力できるようにしてしまいましょう.

main.cpp
#include <iostream>

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

    double y1;
    double y2;
    double x1;
    double x2;
    double diff_x;
    double next_x;

    x1 = x;

    while (1) {

        x2 = x1 + delta;

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

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

        if (fabs(y1 - 0) < error) {
            break;
        }
        x1 = next_x;
    }

    return x1;
};

int main()
{
    double x;
    double y;
    double xini;

    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;



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

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

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

    return 0;
}

と,こんな風に,func2を定義して,solve_by_newton_methodへ渡せば解いてくれますし,ニュートン法の収束条件なども,solve_by_newton_methodの引数として与えました.
これで,そこそこ使い回せるニュートン法ができあがりました.

Why do not you register as a user and use Qiita more conveniently?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
Comments
Sign up for free and join this conversation.
If you already have a Qiita account
Why do not you register as a user and use Qiita more conveniently?
You need to log in to use this function. Qiita can be used more conveniently after logging in.
You seem to be reading articles frequently this month. Qiita can be used more conveniently after logging in.
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away