1
+ '''
2
+ Author: P Shreyas Shetty
3
+ Implementation of Newton-Raphson method for solving equations of kind
4
+ f(x) = 0. It is an iterative method where solution is found by the expression
5
+ x[n+1] = x[n] + f(x[n])/f'(x[n])
6
+ If no solution exists, then either the solution will not be found when iteration
7
+ limit is reached or the gradient f'(x[n]) approaches zero. In both cases, exception
8
+ is raised. If iteration limit is reached, try increasing maxiter.
9
+ '''
10
+
11
+ import math as m
12
+
13
+ def calc_derivative (f , a , h = 0.001 ):
14
+ '''
15
+ Calculates derivative at point a for function f using finite difference
16
+ method
17
+ '''
18
+ return (f (a + h )- f (a - h ))/ (2 * h )
19
+
20
+ def newton_raphson (f , x0 = 0 , maxiter = 100 , step = 0.0001 , maxerror = 1e-6 ,logsteps = False ):
21
+
22
+ a = x0 #set the initial guess
23
+ steps = [a ]
24
+ error = abs (f (a ))
25
+ f1 = lambda x :calc_derivative (f , x , h = step ) #Derivative of f(x)
26
+ for _ in range (maxiter ):
27
+ if f1 (a ) == 0 :
28
+ raise ValueError ("No converging solution found" )
29
+ a = a - f (a )/ f1 (a ) #Calculate the next estimate
30
+ if logsteps :
31
+ steps .append (a )
32
+ error = abs (f (a ))
33
+ if error < maxerror :
34
+ break
35
+ else :
36
+ raise ValueError ("Itheration limit reached, no converging solution found" )
37
+ if logsteps :
38
+ #If logstep is true, then log intermediate steps
39
+ return a , error , steps
40
+ return a , error
41
+
42
+ if __name__ == '__main__' :
43
+ import matplotlib .pyplot as plt
44
+ f = lambda x :m .tanh (x )** 2 - m .exp (3 * x )
45
+ solution , error , steps = newton_raphson (f , x0 = 10 , maxiter = 1000 , step = 1e-6 , logsteps = True )
46
+ plt .plot ([abs (f (x )) for x in steps ])
47
+ plt .xlabel ("step" )
48
+ plt .ylabel ("error" )
49
+ plt .show ()
50
+ print ("solution = {%f}, error = {%f}" % (solution , error ))
0 commit comments