C++でニュートン法を実装して非線形方程式の解を求めたい！

このブログの対象者

・数値計算を１ミリも知らない理系学生

・数値計算に興味がある社会人

・何かでC++を使わないといけない人

今回の趣旨

今回は非線形方程式の解法をC++で実装してみようといったものです。

なんでC++やねん！PythonとかFortranでええやろ！という声が聞こえてきます笑

たしかに最近の科学技術関連の計算はPythonが人気ですよね。
数値計算ならFortranがパフォーマンス上優秀とは聞きます。。。

でも僕の通っている大学ではなぜかC++を使う文化があるんですよねー。
早くシラバスを変更してほしい！！

そんなこんなでC++を使ってきた僕が数値計算のいろはをC++で紹介していきます！！
ちなみに僕は素粒子理論の研究室に入っていて数値計算を使った研究をしてます。

そもそもC++の情報ってCやC#に比べて少ないように思います。
なのでC++を使ってレポート課題をやらなければいけない理系学部生にとって有意義なものになったら最高かなって思います！

ニュートン法の簡単な説明

正式にはNewton-Raphson法といって方程式の解を求める時に便利な数値計算法の一つです。

あくまで微分を含まない形の方程式を求める手法ですので微分方程式の解は求められません。

そのような微分方程式を説きたい人はオイラー法とかルンゲクッタ法で調べるとよかと思います。

この辺もいつかブログにしたいです。(But 時間がない。。。)

解を求めるというのは具体的にどういうことか？

$x^2 + 3x + 2 = 0$

という方程式の解を求めたい場合を考えます。

これは中学数学で簡単に解が求まりますよね。

でも今回はプログラミングで求めてみたい！！という趣旨なので因数分解はしません！！
というか因数分解なんかできないような非線形方程式を扱います！

たとえば

$$sinhx + 3x -2 =0$$

なんかは簡単には解けなさそう。。。ってなりますよね。

ニュートン法はこのようなときも役に立つのですよー。

さっそく始めましょう！

まずは下の図をみてください。

とりあえずこんな感じのグラフを考えてみます。

$f(x)=0$の解がほしい場合、

僕たちが考えるべきなのはオレンジ色の関数とｘ軸の交点の座標なんです。

まず始めに任意の初期値$x_0$を定めます。

その$x_0$の時の関数の接線を考えます。

その接線とx軸の交点を$x_1$とします。

こんな感じで$x_0 \to x_1 \to x_2$で続けていきます。

なんかこれ繰り返したら解が求められそうじゃね？というのがニュートン法です。

でも、これって何回も繰り返しても本当の解って定まらなくね？
結局は近似解しか求まらない？
って思ったそこのあなた！。。。その通りです(><)

数値計算においてぴったり解が求まることなんて滅多にないんですよね笑

なんなら弱点もたくさんあったりします。
有名な弱点を２つ紹介します。。。

極値の近くで接線飛ばされちゃったパターン

これは初期値の値が関数の極値の近くにある場合です。

ご存知の通り極値は傾きが０で、その接線はx軸に平行になります。

なのでx軸に交わらない、または交わっても遠方になってしまうんですよね。

いい感じに極値をまたいでるパターン

こんな感じに接線が無限に繰り返しちゃう場合があるんですよ。

８の字で永遠に繰り返して解が求まりそうになさそうですねー。

これらの問題をニュートン法の初期値問題と言います。

解決策として減速法というのがあります。

接線の方程式

プログラムを書いていく上で重要な接線の方程式はこんな感じです。

$$f(x) = f^\prime (x_0) (x - x_0) + f(x_0)$$

ここから$x = x_1$、$f(x) = 0$として$x_0$の次のx座標である$x_1$を求めます。

$$x_1 = x_0 - \frac{f(x_0)}{ f^\prime (x_0)}$$

これを一般化すると

$$x_{n + 1} = x_n - \frac{f(x_n)}{ f^\prime (x_n)}$$

となりますね。これが前後関係です。

ここまで大丈夫でしょうか？

とりあえずここまでで下準備は完了です。

次からプログラムを書いていきます！！

C++でプログラムを書いてみる

今回は

$$sinhx + 3x -2 =0$$

の解を求めようと思います。

とりあえずの初期値を$x_0 = 1$でやってみます。
（これは完全に適当です。）

無理っぽかったら他の初期値にしてみたり、減速法を使います。

#include <iostream>
#include<math.h>
#include<string>
using namespace std;
int main()
{
    double x;
    double x0;//initial data
    double x1;//next data
    double delta = 0.00001;

    int i = 0;//繰り返した回数

    auto f = [](double x) {
        return sinh(x) + 3*x - 2;//f(x)
    };
    auto g = [](double x) {
        return coshx + 3;//g(x)：微分した関数
    };
    x0 = 1;

    while (i < 100) {//繰り返し回数の上限
        x1 = x0 - f(x0)/g(x0);
        cout << "x=" << x1 << '\n';
        if (fabs(x1 - x0) < delta) {//前後の差分がdelta以下なら終了
            break;
            }
        x0 = x1;//x1を初期値に代入
        i++;//カウント
    }
    return 0;
}

こんな感じのプログラムを書いてみました。

初期値$x_0$からスタートしていき$x_0 \to x_1 \to x_2$といったようにx軸の値が変化していきます。

じゃあ、その終着点はどうなるん？といった話をします。

これはxの数列が$x_0 \to x_1 \to x_2$となるときに数列の前後の差分が0.00001以下になったときに繰り返しを終わるようにプログラムを作っているんですよね。

0.00001のようなところはより小さい数のほうが精度は高く出ます。
このへんは自分の好きに変えても大丈夫です。

でも単精度の有効数字8桁までには納めた方が良さそうです。へい！

これを実行するとコンソール上には以下のような数字が出力されました。

x=0.521206
x=0.494932
x=0.494888
x=0.494888

５回目で求まりましたねー。これはかなりいいです。

最終的に解はx=0.494888って教えてくれました。

妥当性を検証してみる。

$$sinhx + 3x -2 = 0$$

の解として$x = 0.494888$をコンピュータ君が教えてくれたので、これを逆に代入してみて本当に0になるのかを確かめてみます。

実際に代入すると

$$f(0.494888) = 0.0000017 \simeq 2.0 * 10^{-6}$$

となりました。精度的には申し分ない気がします－。やったぜ~！

おわりに

ここまで読んでくれてありがとうございます。

もし、間違えているところがあったり、質問などがあるかたはツイッターのDMなどで気楽に教えてくださいーー

では！お疲れ様です！

[https://twitter.com/izukunKK:embed]

いずのBlog

好きなこと書いていくブログ