誰もが一度は通る道かもしれませんが，数値計算の例として，ニュートン法を書いてみたいと思います．なぜc++かというと，私がc++でいろいろと数値計算をしているからです．数値計算というとやっぱりFortranと言う気もしますが，c++はそこそこ実行速度も速いし，いろいろライブラリも整ってるし，オブジェクト指向はそこそこ便利と言うこともあって，私はC++は数値計算には結構向いていると思います．あともちろん，自分が所属している組織のルールで言語指定っていうのもあるしね．
それで，一応例題としては，
y=x^2-5
で，y=0となるときのxを求めてみます．つまり$\sqrt{5}$（2.236･･･）を求めるってことですね．
（$-\sqrt{5}$については取りあえずおいておきます．）
早速ですが，ニュートン法の第1バージョン
#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;
}
とりあえずで良いならこれで何とか解は求まる．
でも，いろいろ改善したいことはあって，
- 式を直打ちしているので，ニュートン法の部分が使い回せない．しかも同じ式を二回も打っていて無駄だしバグになりやすい．
- 今回はたまたま解に収束してくれたから良いけど，何回もループを回して解にたどりつかなかったらエラーを返して欲しい．
とか，挙げだしたら切りが無い．
これを少しずつクールにして行こう．
まずは，解きたい関数を直打ちするのは良くない．直そう．同じ式を何度も打つのはバグの温床になる．
#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文のあたりをまるまる関数にしてしまおう．
#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文に持ってきたい．
と言うことで次のバージョンは以下の通り．
#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も外部から入力できるようにしてしまいましょう．
#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の引数として与えました．
これで，そこそこ使い回せるニュートン法ができあがりました．