GopherProxy

	tnewton.cpp - numeric - C++ library with numerical algorithms
	git clone git://src.adamsgaard.dk/numeric
	Log
	Files
	Refs
	LICENSE
	---
	tnewton.cpp (1761B)
	---
	1 #include <functional>
	2 #include <armadillo>
	3 #include "qrfunc.h"
	4
	5 using namespace arma;
	6
	7 // Newton's method
	8 vec newton(std::function<vec(vec)> f, vec x_0, vec dx, double eps)
	9 {
	10 vec x = x_0;
	11 int n = x.size();
	12 mat A(n,n);
	13 vec y(n);
	14 vec fx(n);
	15 vec fy(n);
	16 vec df(n);
	17 vec Dx(n);
	18
	19 do {
	20 fx = f(x);
	21 for (int j=0; j<n; ++j) {
	22 x[j] += dx[j];
	23 df = f(x) - fx;
	24
	25 for (int i=0; i<n; ++i)
	26 A(i,j) = df[i]/dx[j];
	27
	28 x[j] -= dx[j];
	29 }
	30
	31 // Perform QR decomposition
	32 QR qr(A);
	33 vec fx_neg = -1.0f*fx;
	34 Dx = qr.backsub(fx_neg);
	35
	36 double lambda = 2.0f;
	37 do {
	38 lambda /= 2.0f;
	39 y = x + Dx * lambda;
	40 fy = f(y);
	41 } while (norm(fy,2.0f) > (1.0f-lambda/2.0f)*norm(fx,2.0f) && lambda …
	42
	43 x = y;
	44 fx = fy;
	45 std::cerr << "Newton: f(x).norm() = " << norm(fx, 2.0f) << '\n';
	46 } while (norm(Dx,2.0f) > norm(dx,2.0f) && norm(fx,2.0f) > eps);
	47
	48 // Return solution
	49 return x;
	50 }
	51
	52 // Newton's method - the user supplies his own derivatives
	53 vec newtonJac(std::function<vec(vec)> f, vec x_0, vec dx, double eps,
	54 mat (*J)(vec))
	55 {
	56 vec x = x_0;
	57 int n = x.size();
	58 vec y(n);
	59 vec fx(n);
	60 vec fy(n);
	61 vec Dx(n);
	62
	63 fx = f(x);
	64
	65 do {
	66
	67 // Get Jacobian matrix
	68 mat Jx = J(x_0);
	69
	70 // Perform QR decomposition
	71 QR qr(Jx);
	72 vec fx_neg = -1.0f*fx;
	73 Dx = qr.backsub(fx_neg);
	74
	75 double lambda = 2.0f;
	76 do {
	77 lambda /= 2.0f;
	78 y = x + Dx * lambda;
	79 fy = f(y);
	80 } while (norm(fy,2.0f) > (1.0f-lambda/2.0f)*norm(fx,2.0f) && lambda …
	81
	82 x = y;
	83 fx = fy;
	84 std::cerr << "Newton: f(x).norm() = " << norm(fx, 2.0f) << '\n';
	85 } while (norm(Dx,2.0f) > norm(dx,2.0f) && norm(fx,2.0f) > eps);
	86
	87 // Return solution
	88 return x;
	89 }