方程式を解くための求根アルゴリズムのひとつにニュートン法があります。
\(x_{n+1} = x_n – \frac{f(x_n)}{f'(x_n)}\) の漸化式を反復的に計算して \(x\) の収束値を算出するアルゴリズムです。
この記事では、C++ でのニュートン法のプログラム例を紹介します。
ニュートン法 プログラム例
早速ですが実装例です。
ニュートン法への入力引数としては次のものが必要になります。
| 引数 | 説明 |
|---|---|
| f | \(f(x)\) の関数 |
| f_prime | \(f'(x)\) の関数 |
| x0 | \(x\) の初期値 |
| torelance | 誤差許容値 |
| max_iter | 最大試行回数 |
漸化式を逐次計算するだけのため、実装はとてもシンプルです。
次のようになります。
C++
#include <cmath>
#include <functional>
double newton_method(
std::function<double(double)> f,
std::function<double(double)> f_prime,
double x0,
double torelance,
std::size_t max_iter) {
double x = x0;
for (std::size_t i = 0; i < max_iter; ++i) {
double fx = f(x);
if (std::abs(fx) < torelance) {
return x;
}
x -= fx / f_prime(x);
}
throw std::runtime_error("最大試行回数に達しました");
}ニュートン法 実行例
以下、\(f(x) = x^2 – 2\) の方程式をニュートン法で計算する例です。
C++
#include <cmath>
#include <functional>
#include <iomanip>
#include <iostream>
double newton_method(
...
}
double f(double x) {
return x * x - 2.0;
}
double f_prime(double x) {
return 2.0 * x;
}
int main() {
auto y = newton_method(f, f_prime, 1.0, 1e-12, 100);
std::cout << std::fixed << std::setprecision(11) << y << std::endl;
// Output: 1.41421356237
return 0;
}
ニュートン法では、\(f(x)\) の他に、\(f'(x)\) を計算するための関数も必要になります。
\(f'(x) = 2x\) の関数をf_primeとして定義します。
実行結果として \(\sqrt{2} = 1.414\cdots\) が得られます。
目次
自動微分によるニュートン法の改良
ニュートン法の改良として、自動微分を活用することで一次微分式 f_prime の定義を不要にすることができます。
以下の記事で解説します。
コメント