Create lab_5

Р3265-69/fedorova_367582/lab5/lab_5

@@ -0,0 +1,310 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import sys
+import math
+
+# Проверка уникальности узлов
+
+def check_unique(x):
+    if len(x) != len(set(x)):
+        raise ValueError("Координаты узлов должны быть уникальными!")
+
+# Ввод данных режим а: с клавиатуры
+
+def input_keyboard():
+    n = int(input("Введите число узлов: "))
+    if n < 2:
+        raise ValueError("Число узлов должно быть не менее двух")
+    x = []
+    y = []
+    for i in range(n):
+        xi = float(input(f"x[{i}] = "))
+        yi = float(input(f"y[{i}] = "))
+        x.append(xi)
+        y.append(yi)
+    x = np.array(x)
+    y = np.array(y)
+    check_unique(x)
+    return x, y
+
+# Ввод данных режим b: из файла
+
+def input_file(filename):
+    data = np.loadtxt(filename)
+    x = data[:,0]
+    y = data[:,1]
+    check_unique(x)
+    return x, y
+
+# Режим c: выбранная пользователем функция
+functions = {
+    '1': np.sin,
+    '2': np.cos,
+    '3': np.exp,
+    '4': lambda x: x**2
+}
+
+def input_function():
+    print("Выберите функцию:")
+    print("1: sin(x)")
+    print("2: cos(x)")
+    print("3: exp(x)")
+    print("4: x^2")
+    key = input("Ваш выбор: ")
+    f = functions.get(key)
+    if f is None:
+        raise ValueError("Неверный выбор функции")
+    a = float(input("Левая граница интервала a = "))
+    b = float(input("Правая граница интервала b = "))
+    n = int(input("Число узлов n >= 2: "))
+    if n < 2:
+        raise ValueError("Число узлов должно быть не менее двух")
+    x = np.linspace(a, b, n)
+    y = f(x)
+    check_unique(x)
+    return x, y, f
+
+# Построение таблицы конечных разностей
+
+def finite_diff_table(x, y):
+    n = len(x)
+    table = np.zeros((n, n))
+    table[:,0] = y
+    for j in range(1, n):
+        for i in range(n - j):
+            table[i,j] = table[i+1,j-1] - table[i,j-1]
+    return table
+
+# Многочлен Лагранжа
+
+def lagrange(xi, x, y):
+    n = len(x)
+    total = 0.0
+    for i in range(n):
+        term = y[i]
+        for j in range(n):
+            if j != i:
+                term *= (xi - x[j]) / (x[i] - x[j])
+        total += term
+    return total
+
+# Многочлен Ньютона с разделенными разностями
+
+def newton_divided(xi, x, y):
+    n = len(x)
+    dd = np.zeros((n, n))
+    dd[:,0] = y
+    for j in range(1, n):
+        for i in range(n - j):
+            dd[i,j] = (dd[i+1,j-1] - dd[i,j-1]) / (x[i+j] - x[i])
+    result = dd[0,0]
+    for j in range(1, n):
+        prod = 1.0
+        for k in range(j):
+            prod *= (xi - x[k])
+        result += dd[0,j] * prod
+    return result
+
+# Многочлен Ньютона с конечными разностями
+
+def newton_finite(xi, x, y):
+    n = len(x)
+    h = x[1] - x[0]
+    table = finite_diff_table(x, y)
+    result = y[0]
+    for i in range(1, n):
+        coeff = table[0,i] / math.factorial(i)
+        t = (xi - x[0]) / h
+        term = 1
+        for k in range(i):
+            term *= (t - k)
+        result += coeff * term
+    return result
+
+# Схема Стирлинга (центральная)
+
+def stirling(xi, x, y):
+    if not is_uniform(x):
+        raise ValueError("Узлы должны быть равноотстоящими")
+    n = len(x)
+    if n % 2 == 0: return None
+    h = x[1] - x[0]
+    # Находим ближайший центральный узел
+    m = np.argmin(np.abs(x - xi))
+    t = (xi - x[m]) / h
+    table = finite_diff_table(x, y)
+    # Стартовое значение среднее двух центральных
+    result = (y[m] + y[m-1]) / 2
+    for i in range(1, min(m, table.shape[1] // 2 + 1)):
+        if (2 * i - 1) >= table.shape[1]:
+            break
+        if m - i - 1 < 0:
+            break
+        delta = (table[m-i, 2*i-1] + table[m-i-1, 2*i-1]) / 2
+        term = delta
+        for k in range(1, 2*i):
+            term *= (t**2 - k**2)
+        result += term / math.factorial(2*i)
+    return result
+
+# Схема Бесселя (центральная)
+
+def bessel(xi, x, y):
+    if not is_uniform(x):
+        raise ValueError("Узлы должны быть равноотстоящими")
+    n = len(x)
+    if n % 2: return None
+    h = x[1] - x[0]
+    n = len(x)
+
+    # Найдём индекс m такой, что xi лежит ближе к середине x[m] и x[m+1]
+    m = np.searchsorted(x, xi) - 1
+    m = max(1, min(m, n - 3))  # чтобы не выйти за границы при обращении к table
+
+    x_center = (x[m] + x[m + 1]) / 2
+    t = (xi - x_center) / h
+
+    table = finite_diff_table(x, y)
+
+    # Начальное приближение — среднее значение y[m] и y[m+1]
+    result = (y[m] + y[m + 1]) / 2
+
+    # Первая разность
+    result += t * table[m, 1]
+
+    # Цикл по i — добавление чётных и нечётных членов
+    for i in range(1, n // 2):
+        if m - i < 0 or m - i - 1 < 0 or m + i + 1 >= n:
+            break  # защита от выхода за границы
+
+        # чётная разность (средняя из Δ^{2i})
+        delta_even = (table[m - i, 2 * i] + table[m - i - 1, 2 * i]) / 2
+        prod_even = 1
+        for k in range(1, i + 1):
+            prod_even *= (t ** 2 - (k - 0.5) ** 2)
+        result += delta_even * prod_even / math.factorial(2 * i)
+
+        # нечётная разность (Δ^{2i+1})
+        delta_odd = table[m - i, 2 * i + 1]
+        prod_odd = t
+        for k in range(1, i + 1):
+            prod_odd *= (t ** 2 - k ** 2)
+        result += delta_odd * prod_odd / math.factorial(2 * i + 1)
+
+    return result
+
+# функция проверки равномерности узлов
+def is_uniform(x):
+    if len(x) < 2:
+        return True
+    h = x[1] - x[0]
+    return np.allclose(np.diff(x), h)
+
+# Построение графиков
+
+def plot_interpolation(x, y, f_exact=None, x0=None):
+    xs = np.linspace(min(x), max(x), 300)
+    plt.figure(figsize=(8, 6))
+
+    if f_exact is not None:
+        plt.plot(xs, f_exact(xs), label='Исходная функция')
+
+    # Сортируем узлы для корректного отображения
+    sorted_idx = np.argsort(x)
+    x_sorted = x[sorted_idx]
+    y_sorted = y[sorted_idx]
+
+    # Строим интерполяционные кривые
+    ys_newton = [newton_divided(xx, x_sorted, y_sorted) for xx in xs]
+    plt.plot(xs, ys_newton, linestyle='--', label='Полином Ньютона')
+
+    ys_lagr = [lagrange(xx, x_sorted, y_sorted) for xx in xs]
+    plt.plot(xs, ys_lagr, linestyle=':', label='Полином Лагранжа')
+
+    # Добавим узлы интерполяции
+    plt.scatter(x_sorted, y_sorted, color='red', label='Узлы интерполяции')
+
+    # Добавим точку x0 и значения в ней
+    if x0 is not None:
+        y_newton = newton_divided(x0, x_sorted, y_sorted)
+        y_lagr = lagrange(x0, x_sorted, y_sorted)
+
+        plt.scatter(x0, y_newton, color='blue', marker='x', s=100, label='Newton в x0')
+        plt.scatter(x0, y_lagr, color='green', marker='^', s=100, label='Lagrange в x0')
+
+        if f_exact is not None:
+            y_real = f_exact(x0)
+            plt.scatter(x0, y_real, color='black', marker='o', s=100, label='f(x0)')
+
+    plt.legend()
+    plt.xlabel('x')
+    plt.ylabel('y')
+    plt.title('Интерполяция')
+    plt.grid(True)
+    plt.show()
+
+
+# Основная функция
+
+def main():
+    print("Режимы ввода данных:")
+    print("a: с клавиатуры")
+    print("b: из файла")
+    print("c: по выбранной функции")
+    mode = input("Выберите режим (a/b/c): ")
+    try:
+        if mode == 'a':
+            x, y = input_keyboard()
+            f = lambda z: None
+        elif mode == 'b':
+            filename = input("Имя файла с данными: ")
+            x, y = input_file(filename)
+            f = lambda z: None
+        elif mode == 'c':
+            x, y, f = input_function()
+        else:
+            raise ValueError("Неверный режим")
+    except Exception as e:
+        print("Ошибка ввода данных:", e)
+        sys.exit(1)
+
+    # Сортируем данные по x
+    sort_idx = np.argsort(x)
+    x = x[sort_idx]
+    y = y[sort_idx]
+
+    # Таблица конечных разностей
+    diff_table = finite_diff_table(x, y)
+    print("\nТаблица конечных разностей:")
+    print(diff_table)
+
+    # Запрос аргумента для интерполяции
+    xi = float(input("\nВведите значение аргумента xi для расчета: "))
+    print("\nМетод Лагранжа:", lagrange(xi, x, y))
+    print("Метод Ньютона (разделенные разности):", newton_divided(xi, x, y))
+    print("Метод Ньютона (конечные разности):", newton_finite(xi, x, y))
+
+    # Проверка и вывод для Стирлинга
+    stirling_res = stirling(xi, x, y)
+    if stirling_res is not None:
+        print("Схема Стирлинга:", stirling_res)
+    else:
+        print("Схема Стирлинга: не применима (требуются равноотстоящие узлы и нечетное n)")
+
+    # Проверка и вывод для Бесселя
+    bessel_res = bessel(xi, x, y)
+    if bessel_res is not None:
+        print("Схема Бесселя:", bessel_res)
+    else:
+        print("Схема Бесселя: не применима (требуются равноотстоящие узлы и четное n)")
+
+    # Вызываем построение графиков всегда
+    if mode == 'c':
+        plot_interpolation(x, y, f, xi)
+    else:
+        plot_interpolation(x, y, None, xi)
+
+
+if __name__ == "__main__":
+    main()
+