Permalink
Browse files

Add Python code by Zahir Jacobs

  • Loading branch information...
Andreas Pauley
Andreas Pauley committed Sep 2, 2010
1 parent fe2c711 commit 334c86378587bf775d8bf849ef5b64eb1a0b8c70
Showing with 158 additions and 0 deletions.
  1. +158 −0 Chapter1/Python/numerical-methods.py
@@ -0,0 +1,158 @@
+# This code was written by Zahir Jacobs:
+# http://zahirj.wordpress.com/2009/04/04/complete-listing-of-python-code-for-selected-root-finding-methods/
+
+import math
+from math import *
+
+def deg_to_rad(x):
+ return x*(math.pi/180.0)
+
+def rad_to_deg(x):
+ return x * (180.0/float(math.pi))
+
+def isdiffsign(a,b):
+ return math.copysign(a,b) != a
+
+def swapPoints(x):
+ s = []
+ s = x
+ s.sort()
+ f = s[1]
+ sn = s[2]
+ t = s[0]
+ s[0] = f
+ s[1] = sn
+ s[2] = t
+ return s
+
+def sub_dict(somedict, somekeys, default=math):
+ return dict([ (k, somedict.get(k, default)) for k in somekeys ])
+
+def evaluate(func,x):
+ safe_list = ['math','acos', 'asin', 'atan', 'atan2', 'ceil', 'cos', 'cosh', 'de grees', 'e', 'exp', 'fabs', 'floor', 'fmod', 'frexp', 'hypot', 'ldexp', 'log', 'log10', 'modf', 'pi', 'pow', 'radians', 'sin', 'sinh', 'sqrt', 'tan', 'tanh'] #use the list to filter the local namespace
+ safe_dict = sub_dict(locals(),safe_list)
+ safe_dict['abs'] = abs
+ safe_dict['x']=x
+ return eval(compile(func,"",'eval'),{"__builtins__":{}},safe_dict)
+
+def MullersMethod(func,a,b,r, maxSteps = 30, tolerance = 0.0001):
+ print_header("Muller's method",func)
+ x = [a,b,r]
+ for loopCount in range(maxSteps):
+ x = swapPoints(x)
+ y = evaluate(func,float(x[0])),evaluate(func,float(x[1])),evaluate(func,float(x[2]))
+ #y = map(func,x)
+ h1 = x[1]-x[0]
+ h2 = x[0]-x[2]
+ lam = h2/h1
+ c = y[0]
+ a = (lam*y[1] - y[0]*((1.0+lam))+y[2])/(lam*h1**2.0*(1+lam))
+ b = (y[1] - y[0] - a*((h1)**2.0))/(h1)
+ if b > 0:
+ root = x[0] - ((2.0*c)/(b+ (b**2 - 4.0*a*c)**0.5))
+ else:
+ root = x[0] - ((2.0*c)/(b- (b**2 - 4.0*a*c)**0.5))
+ print "a = %.5f b = %.5f c = %.5f root = %.5f " % (a,b,c,root)
+ print "Current approximation is %.9f" % root
+ if abs(evaluate(func,float(root))) > x[0]:
+ x = [x[1],x[0],root]
+ else:
+ x = [x[2],x[0],root]
+ x = swapPoints(x)
+ print_end()
+
+def biSection(func, a,b,maxSteps = 30, tolerance = 0.0001):
+ print_header("Bisection Method", func)
+ initial = evaluate(func,float(a))
+ for loopCounter in range(maxSteps):
+ midPoint = a + (b-a)/2.0
+ result = evaluate(func,float(midPoint))
+ print "Accuracy is within %.9f " % (abs(b-a)/2.0)
+ if (result == 0 or ( (abs(b-a)/2.0) > 0)):
+ a = midPoint
+ initial = result
+ else:
+ b = midPoint
+ print_end()
+
+def double(x):
+ return 2*x
+
+def Secant(func,a,b,maxSteps=30, tolerance = 0.0001):
+ print_header("secant method",func)
+ if (abs(evaluate(func,float(a))) < abs(evaluate(func,float(b)))):
+ t = b
+ b = a
+ a = t
+ print "a = %.9f b = %.9f" % (a,b)
+ for loopCount in range(maxSteps):
+ p = b - (evaluate(func,float(a)) * ((a-b)/(evaluate(func,float(a))-evaluate(func,float(b)))))
+ print "Current approximation is %.9f" % p
+ if (abs(evaluate(func,float(p))) < tolerance):
+ print "Root is %.9f (%d iterations)" % (p,loopCount+1)
+ return
+ a = b
+ b = p
+ print "Root find stopped at %.9f" % p
+ print_end()
+
+def getR(x1,x2,x3):
+ return (x1*x3-(x2**2))/(x1-2*x2+x3)
+def isConvergent(f, df,x):
+ return ((f(x)*df(x))/(df(x))**2) < 1
+def regulaFalsi(func, a,b, maxSteps = 30, tolerance = 0.0001):
+ print_header("regula falsi (false position)",func)
+ p = 0.0
+ for loopCount in range(maxSteps):
+ p = b - (evaluate(func,float(b)) * ((a-b)/(evaluate(func,float(a))-evaluate(func,float(b)))))
+ print "Current approximation is %.9f" % p
+ if (math.copysign(evaluate(func,float(a)),evaluate(func,float(b))) != evaluate(func,float(a))):
+ b = p
+ else:
+ a = p
+ if (abs(evaluate(func,float(p))) < tolerance):
+ print "Root is %.9f (%d iterations)" % (p,int(loopCount))
+ return
+ print "Root find cancelled at %.9f" % p
+ print_end()
+
+def print_header(t,f):
+ print "\n"
+ print "-" * 50
+ print "Using %s to solve %s" % (t,f)
+ print "-" * 50
+ print "\n"
+
+def print_end():
+ print "-" * 50
+
+def fixedPoint(func, initialApproximation, maxSteps=30, tolerance=0.0001):
+ print_header("fixed point iteration",func)
+ p = initialApproximation
+ loopCounter = 1
+ for loopCounter in range(maxSteps):
+ oldP = evaluate(func,float(p))
+ print "Current approximation is %.9f" % oldP
+ if (abs(p - oldP) < tolerance):
+ print "Approximate root is %.9f (found in %d steps)" % (oldP,int(loopCounter))
+ return
+ p = oldP
+ print_end()
+
+def newtonRaphson(func, derFunc, initialApproximation, maxSteps = 30, tolerance= 0.0001):
+ print_header("newton/raphson",func + " with " + defFunc + " as derivative")
+ for loopCounter in range(maxSteps):
+ p = initialApproximation - (evaluate(func,float(initialApproximation))/evaluate(derFunc,float(initialApproximation)))
+ if (abs(p - initialApproximation) < tolerance):
+ print "Approximate root is %.9f (found in %i steps)" %(p,loopCounter)
+ break
+ print "Current approximation %.9f" % initialApproximation
+ initialApproximation = p
+ print_end()
+
+if __name__ == '__main__':
+ try:
+ print "in __main__"
+ MullersMethod("math.cos(x) - x",0.0,0.5,1.0)
+ finally:
+ print "Fully done."

0 comments on commit 334c863

Please sign in to comment.