Different method to find roots of a function

Author: only-dev-ops

aitkens_method.py

@@ -0,0 +1,19 @@
+import math
+e = math.e
+pi = math.pi
+
+accelerate = lambda x0, x1, x2: (x0*x2 - x1*x1)/(x0 + x2 - 2*x1)
+
+def aitkens(f, p, n):
+    initial_values = [0] * (n+2)
+    for i in range(n+2):
+        initial_values[i] = f(p)
+        p = initial_values[i]
+    print(initial_values)
+    for i in range(n):
+        result = accelerate(initial_values[i], initial_values[i+1], initial_values[i+2])
+        print(f'P{i+1}:\t{result:.7f}')
+
+if __name__ == '__main__':
+    f = lambda x: math.cos(x)
+    aitkens(f, 0.5, 5)

bisection_method.py

@@ -0,0 +1,26 @@
+import math
+
+def bisect(f,a,b, TOL):
+    if (f(a)*f(b) >= 0):
+        print("You have not assumed the right interval.")
+        return
+    mid = a
+    iterations = 0
+    while ((b-a) >= TOL):
+        mid = (a+b)/2
+        result = f(mid)
+        # Check if middle point is root
+        if (result == 0):
+            break
+        # Decide the side to repeat the steps
+        if (result*f(a) < 0):
+            b = mid
+        else:
+            a = mid
+        iterations += 1
+        print(f'P{iterations}: {mid:.5f}\tResult: {result}')
+    print(f'Root: {mid:.7f}\nIterations: {iterations}')
+
+if __name__ == '__main__':
+    f = lambda x: 0.331 - math.asin(x) - (x*(1-x*x)**0.5)
+    bisect(f, 0, 1, 0.01)

fixed_point.py

@@ -0,0 +1,16 @@
+import math
+
+def fixed_point_iteration(f, p, TOL):
+    error = 1
+    iterations = 0
+    while error > TOL:
+        p_new = f(p)
+        error = abs(p_new - p)
+        p = p_new
+        iterations += 1
+        print(f'p{iterations} = {p: 0.5f}')
+    print(f'Root: {p}\nIterations: {iterations}')
+
+if __name__ == '__main__':
+    f = lambda x: (math.sin(x) + math.cos(x))/2
+    fixed_point_iteration(f, 0, 1e-5)

muellers_method.py

@@ -0,0 +1,38 @@
+MAX_ITERATIONS = 100
+
+def muellers(f, p0, p1, p2, TOL):
+    iterations = 3
+    error = 1
+    print(f'Iteration\t\tPi')
+    while error > TOL and iterations < MAX_ITERATIONS:
+
+        h1 = p1-p0
+        h2 = p2-p1
+        del1 = (f(p1)-f(p0))/h1
+        del2 = (f(p2)-f(p1))/h2
+        d = (del2-del1)/(h2+h1)
+        c = f(p2)
+
+        b = del2 + h2*d
+        discriminant = (b*b - 4*c*d)**0.5
+
+        E = b+discriminant if abs(b-discriminant) < abs(b+discriminant) else b-discriminant
+
+        proximity = -2*c/E
+        p = p2+proximity
+
+        print(f'{iterations:^5}\t\t{p:.5f}')
+
+        error = abs(proximity)
+
+        p0 = p1
+        p1 = p2
+        p2 = p
+        iterations += 1
+
+    print("Incorrect initial points. Try again." if iterations > MAX_ITERATIONS else "Exiting...")
+
+if __name__ == '__main__':
+    #f = lambda x: x**4 + 2*x**2 - x - 3
+    f = lambda x: x**5 - x**4 + 2*x**3 - 3*x**2 + x - 4
+    muellers(f, 1, 0, -1, 1e-5)

newton_raphson.py

@@ -0,0 +1,41 @@
+import math
+e = math.e
+pi = math.pi
+
+def newton_raphson(f, df, x, TOL):
+    error = 1
+    iterations = 0
+    while error > TOL:
+        new_x = x - f(x)/df(x)
+        error = abs(new_x - x)
+        x = new_x
+        iterations += 1
+        print(f'x{iterations}: {x}')
+    print(f"Newton's Estimate = {x:.15f}\nIterations: {iterations}")
+
+def modified_newton(f, df, ddf, x, TOL):
+    error = 1
+    iterations = 0
+    while error > TOL:
+        f_x = f(x)
+        d_x = df(x)
+        new_x = x - (f_x*d_x)/(d_x*d_x - f_x*ddf(x))
+        error = abs(new_x - x)
+        x = new_x
+        iterations += 1
+        print(f'x{iterations}: {x}')
+    print(f"Modified Newton Estimate = {x:.15f}\nIterations: {iterations}")
+
+def fibonacci_estimate():
+    values = [1, 22, 7, 42, 33, 4, 40]
+    total = 0
+    for i in range(len(values)):
+        total += values[i]/(60**i)
+    print(f"Fibonacci's Estimate = {total:.15f}")
+
+if __name__ == '__main__':
+    f = lambda x: x**4 - 2*x**3 - 12*x*x + 16*x - 40
+    df = lambda x: 4*x**3 - 6*x*x - 24*x + 16
+    ddf = lambda x: 12*x*x - 12*x - 24
+    newton_raphson(f, df, 1, 1e-5)
+    modified_newton(f, df, ddf, 1, 1e-5)

secant_method.py

@@ -0,0 +1,18 @@
+import math
+
+def secant(f, x1, x2, TOL):
+    error = 1
+    iterations = 0
+    print(f'x0: {x1}\nx1: {x2}')
+    while error > TOL:
+        x_new = x2 - (f(x2)*(x2-x1))/(f(x2)-f(x1))
+        error = abs(x_new - x2)
+        x1 = x2
+        x2 = x_new
+        iterations += 1
+        print(f'x{iterations+1}: {x2}')
+    print(f'Result: {x2}\nIterations: {iterations}')
+
+if __name__ == '__main__':
+    f = lambda x: x - math.cos(x)
+    secant(f, 0, math.pi/2, 1e-4)

steffensens_method.py

@@ -0,0 +1,20 @@
+import math
+
+accelerate = lambda x0, x1, x2: (x0*x2 - x1*x1)/(x0 + x2 - 2*x1)
+
+def steffensens(f, p, TOL):
+    error = 1
+    iterations = 0
+    while error > TOL:
+        p1 = f(p)
+        p2 = f(p1)
+        p_new = accelerate(p, p1, p2)
+        error = abs(p_new-p)
+        iterations += 1
+        print(f'P0: {p:.7f}\nP1: {p1:.7f}\nP2: {p2:.7f}\nP-hat_{iterations}: {p_new:.7f}\n')
+        p = p_new
+    print(f'Final value: {p}\nIterations: {iterations}')
+
+if __name__ == '__main__':
+    f = lambda x: 5**-x
+    steffensens(f, 0, 5e-2)