numericals/bisection.py at master · oilshit/numericals · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
import math
import matplotlib.pyplot as plt

def f(x):
    return math.exp(x) - (5 * (x ** 2))

x = [i * 0.05 for i in range(-50, 50 + 1)]
y = [f(i * 0.05) for i in range(-50, 50 + 1)]

plt.plot(x, y, label="f(x) = e^x - 5x^2")

def bisection(a, b, eps):
    a, b = a, b
    iter = 0

    while (abs(a - b) >= eps):
        c = (a + b) / 2
        fa, fc = f(a), f(c)

        if (fa * fc < 0):
            b = c

        else:
            a = c

        iter += 1

    return c, iter

def bisection_plot(a, b, eps) -> None:
    a, b = a, b
    iter = 0

    plt.xlim(-1 * abs((a - b) / 2), abs(a - b))
    plt.ylim(f(abs(a - b)) * -1, f(abs(a - b)))

    plt.plot([a, a], [0, f(a)], color="black", linestyle="dashed")
    plt.plot([b, b], [0, f(b)], color="black", linestyle="dashed")

    while (abs(a - b) >= eps):
        c = (a + b) / 2
        fa, fc = f(a), f(c)

        plt.plot(
            [c, c], [0, f(c)],
            color="black",
            linestyle="dashed",
            label=("iter " + str(iter))
        )

        if (fa * fc < 0):
            b = c

        else:
            a = c

        iter += 1

    plt.grid()
    plt.xlabel("x")
    plt.ylabel("y = f(x)")

    plt.scatter(c, f(c), color="orange", label="solution")

    plt.show()

    return

a, b = 0, 1
eps = 0.00001
eps2 = 0.000001

solution, iter = bisection(a, b, eps)
print("Using bisection method, solution of e^x - 5x^2 =",
      round(solution, 6),
      "with", iter, "iterations")

bisection_plot(a, b, eps)