From bf1e581ce9cf21776666e9c968eaa96a07eab930 Mon Sep 17 00:00:00 2001
From: Your Name <guilhem.lavaux@gmail.com>
Date: Wed, 9 Feb 2011 14:31:45 -0500
Subject: [PATCH] Added newton's method

---
 sample/testNewton.cpp | 27 +++++++++++++++++++++++++++
 src/newton.hpp        | 28 ++++++++++++++++++++++++++++
 2 files changed, 55 insertions(+)
 create mode 100644 sample/testNewton.cpp
 create mode 100644 src/newton.hpp

diff --git a/sample/testNewton.cpp b/sample/testNewton.cpp
new file mode 100644
index 0000000..27415b8
--- /dev/null
+++ b/sample/testNewton.cpp
@@ -0,0 +1,27 @@
+#include <iostream>
+#include <cmath>
+
+using namespace std;
+
+#include "newton.hpp"
+
+using namespace CosmoTool;
+
+struct SimplePolynom
+{
+  double eval(double x) { return x*x*x+2*x-1; }
+  double derivative(double x) { return 3*x*x+2; }
+};
+
+int main()
+{
+  SimplePolynom poly;
+  double r;
+
+  double solution = -2*pow(2./(3*(9+sqrt(177.))), 1./3) + (pow(0.5*(9+sqrt(177.))/9, 1./3));
+
+  r = newtonSolver(3.0, poly);
+
+  cout << "Result = " << r << " delta = " << r-solution << endl;
+  return 0;
+}
diff --git a/src/newton.hpp b/src/newton.hpp
new file mode 100644
index 0000000..ca39804
--- /dev/null
+++ b/src/newton.hpp
@@ -0,0 +1,28 @@
+#ifndef _COSMOTOOL_NEWTON_HPP
+#define _COSMOTOOL_NEWTON_HPP
+
+namespace CosmoTool
+{
+  template<typename T, typename FunT>
+  T newtonSolver(T x0, FunT function, double residual = 1e-3)
+  {
+    T x, xold = x0;
+    T f_x = function.eval(x0);
+    T df_x = function.derivative(x0);
+    
+    x = xold - f_x/df_x;
+
+    while (abs(xold-x) > residual)
+      {	
+	xold = x;
+	f_x = function.eval(x);
+	df_x = function.derivative(x);
+	x = xold - f_x/df_x;
+      }
+	   
+    return x;
+  }
+
+};
+
+#endif