diff --git "a/\320\2403212/bondarev_408300/lab1/input.txt" "b/\320\2403212/bondarev_408300/lab1/input.txt" new file mode 100644 index 0000000..285ef82 --- /dev/null +++ "b/\320\2403212/bondarev_408300/lab1/input.txt" @@ -0,0 +1,8 @@ +3 +4 1 1 +1 5 1 +1 1 6 +6 7 8 + + + diff --git "a/\320\2403212/bondarev_408300/lab1/main.py" "b/\320\2403212/bondarev_408300/lab1/main.py" new file mode 100644 index 0000000..7293a83 --- /dev/null +++ "b/\320\2403212/bondarev_408300/lab1/main.py" @@ -0,0 +1,179 @@ +import numpy as np +import random + + +def read_data_from_keyboard(): + print("Введите размерность матрицы n (1-20): ") + n = int(input().strip()) + if n > 20 or n <= 0: + raise ValueError("Недопустимое значение n.") + + A = [] + print(f"Введите матрицу A ({n}x{n} построчно, числа через пробел):") + for _ in range(n): + row = list(map(float, input().split())) + if len(row) != n: + raise ValueError("Неверное количество элементов в строке.") + A.append(row) + + print(f"Введите вектор b ({n} чисел через пробел):") + b = list(map(float, input().split())) + if len(b) != n: + raise ValueError("Неверное количество элементов в векторе b.") + + return n, A, b + + +def read_data_from_file(filename="input.txt"): + with open(filename, 'r') as f: + lines = [line.strip() for line in f if line.strip()] + + idx = 0 + n = int(lines[idx]) + idx += 1 + A = [] + for _ in range(n): + A.append(list(map(float, lines[idx].split()))) + idx += 1 + + b = list(map(float, lines[idx].split())) + return n, A, b + + +def generate_random_system(n): + A = [] + for i in range(n): + row = [random.uniform(-10, 10) for _ in range(n)] + diag_temp = row[i] + row[i] = 0 + sum_other = sum(abs(x) for x in row) + sign = random.choice([-1, 1]) + row[i] = sign * (sum_other + random.uniform(1, 5)) + A.append(row) + + b = [random.uniform(-10, 10) for _ in range(n)] + return np.array(A, dtype=float), np.array(b, dtype=float) + + +def check_and_make_diagonally_dominant(A, b): + n = len(A) + for i in range(n): + max_row = i + max_val = abs(A[i][i]) + for j in range(i, n): + current = abs(A[j][i]) + if current > max_val: + max_val = current + max_row = j + + if max_val == 0: + return False + + A[i], A[max_row] = A[max_row], A[i] + b[i], b[max_row] = b[max_row], b[i] + return True + + +def simple_iterations(A, b, eps=1e-6, max_iter=1000): + n = len(A) + D = np.diag(A) + if any(d == 0 for d in D): + raise ValueError("На диагонали есть нулевые элементы после перестановки.") + + B = -np.array(A) / D[:, None] + np.fill_diagonal(B, 0) + c = np.array(b) / D + + x = np.zeros(n) + for it in range(1, max_iter + 1): + x_new = B @ x + c + if np.linalg.norm(x_new - x, np.inf) < eps: + return x_new, it, (x_new - x) + x = x_new + return x, max_iter, (x - x_new) + + +def print_results(A, b, x, iterations, diff, eps): + print("\nРезультаты:") + print(f"Количество итераций: {iterations}") + print("Вектор решения:") + print(np.array2string(x, precision=9, suppress_small=True)) + + print("\nВектор погрешностей (последняя итерация):") + print(np.array2string(diff, precision=2, suppress_small=True)) + + residual = A @ x - b + print("\nВектор невязки (A*x - b):") + print(np.array2string(residual, precision=2, suppress_small=True)) + + print(f"\nПроверка точности (||Δx|| < {eps}): {np.linalg.norm(diff, np.inf):.2e}") + + +def compare_with_numpy(A, b): + try: + x_lib = np.linalg.solve(A, b) + print("\nРешение с использованием numpy.linalg.solve:") + print(np.array2string(x_lib, precision=9, suppress_small=True)) + return x_lib + except np.linalg.LinAlgError: + print("\nМатрица вырождена, библиотечное решение невозможно.") + return None + + +def main(): + print("Лабораторная работа: Метод простых итераций") + print("=" * 50) + + source = input("Ввод данных:\n1 - Клавиатура\n2 - Файл\n3 - Случайная генерация\nВыбор: ").strip() + + if source == '1': + n, A, b = read_data_from_keyboard() + A = np.array(A, dtype=float) + b = np.array(b, dtype=float) + elif source == '2': + n, A, b = read_data_from_file() + A = np.array(A, dtype=float) + b = np.array(b, dtype=float) + elif source == '3': + print("Введите размерность матрицы n (1-20): ") + n = int(input().strip()) + if n > 20 or n <= 0: + print("Недопустимое значение n!") + return + A, b = generate_random_system(n) + print("\nСгенерированная матрица A:") + print(np.array2string(A, precision=2, suppress_small=True)) + print("\nСгенерированный вектор b:") + print(np.array2string(b, precision=2, suppress_small=True)) + else: + print("Неверный ввод!") + return + + if not check_and_make_diagonally_dominant(A, b): + print("Невозможно достичь диагонального преобладания!") + return + + eps = float(input("Введите точность (например, 1e-6): ").strip() or 1e-6) + + try: + x, iterations, diff = simple_iterations(A, b, eps) + except ValueError as e: + print(f"Ошибка: {str(e)}") + return + + print_results(A, b, x, iterations, diff, eps) + x_lib = compare_with_numpy(A, b) + + if x_lib is not None: + diff = x - x_lib + print("\nРазница с библиотечным решением:") + print(f"Максимальная: {np.abs(diff).max():.2e}") + print(f"Средняя: {np.abs(diff).mean():.2e}") + print("\nОбъяснение: Различия вызваны:") + print("- Ограниченным числом итераций метода") + print("- Накоплением ошибок округления") + print("- Особенностями метода (диагональное преобладание)") + + +if __name__ == "__main__": + main() \ No newline at end of file diff --git "a/\320\2403212/bondarev_408300/lab1/\320\222\321\213\321\207\320\274\320\260\321\202_\320\233\320\2401_\320\221\320\276\320\275\320\264\320\260\321\200\320\265\320\262\320\220\320\234.pdf" "b/\320\2403212/bondarev_408300/lab1/\320\222\321\213\321\207\320\274\320\260\321\202_\320\233\320\2401_\320\221\320\276\320\275\320\264\320\260\321\200\320\265\320\262\320\220\320\234.pdf" new file mode 100644 index 0000000..c6edfaa Binary files /dev/null and "b/\320\2403212/bondarev_408300/lab1/\320\222\321\213\321\207\320\274\320\260\321\202_\320\233\320\2401_\320\221\320\276\320\275\320\264\320\260\321\200\320\265\320\262\320\220\320\234.pdf" differ diff --git "a/\320\2403212/bondarev_408300/lab2/CalMath lab2 BondarevAM.pdf" "b/\320\2403212/bondarev_408300/lab2/CalMath lab2 BondarevAM.pdf" new file mode 100644 index 0000000..fdd88d9 Binary files /dev/null and "b/\320\2403212/bondarev_408300/lab2/CalMath lab2 BondarevAM.pdf" differ diff --git "a/\320\2403212/bondarev_408300/lab2/functions.py" "b/\320\2403212/bondarev_408300/lab2/functions.py" new file mode 100644 index 0000000..1c407fd --- /dev/null +++ "b/\320\2403212/bondarev_408300/lab2/functions.py" @@ -0,0 +1,29 @@ +import math + +def f1(x): + """ 2.74x^3 -1.93x^2 -15.28x -3.72 """ + return 2.74*(x**3) - 1.93*(x**2) - 15.28*x - 3.72 + +def df1(x): + """ Производная для f1 """ + return 3*2.74*(x**2) - 2*1.93*x - 15.28 + +def f2(x): + """ sin(x) - 0.5x """ + return math.sin(x) - 0.5*x + +def df2(x): + return math.cos(x) - 0.5 + +def f3(x): + """ e^x - x^2 + 1 """ + return math.e**x - x**2 + 1 + +def df3(x): + return math.e**x - 2*x + +FUNCTIONS = { + "1": (f1, df1, "f1(x)=2.74x^3 -1.93x^2 -15.28x -3.72"), + "2": (f2, df2, "f2(x)=sin(x) - 0.5x"), + "3": (f3, df3, "f3(x)=e^x - x^2 +1") +} diff --git "a/\320\2403212/bondarev_408300/lab2/graphs.py" "b/\320\2403212/bondarev_408300/lab2/graphs.py" new file mode 100644 index 0000000..a812890 --- /dev/null +++ "b/\320\2403212/bondarev_408300/lab2/graphs.py" @@ -0,0 +1,40 @@ +import numpy as np +import matplotlib +matplotlib.use('TkAgg') # Для некоторых сред IDE (PyCharm и т.п.) +import matplotlib.pyplot as plt + +def f(x): + return 2.74*(x**3) - 1.93*(x**2) - 15.28*x - 3.72 + +x_min, x_max = -4, 4 + +xs = np.linspace(x_min, x_max, 400) +ys = [f(xx) for xx in xs] + +plt.figure(figsize=(8,6)) +plt.title("График f(x) = 2.74x^3 - 1.93x^2 - 15.28x - 3.72; \n Интервал x ∈ [-4, 4]") + +# Основная кривая: +plt.plot(xs, ys, label="f(x)") +# Горизонтальная ось x (ось f(x)=0): +plt.axhline(0, color='k', linewidth=1) + +# Ограничение диапазонов осей: +plt.xlim(x_min, x_max) +# y-границы подберем автоматически или умножим при желании: +# plt.ylim(min(ys)*1.1, max(ys)*1.1) + +# Подписи осей: +plt.xlabel("x in [-4,4]") +plt.ylabel("f(x)") + +# Несколько «контрольных» точек с подписями: +test_points = [-3, -2, -1, 0, 1, 2, 3] +for tp in test_points: + val = f(tp) + plt.plot(tp, val, 'ro') # 'r'=red,'o'=marker + plt.text(tp, val, f"({tp}, {val:.1f})", fontsize=8) + +plt.grid(True) +plt.legend() +plt.show() diff --git "a/\320\2403212/bondarev_408300/lab2/input_equation.txt" "b/\320\2403212/bondarev_408300/lab2/input_equation.txt" new file mode 100644 index 0000000..60eb866 --- /dev/null +++ "b/\320\2403212/bondarev_408300/lab2/input_equation.txt" @@ -0,0 +1,2 @@ +2 3 +0.01 diff --git "a/\320\2403212/bondarev_408300/lab2/input_system.txt" "b/\320\2403212/bondarev_408300/lab2/input_system.txt" new file mode 100644 index 0000000..9a1f866 --- /dev/null +++ "b/\320\2403212/bondarev_408300/lab2/input_system.txt" @@ -0,0 +1,2 @@ +1 2 +0.001 diff --git "a/\320\2403212/bondarev_408300/lab2/io_utils.py" "b/\320\2403212/bondarev_408300/lab2/io_utils.py" new file mode 100644 index 0000000..1a50c0e --- /dev/null +++ "b/\320\2403212/bondarev_408300/lab2/io_utils.py" @@ -0,0 +1,70 @@ + +def read_ab_eps_keyboard(): + print("Введите a,b:") + line= input("> ").strip() + a_str,b_str= line.split() + a,b= float(a_str), float(b_str) + if a>b: + a,b= b,a + print("Введите eps:") + eps= float(input("> ").strip()) + return a,b,eps + +def read_ab_eps_file(filename="input_equation.txt"): + try: + with open(filename,"r",encoding="utf-8") as f_in: + line1= f_in.readline().strip() + a_str,b_str= line1.split() + a,b= float(a_str), float(b_str) + if a>b: + a,b= b,a + eps_str= f_in.readline().strip() + eps= float(eps_str) + return a,b,eps + except: + return None,None,None + +def read_xy_eps_keyboard(): + print("Введите x0,y0:") + line= input("> ").strip() + x0_str,y0_str= line.split() + x0,y0= float(x0_str), float(y0_str) + print("Введите eps:") + eps= float(input("> ").strip()) + return x0,y0,eps + +def read_xy_eps_file(filename="input_system.txt"): + try: + with open(filename,"r",encoding="utf-8") as f_in: + line1= f_in.readline().strip() + x0_str,y0_str= line1.split() + x0,y0= float(x0_str), float(y0_str) + eps_str= f_in.readline().strip() + eps= float(eps_str) + return x0,y0,eps + except: + return None,None,None + +def output_root(root, fx, iters): + print("Куда вывести результат? (1-экран, 2-файл)") + ans= input("> ").strip() + if ans=="2": + with open("output.txt","a",encoding="utf-8") as f_out: + f_out.write(f"Корень x={root:.6f}, f(x)={fx:.3e}, итераций={iters}\n") + print("Записано в output.txt") + else: + print(f"Корень x={root:.6f}, f(x)={fx:.3e}, итераций={iters}") + +def output_system_solution(X,Y, iters, F1, F2): + r1= abs(F1(X,Y)) + r2= abs(F2(X,Y)) + print("Куда вывести результат? (1-экран, 2-файл)") + ans= input("> ").strip() + if ans=="2": + with open("output.txt","a",encoding="utf-8") as f_out: + f_out.write(f"Решение системы: x={X:.6f}, y={Y:.6f}, итераций={iters}\n") + f_out.write(f"Вектор погрешностей: |F1|={r1:.3e}, |F2|={r2:.3e}\n") + print("Записано в output.txt") + else: + print(f"Решение системы: x={X:.6f}, y={Y:.6f}, итераций={iters}") + print(f"Вектор погрешностей: |F1|={r1:.3e}, |F2|={r2:.3e}") diff --git "a/\320\2403212/bondarev_408300/lab2/main.py" "b/\320\2403212/bondarev_408300/lab2/main.py" new file mode 100644 index 0000000..1ac4ef2 --- /dev/null +++ "b/\320\2403212/bondarev_408300/lab2/main.py" @@ -0,0 +1,166 @@ + +from functions import FUNCTIONS +from methods import ( + bisection, secant, newton_method, simple_iteration, + check_interval_for_root, check_single_iter_convergence +) +from systems import ( + SYSTEMS, newton_system, simple_iter_system, + check_system_iter_convergence +) +from plot_utils import plot_function, plot_system +from io_utils import ( + read_ab_eps_keyboard, read_ab_eps_file, + read_xy_eps_keyboard, read_xy_eps_file, + output_root, output_system_solution +) + +def main(): + print("=== Программа: решить одиночное уравнение или систему ===") + print("1) Одиночное уравнение") + print("2) Система (3 варианта)") + choice= input("> ").strip() + + if choice=="1": + print("Выберите уравнение:") + for key,(f,df,desc) in FUNCTIONS.items(): + print(f"{key}) {desc}") + fch= input("> ").strip() + if fch not in FUNCTIONS: + fch="1" + f, df, desc= FUNCTIONS[fch] + print(f"Вы выбрали:", desc) + + print("Откуда брать [a,b,eps]? (1 - клавиатура, 2 - файл)") + ans= input("> ").strip() + if ans=="2": + a,b,eps= read_ab_eps_file("input_equation.txt") + if a is None: + print("Не удалось прочесть файл, переходим к клавиатуре.") + a,b,eps= read_ab_eps_keyboard() + else: + a,b,eps= read_ab_eps_keyboard() + + c= check_interval_for_root(f,a,b) + if c==0: + print("На [a,b] нет корня.") + return + elif c>1: + print("На [a,b] несколько корней.") + return + + print("Выберите метод:") + print("1) Бисекция\n2) Секущие\n3) Ньютона\n4) Итерация") + mm= input("> ").strip() + if mm not in ["1","2","3","4"]: + mm="1" + + try: + if mm=="1": + root,iters= bisection(f,a,b, eps=eps) + elif mm=="2": + root,iters= secant(f,a,b, eps=eps) + elif mm=="3": + x0= (a+b)/2 + root,iters= newton_method(f, df, x0, eps=eps) + else: + ok,mv= check_single_iter_convergence(f, df, a,b, alpha=0.01) + if not ok: + print(f"Внимание: max|phi'(x)|={mv:.3f}>=1 => может не сойтись. Продолжить? (y/n)") + c2= input("> ").strip().lower() + if c2!="y": + return + x0= (a+b)/2 + root,iters= simple_iteration(f, df, a,b, x0, eps=eps, alpha=0.01) + except Exception as e: + print("Ошибка вычисления:", e) + return + + fx= f(root) + output_root(root, fx, iters) + + print("Построить график f(x) на [a,b]? (y/n)") + ga= input("> ").strip().lower() + if ga=="y": + plot_function(f,a,b) + + elif choice=="2": + print("Выберите систему:") + for key,info in SYSTEMS.items(): + print(f"{key}) {info['name']}") + sch= input("> ").strip() + if sch not in SYSTEMS: + sch="1" + sys_data= SYSTEMS[sch] + F1= sys_data["F1"] + F2= sys_data["F2"] + dF1dx= sys_data["dF1dx"] + dF1dy= sys_data["dF1dy"] + dF2dx= sys_data["dF2dx"] + dF2dy= sys_data["dF2dy"] + print(f"Выбрана:", sys_data["name"]) + + print("Построить график F1=0, F2=0 сейчас? (y/n)") + ga= input("> ").strip().lower() + if ga=="y": + print("Введите x_min,x_max (напр. -3 3):") + line1= input("> ").strip() + xm_str,xM_str= line1.split() + x_min,x_max= float(xm_str), float(xM_str) + + print("Введите y_min,y_max (напр. -3 3):") + line2= input("> ").strip() + ym_str,yM_str= line2.split() + y_min,y_max= float(ym_str), float(yM_str) + + + plot_system(F1, F2, x_min, x_max, y_min, y_max, []) + + + print("Откуда брать x0,y0,eps? (1 - клавиатура, 2 - файл)") + ans= input("> ").strip() + if ans=="2": + x0,y0,eps= read_xy_eps_file("input_system.txt") + if x0 is None: + print("Не удалось прочесть, переходим к клавиатуре.") + x0,y0,eps= read_xy_eps_keyboard() + else: + x0,y0,eps= read_xy_eps_keyboard() + + print("Метод? (6 - Ньютона, 7 - итерация)") + mm= input("> ").strip() + if mm not in ["6","7"]: + mm="6" + + try: + if mm=="6": + X,Y,iters= newton_system(F1,F2, dF1dx,dF1dy,dF2dx,dF2dy, + x0,y0, eps=eps) + else: + print("Для проверки итерации укажите [x_min,x_max,y_min,y_max], напр. -3 3 -3 3") + line3= input("> ").strip() + xmi_str,xma_str, ymi_str,yma_str= line3.split() + xmi,xma= float(xmi_str), float(xma_str) + ymi,yma= float(ymi_str), float(yma_str) + + ok, val= check_system_iter_convergence(dF1dx,dF1dy, dF2dx,dF2dy, + xmi,xma, ymi,yma, + alpha=0.01) + if not ok: + print(f"Внимание: max row-sum={val:.3f}>=1 => может не сойтись. Продолжить? (y/n)") + c2= input("> ").strip().lower() + if c2!="y": + return + + X,Y,iters= simple_iter_system(F1,F2, x0,y0, alpha=0.01, eps=eps) + + output_system_solution(X,Y, iters, F1,F2) + + except Exception as e: + print("Ошибка при решении системы:", e) + + else: + print("Неверный выбор. Завершаем.") + +if __name__=="__main__": + main() diff --git "a/\320\2403212/bondarev_408300/lab2/methods.py" "b/\320\2403212/bondarev_408300/lab2/methods.py" new file mode 100644 index 0000000..b872a21 --- /dev/null +++ "b/\320\2403212/bondarev_408300/lab2/methods.py" @@ -0,0 +1,94 @@ +import math +import numpy as np + +def sign_changes_count(f, a, b, steps=2000): + xs = np.linspace(a, b, steps+1) + prev_sign = None + count = 0 + for x in xs: + val = f(x) + s = 1 if val>0 else (-1 if val<0 else 0) + if prev_sign is not None and s*prev_sign<0: + count += 1 + if s != 0: + prev_sign = s + return count + +def check_interval_for_root(f, a, b): + c = sign_changes_count(f, a, b) + return c + +def bisection(f, a, b, eps=1e-3, max_iter=100): + fa = f(a) + fb = f(b) + if fa*fb > 0: + raise ValueError("На [a,b] нет гарантии одного корня (f(a)*f(b)>0).") + + it = 0 + while it < max_iter and abs(b - a) > eps: + it += 1 + x = (a + b)/2 + fx = f(x) + if fa*fx <= 0: + b = x + fb = fx + else: + a = x + fa = fx + return ((a + b)/2, it) + +def secant(f, a, b, eps=1e-3, max_iter=100): + x0, x1 = a, b + for i in range(max_iter): + fx0, fx1 = f(x0), f(x1) + denom = fx1 - fx0 + if abs(denom) < 1e-15: + raise ZeroDivisionError("Секущая выродилась.") + x2 = x1 - fx1*(x1 - x0)/denom + if abs(x2 - x1) < eps: + return (x2, i+1) + x0, x1 = x1, x2 + return (x1, max_iter) + +def newton_method(f, df, x0, eps=1e-3, max_iter=100): + x = x0 + for i in range(max_iter): + fx = f(x) + dfx = df(x) + if abs(dfx) < 1e-15: + raise ZeroDivisionError("Производная ~0 возле x=%.5f" % x) + x_new = x - fx/dfx + if abs(x_new - x) < eps: + return (x_new, i+1) + x = x_new + return (x, max_iter) + + +def phi_simple(x, f, alpha=0.01): + return x - alpha*f(x) + +def simple_iteration(f, df, a, b, x0, eps=1e-3, max_iter=100, alpha=0.01): + x = x0 + for i in range(max_iter): + x_new = phi_simple(x, f, alpha=alpha) + if abs(x_new - x) < eps: + return (x_new, i+1) + x = x_new + return (x, max_iter) + +def check_single_iter_convergence(f, df, a, b, alpha=0.01, steps=100): + """ + Проверяем достаточное условие сходимости для + phi(x)= x - alpha*f(x), т.е. phi'(x)= 1 - alpha*f'(x). + Проверяем max|1 - alpha*f'(x)|<1 на [a,b]. + Возвращаем (ok, max_val). + """ + xs = np.linspace(a, b, steps+1) + max_val = 0.0 + for xx in xs: + dphi = 1 - alpha*df(xx) + val = abs(dphi) + if val > max_val: + max_val = val + is_ok = (max_val < 1) + return (is_ok, max_val) diff --git "a/\320\2403212/bondarev_408300/lab2/output.txt" "b/\320\2403212/bondarev_408300/lab2/output.txt" new file mode 100644 index 0000000..e7af439 --- /dev/null +++ "b/\320\2403212/bondarev_408300/lab2/output.txt" @@ -0,0 +1,6 @@ +Корень x=2.825065, f(x)=-5.120e-01, итераций=7 +Решение системы: x=3.545842, y=-0.769501, итераций=100 +Вектор погрешностей: |F1|=1.417e+00, |F2|=5.810e+00 +Корень x=2.837971, f(x)=2.719e-04, итераций=3 +Решение системы: x=4.545642, y=-0.567501, итераций=79 +Вектор погрешностей: |F1|=1.407e+00, |F2|=4.790e+00 diff --git "a/\320\2403212/bondarev_408300/lab2/plot_utils.py" "b/\320\2403212/bondarev_408300/lab2/plot_utils.py" new file mode 100644 index 0000000..34a3b2e --- /dev/null +++ "b/\320\2403212/bondarev_408300/lab2/plot_utils.py" @@ -0,0 +1,104 @@ +import numpy as np +import matplotlib +matplotlib.use('TkAgg') +import matplotlib.pyplot as plt +import matplotlib.lines as mlines + +def plot_function(f, a, b): + def f(x): + return 2.74*(x**3) - 1.93*(x**2) - 15.28*x - 3.72 + + x_min, x_max = -4, 4 + + xs = np.linspace(x_min, x_max, 400) + ys = [f(xx) for xx in xs] + + plt.figure(figsize=(8,6)) + plt.title("График f(x) = 2.74x^3 - 1.93x^2 - 15.28x - 3.72;\nИнтервал x ∈ [-4, 4]") + + plt.plot(xs, ys, label="f(x)") + plt.axhline(0, color='k', linewidth=1) + + plt.xlim(x_min, x_max) + plt.xlabel("x in [-4,4]") + plt.ylabel("f(x)") + + test_points = [-3, -2, -1, 0, 1, 2, 3] + for tp in test_points: + val = f(tp) + plt.plot(tp, val, 'ro') + plt.text(tp, val, f"({tp}, {val:.1f})", fontsize=8) + + plt.grid(True) + plt.legend() + plt.show() + + +def plot_function_with_root(f, a, b, root): + xs = np.linspace(a, b, 400) + ys = [f(xx) for xx in xs] + + plt.figure(figsize=(8,6)) + plt.title(f"График f(x) на [{a},{b}] с отмеченным корнем") + plt.plot(xs, ys, label="f(x)") + plt.axhline(0, color='k', linewidth=1) + + plt.xlim(a, b) + min_y = min(ys) + max_y = max(ys) + margin = 0.1*(max_y - min_y) if max_y>min_y else 1 + plt.ylim(min_y - margin, max_y + margin) + + plt.xlabel(f"x in [{a},{b}]") + plt.ylabel("f(x)") + + fxr = f(root) + plt.plot(root, fxr, 'ro') + plt.text(root, fxr, f" root=({root:.3f}, {fxr:.3f})", fontsize=9) + + plt.grid(True) + plt.legend() + plt.show() + + +def plot_system(F1, F2, + x_min, x_max, + y_min, y_max, + solutions): + + xs = np.linspace(x_min, x_max, 300) + ys = np.linspace(y_min, y_max, 300) + X, Y = np.meshgrid(xs, ys) + + Z1 = np.zeros_like(X) + Z2 = np.zeros_like(X) + + rows, cols = X.shape + for i in range(rows): + for j in range(cols): + xx = X[i,j] + yy = Y[i,j] + Z1[i,j] = F1(xx, yy) + Z2[i,j] = F2(xx, yy) + + plt.figure(figsize=(8,6)) + plt.title(f"График системы F1=sin(x+1)-y=1.2; \nF2=2x+cos(y)=2; \nx∈[{x_min},{x_max}], y∈[{y_min},{y_max}]") + + plt.xlim(x_min, x_max) + plt.ylim(y_min, y_max) + + c1 = plt.contour(X, Y, Z1, levels=[0], colors='blue') + c2 = plt.contour(X, Y, Z2, levels=[0], colors='red') + + for (x_s, y_s) in solutions: + plt.plot(x_s, y_s, 'mo') # magenta circle + plt.text(x_s, y_s, f"({x_s:.3f}, {y_s:.3f})", fontsize=9) + + lineF1 = mlines.Line2D([], [], color='blue', label='F1=0') + lineF2 = mlines.Line2D([], [], color='red', label='F2=0') + plt.legend(handles=[lineF1, lineF2], loc='upper right') + + plt.xlabel(f"x in [{x_min},{x_max}]") + plt.ylabel(f"y in [{y_min},{y_max}]") + plt.grid(True) + plt.show() diff --git "a/\320\2403212/bondarev_408300/lab2/systems.py" "b/\320\2403212/bondarev_408300/lab2/systems.py" new file mode 100644 index 0000000..c821961 --- /dev/null +++ "b/\320\2403212/bondarev_408300/lab2/systems.py" @@ -0,0 +1,131 @@ +import math +import numpy as np + +# ----- СИСТЕМА №1 ----- +def F1_1(x,y): + return math.sin(x+1) - y -1.2 +def F2_1(x,y): + return 2*x + math.cos(y) -2 + +def dF1_1dx(x,y): + return math.cos(x+1) +def dF1_1dy(x,y): + return -1 +def dF2_1dx(x,y): + return 2 +def dF2_1dy(x,y): + return -math.sin(y) + +# ----- СИСТЕМА №2 ----- +def F1_2(x,y): + return x*x + y*y -4 +def F2_2(x,y): + return x - y -1 + +def dF1_2dx(x,y): + return 2*x +def dF1_2dy(x,y): + return 2*y +def dF2_2dx(x,y): + return 1 +def dF2_2dy(x,y): + return -1 + +# ----- СИСТЕМА №3 ----- +def F1_3(x,y): + return math.e**x - y +def F2_3(x,y): + return x + y -2 + +def dF1_3dx(x,y): + return math.e**x +def dF1_3dy(x,y): + return -1 +def dF2_3dx(x,y): + return 1 +def dF2_3dy(x,y): + return 1 + +SYSTEMS = { + "1": { + "name":"System1: { sin(x+1)-y=1.2, 2x+cos(y)=2 }", + "F1": F1_1, "F2": F2_1, + "dF1dx": dF1_1dx, "dF1dy": dF1_1dy, + "dF2dx": dF2_1dx, "dF2dy": dF2_1dy + }, + "2": { + "name":"System2: { x^2+y^2=4, x-y=1 }", + "F1": F1_2, "F2": F2_2, + "dF1dx": dF1_2dx, "dF1dy": dF1_2dy, + "dF2dx": dF2_2dx, "dF2dy": dF2_2dy + }, + "3": { + "name":"System3: { e^x - y=0, x+y=2 }", + "F1": F1_3, "F2": F2_3, + "dF1dx": dF1_3dx, "dF1dy": dF1_3dy, + "dF2dx": dF2_3dx, "dF2dy": dF2_3dy + } +} + +def newton_system(F1,F2, dF1dx,dF1dy,dF2dx,dF2dy, + x0,y0, eps=1e-3, max_iter=100): + x,y= x0,y0 + for i in range(max_iter): + J11= dF1dx(x,y) + J12= dF1dy(x,y) + J21= dF2dx(x,y) + J22= dF2dy(x,y) + det= J11*J22 - J12*J21 + if abs(det)<1e-15: + raise ZeroDivisionError("Якобиан вырожден.") + Fx= -F1(x,y) + Fy= -F2(x,y) + dx= (Fx*J22 - Fy*J12)/det + dy= (J11*Fy - J21*Fx)/det + x_new= x+dx + y_new= y+dy + if math.hypot(dx,dy)< eps: + return (x_new,y_new, i+1) + x,y= x_new,y_new + return (x,y, max_iter) + +def simple_iter_system(F1,F2, x0,y0, alpha=0.01, eps=1e-3, max_iter=100): + x,y= x0,y0 + for i in range(max_iter): + x_new= x - alpha*F1(x,y) + y_new= y - alpha*F2(x,y) + if math.hypot(x_new - x, y_new - y)< eps: + return (x_new,y_new, i+1) + x,y= x_new,y_new + return (x,y, max_iter) + +def check_system_iter_convergence(dF1dx,dF1dy,dF2dx,dF2dy, + x_min,x_max,y_min,y_max, + alpha=0.01, steps=10): + """ + Проверяем max-норму Якобиана phi'(x,y) для + phi_1(x,y)= x - alpha*F1(x,y), + phi_2(x,y)= y - alpha*F2(x,y). + + J_phi= [ 1 - alpha*dF1dx -alpha*dF1dy ] + [ -alpha*dF2dx 1-alpha*dF2dy] + + Берём максимум построчной нормы (строка 1= |j11|+|j12|, строка2=|j21|+|j22|). + Возвращаем (ok, max_val). + """ + xs= np.linspace(x_min, x_max, steps+1) + ys= np.linspace(y_min, y_max, steps+1) + max_val= 0.0 + for xx in xs: + for yy in ys: + j11= 1 - alpha*dF1dx(xx,yy) + j12= -alpha*dF1dy(xx,yy) + j21= -alpha*dF2dx(xx,yy) + j22= 1 - alpha*dF2dy(xx,yy) + row1= abs(j11)+ abs(j12) + row2= abs(j21)+ abs(j22) + loc= max(row1,row2) + if loc> max_val: + max_val= loc + is_ok= (max_val<1) + return (is_ok, max_val) diff --git "a/\320\2403212/bondarev_408300/lab3/Calmathlab3 BondarevAM.pdf" "b/\320\2403212/bondarev_408300/lab3/Calmathlab3 BondarevAM.pdf" new file mode 100644 index 0000000..eba8d46 Binary files /dev/null and "b/\320\2403212/bondarev_408300/lab3/Calmathlab3 BondarevAM.pdf" differ diff --git "a/\320\2403212/bondarev_408300/lab3/__init__.py" "b/\320\2403212/bondarev_408300/lab3/__init__.py" new file mode 100644 index 0000000..994c79c --- /dev/null +++ "b/\320\2403212/bondarev_408300/lab3/__init__.py" @@ -0,0 +1,5 @@ +from .rectangles import left_rect_method, right_rect_method, mid_rect_method +from .trapezoid import trapezoid_method +from .simpson import simpson_method +from .runge import integrate_with_precision +from .improper import improper_integral diff --git "a/\320\2403212/bondarev_408300/lab3/functions.py" "b/\320\2403212/bondarev_408300/lab3/functions.py" new file mode 100644 index 0000000..6a800bf --- /dev/null +++ "b/\320\2403212/bondarev_408300/lab3/functions.py" @@ -0,0 +1,19 @@ +import math + +def func1(x): + return -x**3 - x**2 - 2*x + 1 + +def func2(x): + return x**2 + +def func3(x): + return math.sin(x) + +def func4(x): # ln(x)/sqrt(x) (разрыв в a = 0) + return math.log(x) / math.sqrt(x) + +def func5(x): # 1 / sqrt(1-x) (разрыв в b = 1) + return 1.0 / math.sqrt(1.0 - x) + +def func6(x): # 1/sqrt(|x-0.5|) разрыв в точке 0.5 + return 1.0 / math.sqrt(abs(x - 0.5)) diff --git "a/\320\2403212/bondarev_408300/lab3/improper.py" "b/\320\2403212/bondarev_408300/lab3/improper.py" new file mode 100644 index 0000000..e542d6b --- /dev/null +++ "b/\320\2403212/bondarev_408300/lab3/improper.py" @@ -0,0 +1,67 @@ +import math +from .runge import integrate_with_precision + + +def _is_finite(value): + return math.isfinite(value) + +def _safe_eval(f, x): + try: + return f(x) + except Exception: + return float("inf") + +def detect_singularity(f, a, b, n_test=10): + # край a + if not _is_finite(_safe_eval(f, a + 1e-12)): + return "left", a + # край b + if not _is_finite(_safe_eval(f, b - 1e-12)): + return "right", b + # поиск внутреннего + h = (b - a) / (n_test + 1) + for i in range(1, n_test + 1): + x = a + i*h + if not _is_finite(_safe_eval(f, x)): + return "inner", x + return None, None + +def improper_integral(f, a, b, eps, method, p, n_start=4): + kind, c = detect_singularity(f, a, b) + if kind is None: + return integrate_with_precision(method, f, a, b, eps, p, n_start) + + max_iter = 20 + delta = 1e-3 + tail_val = 0.0 + for _ in range(max_iter): + if kind == "left": + tail_val = method(f, a, a+delta, n_start) + elif kind == "right": + tail_val = method(f, b-delta, b, n_start) + else: # "inner" + tail_val = ( + method(f, c-delta, c, n_start) + + method(f, c, c+delta, n_start) + ) + if math.isfinite(tail_val) and abs(tail_val) < eps/2: + break + delta /= 2.0 + else: + return None, None + + if kind == "left": + I_main, n_main = integrate_with_precision( + method, f, a+delta, b, eps/2, p, n_start) + elif kind == "right": + I_main, n_main = integrate_with_precision( + method, f, a, b-delta, eps/2, p, n_start) + else: + I1, n1 = integrate_with_precision( + method, f, a, c-delta, eps/4, p, n_start) + I2, n2 = integrate_with_precision( + method, f, c+delta, b, eps/4, p, n_start) + I_main = I1 + I2 + n_main = n1 + n2 + + return I_main + tail_val, n_main diff --git "a/\320\2403212/bondarev_408300/lab3/main.py" "b/\320\2403212/bondarev_408300/lab3/main.py" new file mode 100644 index 0000000..f9a0ce8 --- /dev/null +++ "b/\320\2403212/bondarev_408300/lab3/main.py" @@ -0,0 +1,110 @@ +from functions import ( + func1, func2, func3, # «обычные» + func4, func5, func6 # с разрывами +) +from methods import ( + left_rect_method, right_rect_method, mid_rect_method, + trapezoid_method, simpson_method, + integrate_with_precision, improper_integral +) + +def choose_function(): + print("Функции:") + print("1) -x^3 - x^2 - 2x + 1") + print("2) x^2") + print("3) sin(x)") + print("4) ln(x)/sqrt(x) (разрыв в a)") + print("5) 1/sqrt(1-x) (разрыв в b)") + print("6) 1/sqrt(|x-0.5|) (разрыв внутри)") + while True: + try: + k = int(input("Выберите номер функции: ").strip()) + if k == 1: + return func1 + if k == 2: + return func2 + if k == 3: + return func3 + if k == 4: + return func4 + if k == 5: + return func5 + if k == 6: + return func6 + raise ValueError + except ValueError: + print("Неверный ввод.") + +def choose_method(): + print("\nМетоды:") + print("1) Левые прямоугольники") + print("2) Правые прямоугольники") + print("3) Средние прямоугольники") + print("4) Трапеции") + print("5) Симпсон") + while True: + try: + k = int(input("Номер метода: ").strip()) + if k == 1: + return left_rect_method, 1, "Левые прямоугольники" + elif k == 2: + return right_rect_method, 1, "Правые прямоугольники" + elif k == 3: + return mid_rect_method, 2, "Средние прямоугольники" + elif k == 4: + return trapezoid_method, 2, "Метод трапеций" + elif k == 5: + return simpson_method, 4, "Метод Симпсона" + raise ValueError + except ValueError: + print("Неверный ввод.") + +def main(): + print("=== Численное интегрирование (включая несобственные) ===\n") + + f = choose_function() + + while True: + try: + a = float(input("\nНижний предел a: ").strip()) + b = float(input("Верхний предел b: ").strip()) + if a == b: + raise ValueError + break + except ValueError: + print("Введите два различных числа.") + + while True: + try: + eps = float(input("Требуемая точность eps (например 1e-6): ").strip()) + if eps <= 0: + raise ValueError + break + except ValueError: + print("eps должно быть положительным числом.") + + method, p_order, method_name = choose_method() + + impro = input("\nНесобственный интеграл? (y/N) ").strip().lower() == 'y' + + if impro: + res, n_tot = improper_integral(f, a, b, eps, method, p_order, n_start=4) + if res is None: + print("\n>>> Интеграл не существует (расходится).") + return + algo_used = "несобственный (с усечением хвоста)" + else: + res, n_tot = integrate_with_precision(method, f, a, b, eps, p_order, n_start=4) + algo_used = "обычный" + + print("\n-------- Результаты --------") + print(f"Метод: {method_name}") + print(f"Алгоритм: {algo_used}") + print(f"Интервал: [{a}, {b}]") + print(f"eps: {eps}") + print(f"Значение I: {res}") + print(f"Разбиений n: {n_tot}") + print("----------------------------\n") + +if __name__ == "__main__": + main() diff --git "a/\320\2403212/bondarev_408300/lab3/rectangles.py" "b/\320\2403212/bondarev_408300/lab3/rectangles.py" new file mode 100644 index 0000000..6dfa35e --- /dev/null +++ "b/\320\2403212/bondarev_408300/lab3/rectangles.py" @@ -0,0 +1,26 @@ +def left_rect_method(f, a, b, n): + h = (b - a) / n + s = 0.0 + x = a + for _ in range(n): + s += f(x) + x += h + return s * h + +def right_rect_method(f, a, b, n): + h = (b - a) / n + s = 0.0 + x = a + h + for _ in range(n): + s += f(x) + x += h + return s * h + +def mid_rect_method(f, a, b, n): + h = (b - a) / n + s = 0.0 + x = a + h/2 + for _ in range(n): + s += f(x) + x += h + return s * h diff --git "a/\320\2403212/bondarev_408300/lab3/runge.py" "b/\320\2403212/bondarev_408300/lab3/runge.py" new file mode 100644 index 0000000..ab670e6 --- /dev/null +++ "b/\320\2403212/bondarev_408300/lab3/runge.py" @@ -0,0 +1,17 @@ +def integrate_with_precision(method, f, a, b, eps, p, n_start=4): + """ + Δ = |I_{2n} - I_n| / (2^p - 1), где p — порядок метода (1, 2 или 4). + Возвращает кортеж (v, n), + где: + - v — найденное приближение интеграла, + - n — число разбиений, при котором достигнута точность. + """ + n = n_start + while True: + i_n = method(f, a, b, n) + i_2n = method(f, a, b, 2 * n) + delta = abs(i_2n - i_n) / (2 ** p - 1) + if delta < eps: + return i_2n, 2 * n + else: + n *= 2 diff --git "a/\320\2403212/bondarev_408300/lab3/simpson.py" "b/\320\2403212/bondarev_408300/lab3/simpson.py" new file mode 100644 index 0000000..d07f73a --- /dev/null +++ "b/\320\2403212/bondarev_408300/lab3/simpson.py" @@ -0,0 +1,14 @@ +def simpson_method(f, a, b, n): + if n % 2 != 0: + n += 1 # делаем n чётным + + h = (b - a) / n + s1 = 0.0 #sum нечётные индексы + s2 = 0.0 #четные + for i in range(1, n): + x_i = a + i*h + if i % 2 == 0: + s2 += f(x_i) + else: + s1 += f(x_i) + return (h/3) * (f(a) + 4*s1 + 2*s2 + f(b)) diff --git "a/\320\2403212/bondarev_408300/lab3/trapezoid.py" "b/\320\2403212/bondarev_408300/lab3/trapezoid.py" new file mode 100644 index 0000000..6916532 --- /dev/null +++ "b/\320\2403212/bondarev_408300/lab3/trapezoid.py" @@ -0,0 +1,7 @@ +def trapezoid_method(f, a, b, n): + h = (b - a) / n + s = 0.0 + for i in range(1, n): + x_i = a + i*h + s += f(x_i) + return (h / 2) * (f(a) + 2*s + f(b)) diff --git "a/\320\2403212/bondarev_408300/lab4/CalMathlab4 BondarevAM.pdf" "b/\320\2403212/bondarev_408300/lab4/CalMathlab4 BondarevAM.pdf" new file mode 100644 index 0000000..fd41fc5 Binary files /dev/null and "b/\320\2403212/bondarev_408300/lab4/CalMathlab4 BondarevAM.pdf" differ diff --git "a/\320\2403212/bondarev_408300/lab4/approximations.py" "b/\320\2403212/bondarev_408300/lab4/approximations.py" new file mode 100644 index 0000000..f365eb9 --- /dev/null +++ "b/\320\2403212/bondarev_408300/lab4/approximations.py" @@ -0,0 +1,176 @@ +from dataclasses import dataclass +import math +import numpy as np + +__all__ = [ + "ModelResult", "fit_all_models", "choose_best", + "fit_linear", "fit_poly2", "fit_poly3", + "fit_exponential", "fit_logarithmic", "fit_power", +] + +#вспомогательные статистические функции + +def _rms(y, y_hat) -> float: + """Среднеквадратическая ошибка (σ).""" + return float(np.sqrt(np.mean((y - y_hat) ** 2))) + +def _r_squared(y, y_hat) -> float: + """Коэффициент детерминации R².""" + ss_tot = np.sum((y - y.mean()) ** 2) + ss_res = np.sum((y - y_hat) ** 2) + return float(1 - ss_res / ss_tot) if ss_tot else float("nan") + +def _pearson_r(x, y) -> float: + """Коэффициент корреляции Пирсона (для линейной модели).""" + if len(x) < 2: + return float("nan") + sx, sy = x.std(ddof=0), y.std(ddof=0) + cov = np.cov(x, y, ddof=0)[0, 1] + return float(cov / (sx * sy)) if sx and sy else float("nan") + + +#структура результата + +@dataclass +class ModelResult: + name: str #имя + formula: str # шаблон формулы для отчёта + coef: np.ndarray # коэффициенты (a, b, c, …) + y_hat: np.ndarray # значения модели в узлах + sigma: float # RMS‑ошибка + r2: float # коэффициент детерминации + extra: dict # доп. показатели (р, …) + + def coef_str(self, digits: int = 5) -> str: + return ", ".join(f"{c:.{digits}g}" for c in self.coef) + + +#МНК‑аппроксимации + +def fit_linear(x, y) -> ModelResult: + X = np.vstack([np.ones_like(x), x]).T + a, b = np.linalg.lstsq(X, y, rcond=None)[0] + y_hat = a + b * x + return ModelResult( + name="Линейная", + formula="y = {a:.5g} + {b:.5g}·x", + coef=np.array([a, b]), + y_hat=y_hat, + sigma=_rms(y, y_hat), + r2=_r_squared(y, y_hat), + extra={"Pearson r": _pearson_r(x, y)}, + ) + + +def fit_poly2(x, y) -> ModelResult: + X = np.vstack([np.ones_like(x), x, x**2]).T + a, b, c = np.linalg.lstsq(X, y, rcond=None)[0] + y_hat = a + b * x + c * x**2 + return ModelResult( + name="Полином 2‑й ст.", + formula="y = {a:.5g} + {b:.5g}·x + {c:.5g}·x²", + coef=np.array([a, b, c]), + y_hat=y_hat, + sigma=_rms(y, y_hat), + r2=_r_squared(y, y_hat), + extra={}, + ) + + +def fit_poly3(x, y) -> ModelResult: + X = np.vstack([np.ones_like(x), x, x**2, x**3]).T + a, b, c, d = np.linalg.lstsq(X, y, rcond=None)[0] + y_hat = a + b * x + c * x**2 + d * x**3 + return ModelResult( + name="Полином 3‑й ст.", + formula="y = {a:.5g} + {b:.5g}·x + {c:.5g}·x² + {d:.5g}·x³", + coef=np.array([a, b, c, d]), + y_hat=y_hat, + sigma=_rms(y, y_hat), + r2=_r_squared(y, y_hat), + extra={}, + ) + + +def fit_exponential(x, y) -> ModelResult: + if np.any(y <= 0): + raise ValueError("Экспоненциальная модель требует y > 0") + Y = np.log(y) + X = np.vstack([np.ones_like(x), x]).T + ln_a, b = np.linalg.lstsq(X, Y, rcond=None)[0] + a = math.exp(ln_a) + y_hat = a * np.exp(b * x) + return ModelResult( + name="Экспоненциальная", + formula="y = {a:.5g}·e^{b:.5g}·x", + coef=np.array([a, b]), + y_hat=y_hat, + sigma=_rms(y, y_hat), + r2=_r_squared(y, y_hat), + extra={}, + ) + + +def fit_logarithmic(x, y) -> ModelResult: + if np.any(x <= 0): + raise ValueError("Логарифмическая модель требует x > 0") + X = np.vstack([np.ones_like(x), np.log(x)]).T + a, b = np.linalg.lstsq(X, y, rcond=None)[0] + y_hat = a + b * np.log(x) + return ModelResult( + name="Логарифмическая", + formula="y = {a:.5g} + {b:.5g}·ln x", + coef=np.array([a, b]), + y_hat=y_hat, + sigma=_rms(y, y_hat), + r2=_r_squared(y, y_hat), + extra={}, + ) + + +def fit_power(x, y) -> ModelResult: + if np.any(x <= 0) or np.any(y <= 0): + raise ValueError("Степенная модель требует x, y > 0") + Y = np.log(y) + X = np.vstack([np.ones_like(x), np.log(x)]).T + ln_a, b = np.linalg.lstsq(X, Y, rcond=None)[0] + a = math.exp(ln_a) + y_hat = a * x ** b + return ModelResult( + name="Степенная", + formula="y = {a:.5g}·x^{b:.5g}", + coef=np.array([a, b]), + y_hat=y_hat, + sigma=_rms(y, y_hat), + r2=_r_squared(y, y_hat), + extra={}, + ) + + +#сборщик моделей + +_MODEL_FUNCS = [ + fit_linear, + fit_poly2, + fit_poly3, + fit_exponential, + fit_logarithmic, + fit_power, +] + + +def fit_all_models(x, y): +#Возвращает список ModelResult для всех шести моделей. +#Неприменимая (из-за ограничений) модель пропускается. + models = [] + for fit_func in _MODEL_FUNCS: + try: + models.append(fit_func(x, y)) + except ValueError as e: + print(f"→ {fit_func.__name__[4:]} пропущена: {e}") + return models + + +def choose_best(models): + """Модель с минимальной σ.""" + return min(models, key=lambda m: m.sigma) diff --git "a/\320\2403212/bondarev_408300/lab4/data.csv" "b/\320\2403212/bondarev_408300/lab4/data.csv" new file mode 100644 index 0000000..ed2170c --- /dev/null +++ "b/\320\2403212/bondarev_408300/lab4/data.csv" @@ -0,0 +1,11 @@ +0.0,0.0 +0.2,2.396 +0.4,4.680 +0.6,6.374 +0.8,6.810 +1.0,6.000 +1.2,4.685 +1.4,3.470 +1.6,2.542 +1.8,1.879 +2.0,1.412 diff --git "a/\320\2403212/bondarev_408300/lab4/io_utils.py" "b/\320\2403212/bondarev_408300/lab4/io_utils.py" new file mode 100644 index 0000000..6b875cf --- /dev/null +++ "b/\320\2403212/bondarev_408300/lab4/io_utils.py" @@ -0,0 +1,73 @@ +import pathlib +from typing import Tuple, List, Callable +#import matplotlib +#matplotlib.use("TkAgg") +import numpy as np +import pandas as pd + +from approximations import ModelResult + + +#ввод точек +def read_csv(path: pathlib.Path) -> Tuple[np.ndarray, np.ndarray]: + df = pd.read_csv(path, header=None, names=["x", "y"]) + if df.shape[1] != 2: + raise ValueError("CSV должен содержать ровно 2 столбца (x,y)") + return df["x"].to_numpy(float), df["y"].to_numpy(float) + + +def read_console() -> Tuple[np.ndarray, np.ndarray]: + print("Введите пары x y (пустая строка — конец ввода):") + pts = [] + while True: + line = input() + if not line.strip(): + break + try: + pts.append(tuple(map(float, line.split()))) + except Exception: + print("Неверный формат строки. Повторите.") + if not pts: + raise ValueError("Точки не введены") + x, y = zip(*pts) + return np.array(x, float), np.array(y, float) + + +def validate_points(x: np.ndarray, y: np.ndarray, min_n=8, max_n=12): + if not (min_n <= len(x) <= max_n): + raise ValueError(f"Нужно от {min_n} до {max_n} точек, сейчас {len(x)}") + + +#результаты +def save_report(path: pathlib.Path, + x: np.ndarray, + y: np.ndarray, + models: List[ModelResult], + best_name: str, + r2_interpreter: Callable[[float], str]): + with open(path, "w", encoding="utf-8") as f: + print("Аппроксимация методом наименьших квадратов", file=f) + print("=" * 60, file=f) + + # краткая сводка + for m in models: + msg = r2_interpreter(m.r2) + line = (f"{m.name:<16} σ={m.sigma:.6g} R²={m.r2:.6g} " + f"({msg}) coef=[{m.coef_str()}]") + if "Pearson r" in m.extra: + line += f" r={m.extra['Pearson r']:.6g}" + print(line, file=f) + + print("-" * 60, file=f) + print(f"Лучшее приближение → {best_name}", file=f) + + # развёрнутые таблицы для каждой модели + for m in models: + print("\n" + m.name, file=f) + print(" i x_i y_i φ(x_i) ε_i", file=f) + eps = y - m.y_hat + for i, (xi, yi, fi, ei) in enumerate(zip(x, y, m.y_hat, eps), 1): + print(f"{i:2d} {xi:9.5g} {yi:9.5g} {fi:11.5g} {ei:11.5g}", + file=f) + + print(f"Отчёт сохранён в {path.resolve()}") diff --git "a/\320\2403212/bondarev_408300/lab4/main.py" "b/\320\2403212/bondarev_408300/lab4/main.py" new file mode 100644 index 0000000..0848174 --- /dev/null +++ "b/\320\2403212/bondarev_408300/lab4/main.py" @@ -0,0 +1,110 @@ +import argparse +import pathlib +import sys + +import matplotlib.pyplot as plt +import numpy as np + +import approximations as appr +import io_utils as io + +#python main.py -i data.csv +#вспомогательное пояснение R² +def interpret_r2(r2: float) -> str: + if np.isnan(r2): + return "R² не определён" + if r2 >= 0.9: + return "отличное согласие" + if r2 >= 0.75: + return "хорошее согласие" + if r2 >= 0.5: + return "удовлетворительное согласие" + return "слабое согласие" + + +#построение графика +def _plot_models(x, y, models): + xmin, xmax = x.min(), x.max() + dx = 0.05 * (xmax - xmin) #запас + xs = np.linspace(xmin - dx, xmax + dx, 400) + + plt.figure(figsize=(8, 5)) + plt.scatter(x, y, c="black", label="исходные точки") + + for m in models: + if m.name == "Экспоненциальная": + ys = m.coef[0] * np.exp(m.coef[1] * xs) + elif m.name == "Логарифмическая": + ys = m.coef[0] + m.coef[1] * np.log(xs) + ys[np.isinf(ys) | np.isnan(ys)] = np.nan + elif m.name == "Степенная": + ys = m.coef[0] * xs ** m.coef[1] + elif m.name.startswith("Полином 3"): + a, b, c, d = m.coef + ys = a + b*xs + c*xs**2 + d*xs**3 + elif m.name.startswith("Полином 2"): + a, b, c = m.coef + ys = a + b*xs + c*xs**2 + else: # линейная + a, b = m.coef + ys = a + b*xs + plt.plot(xs, ys, label=f"{m.name} (σ={m.sigma:.3g})") + + plt.title("Аппроксимация (МНК)") + plt.xlabel("x"), plt.ylabel("y") + plt.grid(True), plt.legend() + plt.tight_layout() + plt.show() + + +def main(): + parser = argparse.ArgumentParser( + description="ЛР-4: аппроксимация МНК (6 моделей)") + parser.add_argument("-i", "--input", help="CSV-файл с точками x,y") + parser.add_argument("-o", "--output", default="results.txt", + help="файл текстового отчёта (по умолчанию results.txt)") + args = parser.parse_args() + + #ввод точек + try: + if args.input: + x, y = io.read_csv(pathlib.Path(args.input)) + else: + x, y = io.read_console() + io.validate_points(x, y) + except Exception as e: + sys.exit(f"Ошибка ввода точек: {e}") + + print(f"Принято {len(x)} точек.") + + #аппроксимация + models = appr.fit_all_models(x, y) + if not models: + sys.exit("Ни одна модель не применима к этим данным.") + + best = appr.choose_best(models) + + #вывод (коэффициенты + интерпретация R²) + print("\nРезультаты:") + for m in models: + r2_msg = interpret_r2(m.r2) + line = (f"{m.name:<16} σ={m.sigma:.6g} R²={m.r2:.6g} " + f"({r2_msg}) coef=[{m.coef_str()}]") + if "Pearson r" in m.extra: + line += f" r={m.extra['Pearson r']:.6g}" + print(line) + print(f"\nЛучшее приближение → {best.name}") + + # --- полный отчёт в файл (с массивами φ(x_i) и ε_i) + try: + io.save_report(pathlib.Path(args.output), x, y, models, best.name, + r2_interpreter=interpret_r2) + except Exception as e: + print(f"Не удалось сохранить отчёт: {e}") + + #график + _plot_models(x, y, models) + + +if __name__ == "__main__": + main() diff --git "a/\320\2403212/bondarev_408300/lab4/results.txt" "b/\320\2403212/bondarev_408300/lab4/results.txt" new file mode 100644 index 0000000..9ba1467 --- /dev/null +++ "b/\320\2403212/bondarev_408300/lab4/results.txt" @@ -0,0 +1,49 @@ +Аппроксимация методом наименьших квадратов +============================================================ +Линейная σ=2.10117 R²=0.0161184 (слабое согласие) coef=[4.0841, -0.42523] r=-0.126958 +Полином 2‑й ст. σ=0.932766 R²=0.806105 (хорошее согласие) coef=[0.88639, 10.234, -5.3296] +Полином 3‑й ст. σ=0.335814 R²=0.974869 (отличное согласие) coef=[-0.43557, 20.736, -19.1, 4.5901] +------------------------------------------------------------ +Лучшее приближение → Полином 3‑й ст. + +Линейная + i x_i y_i φ(x_i) ε_i + 1 0 0 4.0841 -4.0841 + 2 0.2 2.396 3.9991 -1.6031 + 3 0.4 4.68 3.914 0.76595 + 4 0.6 6.374 3.829 2.545 + 5 0.8 6.81 3.744 3.066 + 6 1 6 3.6589 2.3411 + 7 1.2 4.685 3.5739 1.1111 + 8 1.4 3.47 3.4888 -0.018818 + 9 1.6 2.542 3.4038 -0.86177 +10 1.8 1.879 3.3187 -1.4397 +11 2 1.412 3.2337 -1.8217 + +Полином 2‑й ст. + i x_i y_i φ(x_i) ε_i + 1 0 0 0.88639 -0.88639 + 2 0.2 2.396 2.72 -0.32399 + 3 0.4 4.68 4.1272 0.55277 + 4 0.6 6.374 5.1081 1.2659 + 5 0.8 6.81 5.6626 1.1474 + 6 1 6 5.7907 0.20926 + 7 1.2 4.685 5.4925 -0.80751 + 8 1.4 3.47 4.7679 -1.2979 + 9 1.6 2.542 3.617 -1.075 +10 1.8 1.879 2.0396 -0.16063 +11 2 1.412 0.035937 1.3761 + +Полином 3‑й ст. + i x_i y_i φ(x_i) ε_i + 1 0 0 -0.43557 0.43557 + 2 0.2 2.396 2.9844 -0.58838 + 3 0.4 4.68 5.0967 -0.41666 + 4 0.6 6.374 6.1216 0.2524 + 5 0.8 6.81 6.2795 0.53048 + 6 1 6 5.7907 0.20926 + 7 1.2 4.685 4.8756 -0.1906 + 8 1.4 3.47 3.7544 -0.28441 + 9 1.6 2.542 2.6475 -0.10552 +10 1.8 1.879 1.7752 0.10376 +11 2 1.412 1.3579 0.054105 diff --git "a/\320\2403212/bondarev_408300/lab5/CalMathlab5 BondarevAM.pdf" "b/\320\2403212/bondarev_408300/lab5/CalMathlab5 BondarevAM.pdf" new file mode 100644 index 0000000..d2ac33a Binary files /dev/null and "b/\320\2403212/bondarev_408300/lab5/CalMathlab5 BondarevAM.pdf" differ diff --git "a/\320\2403212/bondarev_408300/lab5/__init__.py" "b/\320\2403212/bondarev_408300/lab5/__init__.py" new file mode 100644 index 0000000..e69de29 diff --git "a/\320\2403212/bondarev_408300/lab5/io_utils.py" "b/\320\2403212/bondarev_408300/lab5/io_utils.py" new file mode 100644 index 0000000..b16a3ff --- /dev/null +++ "b/\320\2403212/bondarev_408300/lab5/io_utils.py" @@ -0,0 +1,60 @@ +import csv +import math +from pathlib import Path +from typing import List, Tuple, Callable + +_BUILTIN_FUNCTIONS = { + "sin": math.sin, + "cos": math.cos, + "exp": math.exp, + "log": math.log, +} + +def _read_keyboard() -> Tuple[List[float], List[float]]: + n = int(input("Сколько узлов? ")) + x, y = [], [] + for i in range(n): + x.append(float(input(f"x[{i}] = "))) + y.append(float(input(f"y[{i}] = "))) + return x, y + +def _read_from_file() -> Tuple[List[float], List[float]]: + name = input("Имя CSV-файла: ").strip() + p = Path(name).expanduser() + if not p.is_absolute(): + p = Path("data") / p + if not p.exists(): + raise FileNotFoundError(f"{p} не найден") + x, y = [], [] + with p.open(newline="") as f: + for row in csv.reader(f, delimiter=";"): + if row: + x.append(float(row[0].replace(",", "."))) + y.append(float(row[1].replace(",", "."))) + return x, y + +def _generate_from_function() -> Tuple[List[float], List[float], Callable[[float], float]]: + f = _BUILTIN_FUNCTIONS[input(f"Функция {list(_BUILTIN_FUNCTIONS)}: ").strip().lower()] + a = float(input("a = ")) + b = float(input("b = ")) + n = int(input("n ≥ 2: ")) + h = (b - a) / (n - 1) + xs = [a + i * h for i in range(n)] + ys = [f(x) for x in xs] + return xs, ys, f + +def acquire_dataset(): + print("1 — клавиатура\n2 — csv-файл\n3 — встроенная функция") + m = input("Выбор: ").strip() + if m == "1": + x, y = _read_keyboard() + return x, y, None + if m == "2": + x, y = _read_from_file() + return x, y, None + if m == "3": + return _generate_from_function() + raise ValueError("неверный выбор") + +def ask_target_x() -> float: + return float(input("x* = ")) diff --git "a/\320\2403212/bondarev_408300/lab5/lagrange.py" "b/\320\2403212/bondarev_408300/lab5/lagrange.py" new file mode 100644 index 0000000..f149c1b --- /dev/null +++ "b/\320\2403212/bondarev_408300/lab5/lagrange.py" @@ -0,0 +1,18 @@ +from typing import Sequence +import numpy as np + +def lagrange(x: Sequence[float], y: Sequence[float], x_target: float) -> float: + x = np.asarray(x, dtype=float) + y = np.asarray(y, dtype=float) + idx = np.where(np.isclose(x_target, x, atol=1e-12, rtol=0))[0] + if idx.size: + return float(y[idx[0]]) + n = len(x) + res = 0.0 + for i in range(n): + li = 1.0 + for j in range(n): + if i != j: + li *= (x_target - x[j]) / (x[i] - x[j]) + res += y[i] * li + return float(res) diff --git "a/\320\2403212/bondarev_408300/lab5/main.py" "b/\320\2403212/bondarev_408300/lab5/main.py" new file mode 100644 index 0000000..7c31c5d --- /dev/null +++ "b/\320\2403212/bondarev_408300/lab5/main.py" @@ -0,0 +1,35 @@ +from utils.io_utils import acquire_dataset, ask_target_x +from utils.plot_utils import plot_interpolation +from interpolation.lagrange import lagrange +from interpolation.newton_divided import newton_divided +from interpolation.newton_finite import newton_finite, finite_difference_table + +def main() -> None: + x_nodes, y_nodes, intrinsic_f = acquire_dataset() + + print("\nТаблица конечных разностей:") + fd_table = finite_difference_table(y_nodes) + for row in fd_table: + print(" ".join(f"{v:12.6g}" if v is not None else " " * 12 for v in row)) + + x_target = ask_target_x() + + approx = { + "lagrange": lagrange(x_nodes, y_nodes, x_target), + "newton_divided": newton_divided(x_nodes, y_nodes, x_target), + "newton_finite": newton_finite(x_nodes, y_nodes, x_target), + } + + print(f"\n≈ f({x_target})") + for name, value in approx.items(): + print(f"{name:>15}: {value:.10g}") + + if intrinsic_f is not None: + exact = intrinsic_f(x_target) + print(f"{'exact':>15}: {exact:.10g}") + print(f"{'abs error':>15}: {max(abs(v - exact) for v in approx.values()):.3e}") + + plot_interpolation(x_nodes, y_nodes, intrinsic_f, approx, x_target) + +if __name__ == "__main__": + main() diff --git "a/\320\2403212/bondarev_408300/lab5/newton_divided.py" "b/\320\2403212/bondarev_408300/lab5/newton_divided.py" new file mode 100644 index 0000000..0fde8c1 --- /dev/null +++ "b/\320\2403212/bondarev_408300/lab5/newton_divided.py" @@ -0,0 +1,29 @@ +from typing import Sequence, List +import numpy as np + +def _dup(x: Sequence[float]) -> None: + a = np.asarray(x, dtype=float) + if np.any(np.isclose(np.diff(np.sort(a)), 0.0, atol=1e-12, rtol=0)): + raise ValueError("duplicate x") + +def _table(x: Sequence[float], y: Sequence[float]) -> List[List[float]]: + _dup(x) + n = len(x) + t = [[0.0] * n for _ in range(n)] + for i in range(n): + t[i][0] = y[i] + for j in range(1, n): + for i in range(n - j): + t[i][j] = (t[i + 1][j - 1] - t[i][j - 1]) / (x[i + j] - x[i]) + return t + +def newton_divided(x: Sequence[float], y: Sequence[float], x_target: float) -> float: + x = np.asarray(x, dtype=float) + tab = _table(x, y) + n = len(x) + res = tab[0][0] + prod = 1.0 + for k in range(1, n): + prod *= (x_target - x[k - 1]) + res += tab[0][k] * prod + return float(res) diff --git "a/\320\2403212/bondarev_408300/lab5/newton_finite.py" "b/\320\2403212/bondarev_408300/lab5/newton_finite.py" new file mode 100644 index 0000000..21da2d5 --- /dev/null +++ "b/\320\2403212/bondarev_408300/lab5/newton_finite.py" @@ -0,0 +1,33 @@ +from typing import Sequence, List +import numpy as np + +def _eq(arr: Sequence[float], tol: float = 1e-9) -> bool: + a = np.asarray(arr, dtype=float) + h = a[1] - a[0] + return np.allclose(np.diff(a), h, atol=tol, rtol=0) + +def finite_difference_table(y: Sequence[float]) -> List[List[float]]: + n = len(y) + t: List[List[float]] = [[None] * n for _ in range(n)] + for i, v in enumerate(y): + t[i][0] = float(v) + for j in range(1, n): + for i in range(n - j): + t[i][j] = t[i + 1][j - 1] - t[i][j - 1] + return t + +def newton_finite(x: Sequence[float], y: Sequence[float], x_target: float) -> float: + x = np.asarray(x, dtype=float) + if len(x) < 2 or not _eq(x): + raise ValueError("non-equidistant nodes") + h = x[1] - x[0] + t = (x_target - x[0]) / h + table = finite_difference_table(y) + res = table[0][0] + fact = 1.0 + prod = 1.0 + for k in range(1, len(x)): + fact *= k + prod *= (t - (k - 1)) + res += (prod / fact) * table[0][k] + return float(res) diff --git "a/\320\2403212/bondarev_408300/lab5/plot_utils.py" "b/\320\2403212/bondarev_408300/lab5/plot_utils.py" new file mode 100644 index 0000000..d31445f --- /dev/null +++ "b/\320\2403212/bondarev_408300/lab5/plot_utils.py" @@ -0,0 +1,33 @@ +import os +import matplotlib +if os.environ.get("PYCHARM_HOSTED") == "1": + matplotlib.use("TkAgg") + +import matplotlib.pyplot as plt +import numpy as np +from typing import Sequence, Callable, Mapping +from interpolation.lagrange import lagrange + +def plot_interpolation( + x_nodes: Sequence[float], + y_nodes: Sequence[float], + intrinsic_f: Callable[[float], float] | None, + approx_f_values: Mapping[str, float], + x_target: float +) -> None: + xs = np.linspace(min(x_nodes), max(x_nodes), 400) + + plt.figure() + if intrinsic_f is not None: + plt.plot(xs, [intrinsic_f(x) for x in xs], label="f(x)") + + plt.plot(xs, [lagrange(x_nodes, y_nodes, xi) for xi in xs], label="interpolant") + plt.scatter(x_nodes, y_nodes, color="black", zorder=5, label="nodes") + + for name, y_val in approx_f_values.items(): + plt.scatter([x_target], [y_val], label=f"{name} @ {x_target:.3g}") + + plt.title("Интерполяция") + plt.legend() + plt.grid(True) + plt.show() diff --git "a/\320\2403212/bondarev_408300/lab5/variant1.csv" "b/\320\2403212/bondarev_408300/lab5/variant1.csv" new file mode 100644 index 0000000..2a89477 --- /dev/null +++ "b/\320\2403212/bondarev_408300/lab5/variant1.csv" @@ -0,0 +1,7 @@ +0.25;1.2557 +0.30;2.1764 +0.35;3.1218 +0.40;4.0482 +0.45;5.9875 +0.50;6.9195 +0.55;7.8359 diff --git "a/\320\2403212/bondarev_408300/lab5/variant2.csv" "b/\320\2403212/bondarev_408300/lab5/variant2.csv" new file mode 100644 index 0000000..e49f060 --- /dev/null +++ "b/\320\2403212/bondarev_408300/lab5/variant2.csv" @@ -0,0 +1,5 @@ +-1.0;0.00000 +-0.5;0.47943 + 0.0;1.00000 + 0.5;1.64872 + 1.0;2.71828 diff --git "a/\320\2403212/bondarev_408300/lab5/variant3.csv" "b/\320\2403212/bondarev_408300/lab5/variant3.csv" new file mode 100644 index 0000000..9275c85 --- /dev/null +++ "b/\320\2403212/bondarev_408300/lab5/variant3.csv" @@ -0,0 +1,8 @@ +0.0;1.0 +0.2;0.9800665778 +0.4;0.9210609940 +0.6;0.8253356149 +0.8;0.6967067093 +1.0;0.5403023059 +1.2;0.3623577545 +1.4;0.1699671429 diff --git "a/\320\2403212/bondarev_408300/lab6/CalMathLab6 BondarevAM.pdf" "b/\320\2403212/bondarev_408300/lab6/CalMathLab6 BondarevAM.pdf" new file mode 100644 index 0000000..e51f787 Binary files /dev/null and "b/\320\2403212/bondarev_408300/lab6/CalMathLab6 BondarevAM.pdf" differ diff --git "a/\320\2403212/bondarev_408300/lab6/equations.py" "b/\320\2403212/bondarev_408300/lab6/equations.py" new file mode 100644 index 0000000..b1421e1 --- /dev/null +++ "b/\320\2403212/bondarev_408300/lab6/equations.py" @@ -0,0 +1,25 @@ +import math + +def f1(x, y): + return x + y + +def exact1(x, y0, x0): + return (y0 + x0 + 1) * math.exp(x - x0) - x - 1 + +def f2(x, y): + return y - x ** 2 + 1 + +def exact2(x, y0, x0): + return (y0 - (x0 ** 2 - 2 * x0 + 2)) * math.exp(x - x0) + x ** 2 - 2 * x + 2 + +def f3(x, y): + return x * y + +def exact3(x, y0, x0): + return y0 * math.exp((x ** 2 - x0 ** 2) / 2) + +ODES = { + 1: ("y' = x + y", f1, exact1, 1), + 2: ("y' = y - x^2 + 1", f2, exact2, 1), + 3: ("y' = x·y", f3, exact3, 1), +} diff --git "a/\320\2403212/bondarev_408300/lab6/main.py" "b/\320\2403212/bondarev_408300/lab6/main.py" new file mode 100644 index 0000000..bccaba4 --- /dev/null +++ "b/\320\2403212/bondarev_408300/lab6/main.py" @@ -0,0 +1,41 @@ +import os +os.environ["MPLBACKEND"] = "TkAgg" +import matplotlib.pyplot as plt +from equations import ODES +from ode_methods import euler, rk4, adams_pc4, adaptive +from utils import print_table + +def main(): + print("Выберите одно из уравнений:") + for k, v in ODES.items(): + print(f"{k}: {v[0]}") + eq = int(input("Номер уравнения: ")) + f = ODES[eq][1] + exact = ODES[eq][2] + x0 = float(input("x0: ")) + xn = float(input("xn: ")) + y0 = float(input("y0: ")) + h0 = float(input("Начальный шаг h: ")) + eps = float(input("Точность ε для одношаговых методов: ")) + xs_e, ys_e, h_e, err_e = adaptive(euler, 1, f, x0, y0, xn, h0, eps) + xs_rk, ys_rk, h_rk, err_rk = adaptive(rk4, 4, f, x0, y0, xn, h0, eps) + xs_ad, ys_ad = adams_pc4(f, x0, y0, xn, h_rk) + err_ad = max(abs(ys_ad[i] - exact(xs_ad[i], y0, x0)) for i in range(len(xs_ad))) + print(f"\nМетод Эйлера: h = {h_e} |ε| ≈ {err_e}") + print(f"Рунге-Кутта 4: h = {h_rk} |ε| ≈ {err_rk}") + print(f"Адамса : h = {h_rk} |ε| ≈ {err_ad}\n") + print_table(xs_ad, ys_e, ys_rk, ys_ad, exact, y0, x0) + plt.plot(xs_e, ys_e, label="Euler") + plt.plot(xs_rk, ys_rk, label="RK4") + plt.plot(xs_ad, ys_ad, label="Adams") + xs_dense = [x0 + i * (xn - x0) / 1000 for i in range(1001)] + ys_dense = [exact(x, y0, x0) for x in xs_dense] + plt.plot(xs_dense, ys_dense, label="Exact") + plt.legend() + plt.xlabel("x") + plt.ylabel("y") + plt.grid(True) + plt.show() + +if __name__ == "__main__": + main() diff --git "a/\320\2403212/bondarev_408300/lab6/ode_methods.py" "b/\320\2403212/bondarev_408300/lab6/ode_methods.py" new file mode 100644 index 0000000..02f97b9 --- /dev/null +++ "b/\320\2403212/bondarev_408300/lab6/ode_methods.py" @@ -0,0 +1,67 @@ +import math + +def euler(f, x0, y0, xn, h): + xs = [x0] + ys = [y0] + x = x0 + y = y0 + while x < xn - 1e-12: + y += h * f(x, y) + x += h + xs.append(x) + ys.append(y) + return xs, ys + +def rk4(f, x0, y0, xn, h): + xs = [x0] + ys = [y0] + x = x0 + y = y0 + while x < xn - 1e-12: + k1 = f(x, y) + k2 = f(x + h / 2, y + h * k1 / 2) + k3 = f(x + h / 2, y + h * k2 / 2) + k4 = f(x + h, y + h * k3) + y += h * (k1 + 2 * k2 + 2 * k3 + k4) / 6 + x += h + xs.append(x) + ys.append(y) + return xs, ys + +def adams_pc4(f, x0, y0, xn, h): + xs, ys = rk4(f, x0, y0, x0 + 3 * h, h) + x = xs[-1] + while x < xn - 1e-12: + i = len(xs) - 1 + f0 = f(xs[i], ys[i]) + f1 = f(xs[i - 1], ys[i - 1]) + f2 = f(xs[i - 2], ys[i - 2]) + f3 = f(xs[i - 3], ys[i - 3]) + y_pred = ys[i] + h * (55 * f0 - 59 * f1 + 37 * f2 - 9 * f3) / 24 + x_next = x + h + f_pred = f(x_next, y_pred) + y_corr = ys[i] + h * (9 * f_pred + 19 * f0 - 5 * f1 + f2) / 24 + x = x_next + xs.append(x) + ys.append(y_corr) + return xs, ys + +def runge_error(method, p, f, x0, y0, xn, h): + _, y1 = method(f, x0, y0, xn, h) + _, y2 = method(f, x0, y0, xn, h / 2) + return abs(y2[-1] - y1[-1]) / (2 ** p - 1) + +def adaptive(method, p, f, x0, y0, xn, h, eps, + max_halvings=20, max_steps=200_000): + for _ in range(max_halvings): + err = runge_error(method, p, f, x0, y0, xn, h) + steps = int((xn - x0) / h) + 1 + if err <= eps or steps > max_steps: + xs, ys = method(f, x0, y0, xn, h) + return xs, ys, h, err + h /= 2 + print(f"⚠️ adaptive: точность {eps} не достигнута после {max_halvings} делений; " + f"err ≈ {err:e}, h = {h:e}") + xs, ys = method(f, x0, y0, xn, h) + return xs, ys, h, err + diff --git "a/\320\2403212/bondarev_408300/lab6/utils.py" "b/\320\2403212/bondarev_408300/lab6/utils.py" new file mode 100644 index 0000000..f7f10eb --- /dev/null +++ "b/\320\2403212/bondarev_408300/lab6/utils.py" @@ -0,0 +1,10 @@ +def print_table(xs, ys_e, ys_rk, ys_ad, exact_func, y0, x0): + header = f"{'x':>10}{'Euler':>15}{'RK4':>15}{'Adams':>15}{'Exact':>15}" + print(header) + for i in range(len(xs)): + x = xs[i] + y_exact = exact_func(x, y0, x0) + y_e = ys_e[i] if i < len(ys_e) else float('nan') + y_rk = ys_rk[i] if i < len(ys_rk) else float('nan') + y_ad = ys_ad[i] + print(f"{x:10.6f}{y_e:15.6f}{y_rk:15.6f}{y_ad:15.6f}{y_exact:15.6f}")