numeric-methods/lab1.py at main · ark2016/numeric-methods · 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
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
import numpy as np
import matplotlib
import matplotlib.pyplot as plt

matplotlib.use('TkAgg')

def f(x: float) -> float:
    return np.log(x)**2 / x

a = 1/np.e
b = np.e

h = (b-a)/32

def print_func(a, b, h):
    print([[x, f(x)] for x in np.arange(a, b, h)])

print_func(a, b, h)

# Generate points
points = [[x, f(x)] for x in np.arange(a, b, h)]

plt.figure(figsize=(10, 6))
plt.plot([point[0] for point in points], [point[1] for point in points], marker='o')
plt.title('Plot of f(x) = (ln(x)^2) / x')
plt.xlabel('x')
plt.ylabel('f(x)')
plt.grid(True)
plt.show()

def thomas_algorithm(a, b, c, d):
    n = len(d)
    c_prime = np.zeros(n)
    d_prime = np.zeros(n)
    c_prime[0] = c[0] / b[0]
    d_prime[0] = d[0] / b[0]

    for i in range(1, n):
        temp = b[i] - a[i] * c_prime[i - 1]
        c_prime[i] = c[i] / temp
        d_prime[i] = (d[i] - a[i] * d_prime[i - 1]) / temp

    x = np.zeros(n)
    x[-1] = d_prime[-1]

    for i in range(n - 2, -1, -1):
        x[i] = d_prime[i] - c_prime[i] * x[i + 1]

    return x

n = len(points) - 1

y = [point[1] for point in points]

a_coeffs = np.ones(n - 1)
b_coeffs = 4 * np.ones(n - 1)
c_coeffs = np.ones(n - 1)
d_coeffs = np.array([(y[i + 2] - 2 * y[i + 1] + y[i]) / h**2 for i in range(n - 1)])

c_0 = 0
c_n = 0
c_i_values = [c_0] + list(thomas_algorithm(a_coeffs, b_coeffs, c_coeffs, d_coeffs)) + [c_n]

a_i = y

b_i = [(y[i+1] - y[i]) / h - h / 3 * (c_i_values[i+1] + 2 * c_i_values[i]) for i in range(n)]

d_i = [(c_i_values[i+1] - c_i_values[i])/(3*h) for i in range(n)]

print("a_i values:", a_i)
print("b_i values:", b_i)
print("c_i values:", c_i_values)
print("d_i values:", d_i)

x_i = [a + (i - .5) * h for i in range(n + 1)]
points_2 = [[x, f(x)] for x in x_i]

print(points_2)

plt.figure(figsize=(10, 6))
plt.plot([point[0] for point in points_2], [point[1] for point in points_2], marker='o', color='red')
plt.title('Plot of f(x) = (ln(x)^2) / x')
plt.xlabel('x')
plt.ylabel('f(x)')
plt.grid(True)
plt.show()

def evaluate_spline(x, x_i, a_i, b_i, c_i, d_i):
    i = np.searchsorted(x_i, x) - 1
    i = max(0, min(i, len(x_i) - 2))

    dx = x - x_i[i]
    value = a_i[i] + b_i[i] * dx + c_i[i] * dx**2 + d_i[i] * dx**3
    return value

x_eval = 1.5
original_value = f(x_eval)
spline_value = evaluate_spline(x_eval, x_i, a_i, b_i, c_i_values, d_i)

print(f"Original function value at x = {x_eval}: {original_value}")
print(f"Spline value at x = {x_eval}: {spline_value}")


x_eval = a
original_value = f(x_eval)
spline_value = evaluate_spline(x_eval, x_i, a_i, b_i, c_i_values, d_i)

print(f"Original function value at x = {x_eval}: {original_value}")
print(f"Spline value at x = {x_eval}: {spline_value}")

x_eval = b
original_value = f(x_eval)
spline_value = evaluate_spline(x_eval, x_i, a_i, b_i, c_i_values, d_i)

print(f"Original function value at x = {x_eval}: {original_value}")
print(f"Spline value at x = {x_eval}: {spline_value}")