From 02f850965d4f36c58b5cb52fc8aba41454cd2f12 Mon Sep 17 00:00:00 2001 From: P-Shreyas-Shetty Date: Mon, 11 Feb 2019 21:45:49 +0530 Subject: [PATCH] Implementation of Newton-Raphson method (#650) Implemented Newton-Raphson method using pure python. Third party library is used only for visualizing error variation with each iteration. --- maths/newton_raphson.py | 50 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 50 insertions(+) create mode 100644 maths/newton_raphson.py diff --git a/maths/newton_raphson.py b/maths/newton_raphson.py new file mode 100644 index 000000000..c08bcedc9 --- /dev/null +++ b/maths/newton_raphson.py @@ -0,0 +1,50 @@ +''' + Author: P Shreyas Shetty + Implementation of Newton-Raphson method for solving equations of kind + f(x) = 0. It is an iterative method where solution is found by the expression + x[n+1] = x[n] + f(x[n])/f'(x[n]) + If no solution exists, then either the solution will not be found when iteration + limit is reached or the gradient f'(x[n]) approaches zero. In both cases, exception + is raised. If iteration limit is reached, try increasing maxiter. + ''' + +import math as m + +def calc_derivative(f, a, h=0.001): + ''' + Calculates derivative at point a for function f using finite difference + method + ''' + return (f(a+h)-f(a-h))/(2*h) + +def newton_raphson(f, x0=0, maxiter=100, step=0.0001, maxerror=1e-6,logsteps=False): + + a = x0 #set the initial guess + steps = [a] + error = abs(f(a)) + f1 = lambda x:calc_derivative(f, x, h=step) #Derivative of f(x) + for _ in range(maxiter): + if f1(a) == 0: + raise ValueError("No converging solution found") + a = a - f(a)/f1(a) #Calculate the next estimate + if logsteps: + steps.append(a) + error = abs(f(a)) + if error < maxerror: + break + else: + raise ValueError("Itheration limit reached, no converging solution found") + if logsteps: + #If logstep is true, then log intermediate steps + return a, error, steps + return a, error + +if __name__ == '__main__': + import matplotlib.pyplot as plt + f = lambda x:m.tanh(x)**2-m.exp(3*x) + solution, error, steps = newton_raphson(f, x0=10, maxiter=1000, step=1e-6, logsteps=True) + plt.plot([abs(f(x)) for x in steps]) + plt.xlabel("step") + plt.ylabel("error") + plt.show() + print("solution = {%f}, error = {%f}" % (solution, error)) \ No newline at end of file