From f4860f706fcfeb0b4cf983da3a1bfc7078d908c3 Mon Sep 17 00:00:00 2001 From: Alina1344 <367582@edu.itmo.ru> Date: Fri, 30 May 2025 15:43:44 +0300 Subject: [PATCH] Create lab_6 MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Добавлена лаба6 Fedorova_367582 --- "\320\2403265-69/fedorova_367582/lab6/lab_6" | 246 +++++++++++++++++++ 1 file changed, 246 insertions(+) create mode 100644 "\320\2403265-69/fedorova_367582/lab6/lab_6" diff --git "a/\320\2403265-69/fedorova_367582/lab6/lab_6" "b/\320\2403265-69/fedorova_367582/lab6/lab_6" new file mode 100644 index 0000000..b18b774 --- /dev/null +++ "b/\320\2403265-69/fedorova_367582/lab6/lab_6" @@ -0,0 +1,246 @@ +import math +import numpy as np +import matplotlib.pyplot as plt +from collections import defaultdict +from tabulate import tabulate + + +# Определение доступных дифференциальных уравнений +equations = [ + { + 'name': "y' = y + (1 + x)y^2 ", + 'order': 1, # Порядок уравнения + 'f': lambda x, y: y + (1 + x) * y ** 2, # Правая часть уравнения + 'exact': lambda x, C: -1 / ((x + 1) * (C + math.log(x + 1))), # Точное решение + 'C': lambda x0, y0: -1 / (y0[0] * (x0 + 1)) - math.log(x0 + 1) # Постоянная C + }, + { + 'name': "y'' = (2x y')/(x^2 + 1) ", + 'order': 2, + 'f': lambda x, y: [y[1], (2 * x * y[1]) / (x ** 2 + 1)], + 'exact': lambda x: x ** 3 + 3 * x + 1, + 'C': lambda x0, y0: None # Точное решение фиксировано + }, + { + 'name': "y' = (x + y)^2 ", + 'order': 1, + 'f': lambda x, y: [(x + y[0]) ** 2], + 'exact': lambda x, C: math.tan(x - C) - x, + 'C': lambda x0, y0: x0 - math.atan(y0[0] + x0) if isinstance(y0, list) else x0 - math.atan(y0 + x0) + } +] + + +def euler_method(f, x0, y0, h, x_end): + x = x0 # Начальная точка + y = np.array(y0, dtype=float) # Начальное значение y + xs, ys = [x], [y.copy()] # Списки для хранения результатов + + while x < x_end: + h_step = min(h, x_end - x) # Корректировка шага в конце интервала + dy = f(x, y) # Вычисление производной + y_next = y + h_step * np.array(dy) # Формула Эйлера + x += h_step # Переход к следующей точке + y = y_next # Обновление y + xs.append(x) # Сохранение x + ys.append(y.copy()) # Сохранение y + + return xs, ys, [] # Возвращаем x, y и пустой список ошибок + + +def improved_euler_method(f, x0, y0, h, x_end): + x = x0 + y = np.array(y0, dtype=float) + xs, ys = [x], [y.copy()] + while x < x_end: + h_step = min(h, x_end - x) + # шаг Эйлера + k1 = np.array(f(x, y)) + # используем значение в следующей точке + k2 = np.array(f(x + h_step, y + h_step * k1)) + # усреднение + y_next = y + h_step * (k1 + k2) / 2 + x += h_step + y = y_next + xs.append(x) + ys.append(y.copy()) + return xs, ys, [] + + +def milne_method(f, x0, y0, h, x_end, exact, C=None): + # Начальные точки методом Рунге-Кутты 4-го порядка + x = x0 + y = np.array(y0, dtype=float) + xs, ys = [x], [y.copy()] + steps = 3 # Для начала счета требуется задать решения в трех первых точках, + # которые можно получить одношаговыми методами + for _ in range(steps): + k1 = h * np.array(f(x, y)) + k2 = h * np.array(f(x + h / 2, y + k1 / 2)) + k3 = h * np.array(f(x + h / 2, y + k2 / 2)) + k4 = h * np.array(f(x + h, y + k3)) + y += (k1 + 2 * k2 + 2 * k3 + k4) / 6 # Формула Рунге-Кутты 4-го порядка + x += h + xs.append(x) + ys.append(y.copy()) + + # Метод Милна + while x < x_end - 1e-9: + h_step = min(h, x_end - x) + # Прогноз (используем 4 предыдущие точки) + fp = [f(xs[-4], ys[-4])[0], + f(xs[-3], ys[-3])[0], + f(xs[-2], ys[-2])[0], + f(xs[-1], ys[-1])[0]] + y_pred = ys[-4] + 4 * h_step / 3 * (2 * fp[1] - fp[2] + 2 * fp[3]) + # Коррекция + fc = f(x + h_step, y_pred) # f(x_{n+1}, y_pred) + y_corr = ys[-2] + h_step / 3 * (fp[2] + 4 * fp[3] + fc[0]) + y_next = y_corr + x += h_step + y = y_next + xs.append(x) + ys.append(y.copy()) + + # Вычисление ошибки относительно точного решения + errs = [] + if exact: + for x_val, y_val in zip(xs, ys): + if C is not None: + exact_val = exact(x_val, C) + else: + exact_val = exact(x_val) + if isinstance(exact_val, (list, np.ndarray)): + err = max(abs(y_val[i] - exact_val[i]) for i in range(len(y_val))) + else: + err = abs(y_val[0] - exact_val) + errs.append(err) + return xs, ys, max(errs) if errs else 0.0 + + +def runge_rule(method, f, x0, y0, h, x_end, p): + xs_h, ys_h, _ = method(f, x0, y0, h, x_end) # Решение с шагом h + xs_h2, ys_h2, _ = method(f, x0, y0, h / 2, x_end) # Решение с шагом h/2 + + min_len = min(len(xs_h), len(xs_h2[::2])) # Для сравнения в одних точках + errors = [] + + for i in range(min_len): + y1 = ys_h[i][0] # Значение с шагом h + y2 = ys_h2[2 * i][0] # Значение с шагом h/2 (берем каждую вторую точку) + error = abs(y2 - y1) / (2 ** p - 1) #Формула Рунге + errors.append(error) + return max(errors) # Максимальная оценка погрешности + + +def main(): + print("Выберите уравнение:") + for i, eq in enumerate(equations, 1): + print(f"{i}. {eq['name']}") + try: + idx = int(input("Номер уравнения: ")) - 1 + if not 0 <= idx < len(equations): + raise ValueError + except ValueError: + print("Некорректный выбор.") + return + + eq = equations[idx] + print(f"Выбрано уравнение: {eq['name']}") + # Ввод начальных условий + x0 = float(input("x0: ")) + if eq['order'] == 1: + y0_input = float(input("y0: ")) + y0 = [y0_input] + else: + y0_val = float(input("y(x0): ")) + y_prime0 = float(input("y'(x0): ")) + y0 = [y0_val, y_prime0] + + x_end = float(input("x_n: ")) + h = float(input("Шаг h: ")) + eps = float(input("Точность eps: ")) + + # Проверка корректности + if h <= 0 or x_end <= x0 or eps <= 0: + print("Некорректные параметры.") + return + + # Точное решение + exact_solution = eq.get('exact') + C = eq.get('C')(x0, y0) if 'C' in eq else None + + # Решение методами + euler_x, euler_y, _ = euler_method(eq['f'], x0, y0, h, x_end) + impr_x, impr_y, _ = improved_euler_method(eq['f'], x0, y0, h, x_end) + milne_x, milne_y, milne_err = milne_method(eq['f'], x0, y0, h, x_end, exact_solution, C) + + # Формирование таблицы + print("\nРезультаты:") + headers = ["x", "Euler", "Improved Euler", "Milne", "Exact"] + data = [] + + all_x = sorted(set(euler_x + impr_x + milne_x)) + for x in all_x: + row = [f"{x:.4f}"] + # Euler + euler_val = next((y[0] for xi, y in zip(euler_x, euler_y) if abs(xi - x) < 1e-9), None) + row.append(f"{euler_val:.4f}" if euler_val is not None else "---") + # Improved Euler + impr_val = next((y[0] for xi, y in zip(impr_x, impr_y) if abs(xi - x) < 1e-9), None) + row.append(f"{impr_val:.4f}" if impr_val is not None else "---") + # Milne + milne_val = next((y[0] for xi, y in zip(milne_x, milne_y) if abs(xi - x) < 1e-9), None) + row.append(f"{milne_val:.4f}" if milne_val is not None else "---") + # Exact + if exact_solution: + try: + exact_val = exact_solution(x, C) + except TypeError: + exact_val = exact_solution(x) + if isinstance(exact_val, (list, np.ndarray)): + exact_val = exact_val[0] + row.append(f"{exact_val:.4f}") + else: + row.append("---") + data.append(row) + + # Вывод таблицы + print(tabulate(data, headers=headers, tablefmt="simple")) + + # Оценка точности + print(f"\nМаксимальная погрешность Милна: {milne_err:.6f}") + + # Оценка по Рунге + euler_runge_err = runge_rule(euler_method, eq['f'], x0, y0, h, x_end, p=1) + impr_runge_err = runge_rule(improved_euler_method, eq['f'], x0, y0, h, x_end, p=2) + + print(f"\nПогрешность Эйлера по правилу Рунге: {euler_runge_err:.6f}") + print(f"Погрешность улучшенного Эйлера по правилу Рунге: {impr_runge_err:.6f}") + + # Графики + plt.figure(figsize=(10, 6)) + + plt.plot(euler_x, [y[0] for y in euler_y], label='Euler', linestyle='--', marker='o') + plt.plot(impr_x, [y[0] for y in impr_y], label='Improved Euler', linestyle='-.', marker='s') + plt.plot(milne_x, [y[0] for y in milne_y], label='Milne', linestyle='-', marker='^') + + # Точное решение (если есть) + if exact_solution: + x_vals = np.linspace(x0, x_end, 200) + y_exact = [exact_solution(x, C) if C is not None else exact_solution(x) for x in x_vals] + if isinstance(y_exact[0], (list, np.ndarray)): + y_exact = [y[0] for y in y_exact] + plt.plot(x_vals, y_exact, '--', label='Exact', color='black') + + plt.xlabel('x') + plt.ylabel('y') + plt.title("Сравнение численных методов решения ОДУ") + plt.legend() + plt.grid(True) + plt.tight_layout() + plt.show() + + +if __name__ == "__main__": + main()