Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- import math
- import matplotlib.pyplot as plt
- # ---------------------------
- # Function and derivative
- # ---------------------------
- def f(h):
- return h**3 - 10*h + 5*math.exp(-h/2) - 2
- def df(h):
- return 3*h**2 - 10 - (5/2)*math.exp(-h/2)
- # ---------------------------
- # Bisection Method
- # ---------------------------
- def bisection(a, b, es):
- # fa, fb = f(a), f(b)
- fa=f(a)
- fb=f(b)
- if fa*fb > 0:
- print("Bisection method fails.")
- return None
- # raise ValueError("No root in interval: f(a) and f(b) have same sign")
- iters = []
- xr_old = None
- k = 1
- while True: # infinite loop until we break
- xr = (a + b) / 2
- fr = f(xr)
- ea = abs((xr - xr_old) / xr) * 100 if xr_old is not None else None
- iters.append((k, a, b, xr, fr, ea))
- if fr == 0:
- break
- # interval update
- if fa * fr < 0:
- # b, fb = xr, fr
- b = xr
- fb = fr
- else:
- # a, fa = xr, fr
- a = xr
- fa = fr
- # stopping criteria
- if (ea is not None and ea <= es):
- break
- xr_old = xr
- k += 1
- return iters
- # ---------------------------
- # False Position Method
- # ---------------------------
- def false_position(a, b,es):
- fa, fb = f(a), f(b)
- if fa*fb > 0:
- print("False position method is not applicable because f(xl) and f(xu) have the same sign.")
- return None
- iters = []
- xr_old = None
- k = 1
- while True: # run until break
- xr = b - fb * (b - a) / (fb - fa)
- fr = f(xr)
- ea = abs((xr - xr_old) / xr) * 100 if xr_old is not None else None
- iters.append((k, a, b, xr, fr, ea))
- if fr == 0:
- break
- # interval update
- if fa * fr < 0:
- # b, fb = xr, fr
- b = xr
- fb = fr
- else:
- # a, fa = xr, fr
- a = xr
- fa = fr
- # stopping criteria
- if (ea is not None and ea <= es):
- break
- xr_old = xr
- k += 1
- return iters
- # ---------------------------
- # Newton–Raphson Method
- # ---------------------------
- def newton(x0, es):
- iters = []
- xr_old = x0
- # x = x0
- k = 1
- while True:
- fx = f(xr_old)
- dfx = df(xr_old)
- if dfx == 0:
- print("Derivative zero. Not applicable")
- return None
- xr = xr_old - fx / dfx
- ea = abs((xr - xr_old) / xr) * 100
- iters.append((k, xr_old, fx, dfx, xr, ea))
- # stopping criteria
- if (ea is not None and ea <= es):
- break
- xr_old = xr
- k += 1
- return iters
- # ---------------------------
- # Secant Method
- # ---------------------------
- def secant(x0, x1, es):
- iters = []
- k = 1
- while True:
- f0, f1 = f(x0), f(x1)
- if f1 - f0 == 0:
- print("Division by zero. Not applicable")
- return None
- x2 = x1 - (f1 * (x1 - x0)) / (f1 - f0)
- ea = abs((x2 - x1) / x2) * 100
- iters.append((k, x0, x1, x2, f(x2), ea))
- # stopping criteria
- if (ea is not None and ea <= es) :
- break
- # x0, x1 = x1, x2
- x0 = x1
- x1 = x2
- k += 1
- return iters
- # ---------------------------
- # Run all methods
- # ---------------------------
- es=0.001
- a=0.1
- b=0.4
- bis = bisection(a,b,es)
- fp = false_position(a,b,es)
- h0 = 1.5
- h1 = 2.0
- newt = newton(h0,es)
- sec = secant(h0,h1,es)
- def r5(x): return float(f"{x:.5f}")
- # ---------------------------
- # Print Results
- # ---------------------------
- print("\n--- Final Roots ---")
- print("Bisection root:", r5(bis[-1][3]), "Iterations:", len(bis))
- print("False-Position root:", r5(fp[-1][3]), "Iterations:", len(fp))
- print("Newton–Raphson root:", r5(newt[-1][4]), "Iterations:", len(newt))
- print("Secant root:", r5(sec[-1][3]), "Iterations:", len(sec))
- # ---------------------------
- # Print Tables
- # ---------------------------
- print("\n--- Bisection Table ---")
- print("Iter\t a\t\t b\t\t xr\t\t f(xr)\t\t ea(%)")
- for row in bis:
- print(f"{row[0]:<4}\t{row[1]:.6f}\t{row[2]:.6f}\t{row[3]:.6f}\t{row[4]:.6e}\t{row[5]}")
- print("\n--- False Position Table ---")
- print("Iter\t a\t\t b\t\t xr\t\t f(xr)\t\t ea(%)")
- for row in fp:
- print(f"{row[0]:<4}\t{row[1]:.6f}\t{row[2]:.6f}\t{row[3]:.6f}\t{row[4]:.6e}\t{row[5]}")
- print("\n--- Newton–Raphson Table ---")
- print("Iter\t x_i\t\t f(x_i)\t\t f'(x_i)\t x_{i+1}\t\t ea(%)")
- for row in newt:
- print(f"{row[0]:<4}\t{row[1]:.6f}\t{row[2]:.6e}\t{row[3]:.6f}\t{row[4]:.6f}\t{row[5]}")
- print("\n--- Secant Table ---")
- print("Iter\t x_{i-1}\t x_i\t\t x_{i+1}\t f(x_{i+1})\t ea(%)")
- for row in sec:
- print(f"{row[0]:<4}\t{row[1]:.6f}\t{row[2]:.6f}\t{row[3]:.6f}\t{row[4]:.6e}\t{row[5]}")
- # ---------------------------
- # Plot Convergence
- # ---------------------------
- def extract_err(iters, idx=-1):
- return [row[idx] for row in iters if row[idx] is not None]
- plt.figure()
- plt.semilogy(range(1, len(extract_err(bis,5))+1), extract_err(bis,5), 'o-', label="Bisection")
- plt.semilogy(range(1, len(extract_err(fp,5))+1), extract_err(fp,5), 'o-', label="False Position")
- plt.semilogy(range(1, len(extract_err(newt,5))+1), extract_err(newt,5), 'o-', label="Newton–Raphson")
- plt.semilogy(range(1, len(extract_err(sec,5))+1), extract_err(sec,5), 'o-', label="Secant")
- plt.xlabel("Iteration")
- plt.ylabel("Approx. relative error (%)")
- plt.title("Convergence Comparison")
- plt.legend()
- plt.grid(True, which="both")
- plt.show()
Advertisement
Add Comment
Please, Sign In to add comment
