以前の記事で C++ でのニュートン法プログラム例を紹介しました。
ニュートン法は、方程式を解くための求根アルゴリズムの一種ですが、対象関数の微分式が必要であるという特徴があります。
別の記事で、二重数による自動微分のプログラム例を紹介しました。
これを応用することで、ニュートン法での微分式の定義を不要とすることができます。
この記事では、自動微分とニュートン法の組み合わせ例を紹介します。
目次
自動微分用ヘッダファイルを準備
自動微分については、基本的に以前紹介したものをそのまま利用します。
ニュートン法で利用するため、以下のようにヘッダファイル化のみ対処しております。
dual.h
#ifndef DUAL_H
#define DUAL_H
#include <cmath>
#include <iostream>
struct Dual {
double a;
double b;
Dual() : a(0.0), b(0.0) {}
Dual(double a) : a(a), b(0.0) {}
Dual(double a, double b) : a(a), b(b) {}
Dual operator+() const { return *this; }
Dual operator-() const { return Dual(-a, -b); }
};
inline Dual operator+(Dual lhs, Dual rhs) {
return Dual{lhs.a + rhs.a, lhs.b + rhs.b};
}
inline Dual operator-(Dual lhs, Dual rhs) {
return lhs + (-rhs);
}
inline Dual operator*(Dual lhs, Dual rhs) {
return Dual{lhs.a * rhs.a, lhs.a * rhs.b + lhs.b * rhs.a};
}
inline Dual operator/(Dual lhs, Dual rhs) {
return Dual{lhs.a / rhs.a, (lhs.b * rhs.a - lhs.a * rhs.b) / (rhs.a * rhs.a)};
}
/* 数学関数の定義等は省略 */
#endif // DUAL_H
なお、+= 演算子の定義不足、数学関数の定義不足などあります。
この記事はニュートン法との組み合わせが目的のため、この点は引き続き完備しておりません。
ニュートン法の改良
自動微分を利用して改良したニュートン法です。
C++
#include "dual.h"
#include <functional>
#include <iomanip>
#include <iostream>
double newton_method(
std::function<Dual(Dual)> f,
double x0,
double torelance,
std::size_t max_iter) {
Dual x{x0, 1.0};
for (std::size_t i = 0; i < max_iter; ++i) {
Dual fx = f(x);
if (std::abs(fx.a) < torelance) {
return x.a;
}
x = x - fx.a / fx.b;
}
throw std::runtime_error("最大試行回数に達しました");
}第1引数である関数 f について、関数の入出力方を double から Dual に変更しています。
これによって、一次微分の計算値を計算できるようになることから、微分関数の定義が不要になります。
以前の記事でのニュートン法との比較
比較のため、以前の記事でのニュートン法も掲載しておきます。
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("最大試行回数に達しました");
}プログラム実行例
ここで定義したニュートン法の実行例です。
一次微分式を定義していないにも関わらず、ニュートン法の収束計算ができることが確認できます。
C++
Dual f(Dual x) {
return x * x - 2.0;
}
int main() {
auto y = newton_method(f, 1.0, 1e-12, 100);
std::cout << std::fixed << std::setprecision(11) << y << std::endl;
// Output: 1.41421356237
return 0;
}
コメント