やりたいこと
代数方程式f(x)=0をなるべく汎用的に解きたい。これの数値計算で最も有名なのはnewton法だろう。
newton法の実装は巷には溢れかえっている。newton法は微分を使うが筋の悪い書き方をしている例だとdouble func(double)、double df(double)をつまり関数自体と導関数を別々に手で実装している。
少し頑張れば微分は数値微分にできる。更に頑張れば、関数自体を関数ポインタか関数オブジェクトにしてしまえるがそれでも可読性は非常に低い。
それをc++11のstd::functionが劇的に改善したので実装してみる。
Differentaitor(微分器)
ニュートン法で解くために数値微分を定義する。数値微分は簡単なもので前方差分、後方差分、中心差分の3種類がありそれぞれ誤差が$dx, dx, dx^2 $である。当然刻み幅$dx$を小さくすれば精度は上がるがあまり小さくし過ぎるとdouble型の丸め誤差が効いてくる。double(倍精度浮動少数)型で15桁あるがこのレポートによると前方と後方で8桁、中心で5桁で誤差が下がらなくなっているのでその値をデフォルト値として固定しておく。
また、ここでは使わないが複素数用の微分も実装しておく。
#pragma once
#include "common.hh"
namespace MyMath
{
namespace Differentaitor{
double central(const std::function<double(double)> &f, double val, double delta=1e-5);
double forward(const std::function<double(double)> &f, double val, double delta=1e-8);
double backward(const std::function<double(double)> &f, double val, double delta=1e-8);
using complex=std::complex<double>;
complex central(const std::function<complex(complex)> &f, const complex &val, const complex &delta=1e-5);
complex forward(const std::function<complex(complex)> &f, const complex &val, const complex &delta=1e-8);
complex backward(const std::function<complex(complex)> &f, const complex &val, const complex &delta=1e-8);
};
};
実装は非常に単純である。
#include "Differentaitor.hh"
using namespace MyMath;
double Differentaitor::central(const std::function<double(double)> &f, double val, double delta){
return (f(val+delta)-f(val-delta))/(2.0*delta);
}
double Differentaitor::forward(const std::function<double(double)> &f, double val, double delta){
return (f(val+delta)-f(val))/delta;
}
double Differentaitor::backward(const std::function<double(double)> &f, double val, double delta){
return (f(val)-f(val-delta))/delta;
}
using complex=std::complex<double>;
complex Differentaitor::central(const std::function<complex(complex)> &f, const complex &val, const complex &delta){
return (f(val+delta)-f(val-delta))/(2.0*delta);
}
complex Differentaitor::forward(const std::function<complex(complex)> &f, const complex &val, const complex &delta){
return (f(val+delta)-f(val))/delta;
}
complex Differentaitor::backward(const std::function<complex(complex)> &f, const complex &val, const complex &delta){
return (f(val)-f(val-delta))/delta;
}
汎用関数としてテンプレート化も考えたが一般的な関数は$f(x)$は引数と返り値が一致するとも限らず(むしろ一致している方が稀かもしれない)。int型などの離散値は微分自体意味がないので実数と複素数の場合のみ型を明示して実装した。
Solver
関数$f(x)$を引数にとって$f(x)=0$を解くものnewton法と2分値検索を実装した。
収束はnewton法が圧倒的に早いが関数系が複雑なときは初期値によって収束せず振動する場合がある（これを大域的安定性がない(低い)という）ので2分値探索も実装しておく、現在のCPUならばよほど難しい方程式か精度を要求しない限り実行時間の差は体感できないと思う。(そんな状況ならnewton法でも不足でもっと良い実装をするべき)
打ち切り誤差はconst double DEFUALT_ERROR=1.0e-8で定義している。倍精度だと数値微分、積分はこのあたりの値が限界らしい。多分stl組み込み関数もこの程度の誤差を持っているはず。
#pragma once
#include "common.hh"
#include "Differentaitor.hh"
namespace MyMath
{
namespace Solver{
double bynarySearch(const std::function<double(double)> &func, double min, double max, double error=DEFUALT_ERROR);
double newton(const std::function<double(double)> &func, double init, double error=DEFUALT_ERROR);
};
};
実装、newton法のほうが短くなった。
#include "Solver.hh"
using namespace std;
using namespace MyMath;
double Solver::bynarySearch(const std::function<double(double)> &f, double min, double max, double error){
if( min>max ) swap(min, max);
if( f(min)*f(max)>0 ){
cout<<"!!! bynarySearch bad initial parameter ["<<min<<":"<<max<<" !!!"<<endl;
return 0.f/0.f;
}
double mid=0.5*(min+max);
double val_min=f(min);
double val_mid=f(mid);
double val_max=f(max);
while( max-min>error || fabs(val_min)>error ){
if( val_mid*val_max>0 ) max=mid, val_max=val_mid;
else min=mid, val_min=val_mid;
mid=0.5*(max+min);
val_mid=f(mid);
}
return mid;
}
double Solver::newton(const std::function<double(double)> &f, double init, double error){
double next=init-f(init)/Differentaitor::central(f, init);
while( fabs(init-next)>error ){
init=next;
next=init-f(init)/Differentaitor::central(f, init);
}
return next;
}
main関数
#include "Solver.hh"
using namespace std;
using namespace MyMath;
int main(){
cout<<"===== Solver test program ====="<<endl;
cout<<" value : std library : bynary search "<<endl;
for( double root=2; root<100; root+=1.0 ){
auto func=[root](double x){ return x*x-root; };
cout<<"sqrt("<<setw(2)<<root<<") : "<<setw(8)<<sqrt(root)<<" : "
<<setw(8)<<Solver::bynarySearch(func, 0, root)<<" : "
<<setw(8)<<Solver::newton(func, root)<<endl;
}
return 0;
}
数値計算で必ず一回はやらされるアレ、平方根の探索。stlライブラリの値とnewton法と2分値探索の結果を列挙。とりあえず結果に違いは見られない。
実装自体は単純なので最低限にすれば十数行程度で済むと思われる。自分も必要になった時ローカルにベタ書きしていたが間違えても致命的なエラーにならないのでバグらせたりしていた。関数化して簡単なテストを通したものを作っておいた。
#include "common.hh"
#include "Differentaitor.hh"
using namespace std;
using namespace MyMath;
int main(){
cout<<"===== Solver test program ====="<<endl;
cout<<"============ Check real value ========================================================================="<<endl;
cout<<" x : cos(x) : forward : backward : central "<<endl;
auto func=[](double x){ return sin(x); };
for( double x=0; x<PI; x+=0.01 ){
cout<<setw(4)<<x<<" : "<<setw(12)<<cos(x)<<" : "<<setw(12)<<Differentaitor::forward(func, x)<<" : "
<<setw(12)<<Differentaitor::backward(func, x)<<" : "<<setw(12)<<Differentaitor::central(func, x)<<endl;
}
cout<<"============ Check complex value ========================================================================="<<endl;
cout<<" z : cos(z) : forward : backward : central "<<endl;
auto func2=[](complex<double> x){ return sin(x); };
for( complex<double> z=0; norm(z)<PI; z+=0.01*complex<double>(1, 1) ){
cout<<setw(11)<<z<<" : "<<setw(25)<<cos(z)<<" : "<<setw(25)<<Differentaitor::forward(func2, z)<<" : "
<<setw(25)<<Differentaitor::backward(func2, z)<<" : "<<setw(25)<<Differentaitor::central(func2, z)<<endl;
}
return 0;
}
こちらは微分のチェック用、$sin(x)$の微分と$cos(x)$の値を比較している複素微分も問題なかった。
##最後に
ラムダ式とauto、functionの組み合わせは強力、ファンクタを定義する必要もないし関数の引数を標準ライブラリにできるので気が楽
、引数に型も明示できるので可読性もかなり高くなると思う。
自分用ライブラリにするつもりなので全体をnamespace MyMath;で包んでいる。ただnewton法を使うならかなり過多な実装をしている。