upd

archive/labs/lab1/code/main.py

@@ -0,0 +1,189 @@
+import matplotlib.pyplot as plt
+import math
+import numpy as np
+from scipy.integrate import solve_ivp
+from scipy.optimize import fsolve
+
+
+class ODESolver:
+    def euler_method(self, f, t_span, y0, h):
+        t0, t_end = t_span
+        N = int((t_end - t0) / h)
+        t = np.linspace(t0, t_end, N + 1)
+        y = np.zeros(N + 1)
+        y[0] = y0
+
+        for i in range(N):
+            y[i + 1] = y[i] + h * f(t[i], y[i])
+
+        return t, y
+
+    def explicit_trapezoid(self, f, t_span, y0, h):
+        t0, t_end = t_span
+        N = int((t_end - t0) / h)
+        t = np.linspace(t0, t_end, N + 1)
+        y = np.zeros(N + 1)
+        y[0] = y0
+
+        for i in range(N):
+            k1 = f(t[i], y[i])
+            k2 = f(t[i] + h, y[i] + h * k1)
+            y[i + 1] = y[i] + (h / 2) * (k1 + k2)
+
+        return t, y
+
+    def rk4_method(self, f, t_span, y0, h):
+        t0, t_end = t_span
+        N = int((t_end - t0) / h)
+        t = np.linspace(t0, t_end, N + 1)
+        y = np.zeros(N + 1)
+        y[0] = y0
+
+        for i in range(N):
+            k1 = f(t[i], y[i])
+            k2 = f(t[i] + h/2, y[i] + h/2 * k1)
+            k3 = f(t[i] + h/2, y[i] + h/2 * k2)
+            k4 = f(t[i] + h, y[i] + h * k3)
+            y[i + 1] = y[i] + h/6 * (k1 + 2*k2 + 2*k3 + k4)
+
+        return t, y
+
+    def backward_euler(self, f, t_span, y0, h):
+        t0, t_end = t_span
+        N = int((t_end - t0) / h)
+        t = np.linspace(t0, t_end, N + 1)
+        y = np.zeros(N + 1)
+        y[0] = y0
+
+        for i in range(N):
+            t_next = t[i + 1]
+
+            def equation(z):
+                return z - h * f(t_next, z) - y[i]
+
+            y[i + 1] = fsolve(equation, y[i])[0]
+
+        return t, y
+
+    def implicit_trapezoid(self, f, t_span, y0, h):
+        t0, t_end = t_span
+        N = int((t_end - t0) / h)
+        t = np.linspace(t0, t_end, N + 1)
+        y = np.zeros(N + 1)
+        y[0] = y0
+
+        for i in range(N):
+            t_n = t[i]
+            t_np1 = t[i + 1]
+            y_n = y[i]
+
+            def equation(z):
+                return z - y_n - (h/2) * (f(t_n, y_n) + f(t_np1, z))
+
+            y[i + 1] = fsolve(equation, y_n)[0]
+
+        return t, y
+
+def problem_1_4():
+    print("=== ЗАДАЧА 1.4: y' = y² - t, y(0) = 1 ===")
+
+    def f(t, y):
+        return y**2 - t
+
+    t_span = [0, 1.0]
+    y0 = 1
+    h = 0.25
+
+    solver = ODESolver()
+
+    t_euler, y_euler = solver.euler_method(f, t_span, y0, h)
+    t_trap, y_trap = solver.explicit_trapezoid(f, t_span, y0, h)
+    t_rk4, y_rk4 = solver.rk4_method(f, t_span, y0, h)
+
+    sol_ref = solve_ivp(f, t_span, [y0], method='RK45',
+                       rtol=1e-10, atol=1e-12, dense_output=True)
+    t_ref = np.linspace(t_span[0], t_span[1], 100)
+    y_ref = sol_ref.sol(t_ref)[0]
+
+    y_ref_euler = sol_ref.sol(t_euler)[0]
+    y_ref_trap = sol_ref.sol(t_trap)[0]
+    y_ref_rk4 = sol_ref.sol(t_rk4)[0]
+
+    err_euler = np.linalg.norm(y_euler - y_ref_euler)
+    err_trap = np.linalg.norm(y_trap - y_ref_trap)
+    err_rk4 = np.linalg.norm(y_rk4 - y_ref_rk4)
+
+    print(f"Ошибка метода Эйлера: {err_euler:.2e}")
+    print(f"Ошибка явного метода трапеций: {err_trap:.2e}")
+    print(f"Ошибка метода РК4: {err_rk4:.2e}")
+
+    plt.figure(figsize=(12, 8))
+    plt.plot(t_ref, y_ref, 'k-', linewidth=5, label='Эталонное решение (solve_ivp)')
+    plt.plot(t_euler, y_euler, 'bo-', markersize=10, label='Метод Эйлера')
+    plt.plot(t_trap, y_trap, 'c*-',  markersize=15, label='Явный метод трапеций')
+    plt.plot(t_rk4, y_rk4, 'rs-', markersize=5, label='Метод РК4')
+
+    plt.xlabel('t')
+    plt.ylabel('y(t)')
+    plt.title('Задача 1.4: y\' = y² - t, y(0) = 1')
+    plt.legend()
+    plt.grid(True, alpha=0.3)
+    plt.show()
+
+    return t_euler, y_euler, t_trap, y_trap, t_rk4, y_rk4, t_ref, y_ref
+
+
+if __name__ == "__main__":
+    result1 = problem_1_4()
+
+def problem_2_4():
+    print("\n=== ЗАДАЧА 2.4: y' = -2000y + t, y(0) = 0 ===")
+
+    def f(t, y):
+        return -2000 * y + t
+
+    t_span = [0, 1]
+    y0 = 0
+    h = 1
+
+    solver = ODESolver()
+
+    t_euler, y_euler = solver.euler_method(f, t_span, y0, h)
+    t_beuler, y_beuler = solver.backward_euler(f, t_span, y0, h)
+    t_itrap, y_itrap = solver.implicit_trapezoid(f, t_span, y0, h)
+
+    sol_ref = solve_ivp(f, t_span, [y0], method='Radau',
+                       rtol=1e-10, atol=1e-12, dense_output=True)
+    t_ref = np.linspace(t_span[0], t_span[1], 100)
+    y_ref = sol_ref.sol(t_ref)[0]
+
+    y_ref_euler = sol_ref.sol(t_euler)[0]
+    y_ref_beuler = sol_ref.sol(t_beuler)[0]
+    y_ref_itrap = sol_ref.sol(t_itrap)[0]
+
+    err_euler = np.linalg.norm(y_euler - y_ref_euler)
+    err_beuler = np.linalg.norm(y_beuler - y_ref_beuler)
+    err_itrap = np.linalg.norm(y_itrap - y_ref_itrap)
+
+    print(f"Ошибка явного метода Эйлера: {err_euler:.2e}")
+    print(f"Ошибка неявного метода Эйлера: {err_beuler:.2e}")
+    print(f"Ошибка неявного метода трапеций: {err_itrap:.2e}")
+
+    plt.figure(figsize=(12, 8))
+    plt.plot(t_ref, y_ref, 'k-', linewidth=5, label='Эталонное решение (solve_ivp)')
+    plt.plot(t_euler, y_euler, 'bo-', markersize=4, label='Явный метод Эйлера')
+    plt.plot(t_beuler, y_beuler, 'g^-', markersize=4, label='Неявный метод Эйлера')
+    plt.plot(t_itrap, y_itrap, 'm*-',linewidth=2.2, markersize=4, label='Неявный метод трапеций')
+
+    plt.xlabel('t')
+    plt.ylabel('y(t)')
+    plt.title('Задача 2.4: y\' = -2000y + t, y(0) = 0 (жесткая система)')
+    plt.legend()
+    plt.grid(True, alpha=0.3)
+    plt.show()
+
+    return t_euler, y_euler, t_beuler, y_beuler, t_itrap, y_itrap, t_ref, y_ref
+
+
+if __name__ == "__main__":
+    result2 = problem_2_4()

archive/labs/lab1/report/res.typ

@@ -0,0 +1,643 @@
+#set text(size: 1.3em)
+#show link: underline
+#set page(footer: context {
+  if counter(page).get().first() > 1 [
+    #align(center)[
+      #counter(page).display("1")
+    ]
+  ]
+  if counter(page).get().first() == 1 [
+    #align(center)[
+      Санкт-Петербург \ 2025
+    ]
+  ]
+})
+
+#set page(header: context {
+  if counter(page).get().first() == 1 [
+    #align(center)[
+      Санкт-Петербургский национальный исследовательский университет информационных технологий, механики и оптики
+    ]
+  ]
+})
+
+
+
+#show raw.where(block: false): box.with(
+  fill: luma(240),
+  inset: (x: 3pt, y: 0pt),
+  outset: (y: 3pt),
+  radius: 2pt,
+)
+
+#show raw.where(block: true): block.with(
+  fill: luma(240),
+  inset: 10pt,
+  radius: 4pt,
+)
+
+// title
+
+#for _ in range(5) { linebreak() }
+
+#align(center)[Расчетно графическая работа №1]
+
+#for _ in range(15) { linebreak() }
+
+
+#align(right)[Выполнили:]
+#align(right)[Левахин Лев]
+#align(right)[Останин Андрей]
+#align(right)[Дощенников Никита]
+#align(right)[Группы: К3221, К3240]
+#align(right)[Проверил:]
+#align(right)[Владимир Владимирович Беспалов]
+
+#pagebreak()
+
+=== Цель
+
+Изучение и практическое применение численных методов решения обыкновенных дифференциальных уравнений, сравнение их точности и устойчивости на различных типах задач.
+
+=== Задачи
+
+1. Реализовать следующие численные методы решения ОДУ:
+
+    - Явный метод Эйлера (1-го порядка)
+    - Явный метод трапеций (улучшенный Эйлер, 2-го порядка)
+    - Классический метод Рунге-Кутты 4-го порядка (RK4)
+    - Неявный метод Эйлера (1-го порядка)
+    - Неявный метод трапеций (2-го порядка)
+
+
+2. Решить две задачи:
+
+    - Нежесткое нелинейное уравнение типа Риккати (`1.4`): исследование точности и сходимости явных методов
+    - Жесткое линейное уравнение (`2.4`): исследование устойчивости неявных методов
+
+
+3. Провести анализ погрешностей и построить графики решений
+4. Сформулировать выводы о применимости различных методов
+
+=== Теоретические основы
+
+==== Постановка задачи Коши
+
+Рассматривается начальная задача для обыкновенного дифференциального уравнения первого порядка:
+
+$
+frac(d y, d x) eq f(t, y), space.quad y(t_0) eq y_0
+$
+
+где
+
+- $t$ - независимая переменная (время)
+- $y(t)$ - искомая функция
+- $f(t, y)$ - правая часть уравнения
+- $y_0$ - начальное условие
+
+==== Классификация численных методов
+
+*Явные методы* - значение решения на следующем шаге вычисляется непосредственно из известных значений:
+
+- Просты в реализации
+- Требуют малого шага для жестких систем
+- Примеры: явный Эйлер, RK4
+
+*Неявные методы* - значение решения на следующем шаге находится из неявного уравнения:
+
+- Требуют решения нелинейных уравнений на каждом шаге
+- Обладают лучшей устойчивостью для жестких систем
+- Примеры: неявный Эйлер, неявная трапеция
+
+=== Описание реализованных методов
+
+==== Явный метод Эйлера
+
+$
+y_(n + 1) eq y_n plus h dot f(t_n, y_n)
+$
+
+Имеет первый порядок точности $O(h)$.
+
+```python
+def euler_method(self, f, t_span, y0, h):
+    t0, t_end = t_span
+    N = int((t_end - t0) / h)
+    t = np.linspace(t0, t_end, N + 1)
+    y = np.zeros(N + 1)
+    y[0] = y0
+
+    for i in range(N):
+        y[i + 1] = y[i] + h * f(t[i], y[i])
+
+    return t, y
+```
+
+Метод использует линейную аппроксимацию решения на основе производной в начале интервала. Это самый простой, но наименее точный метод.
+
+==== Явный метод трапеций
+
+$
+k_1 eq f(t_n, y_n) \
+k_2 eq f(t_n plus h, y_n plus h dot k_1) \
+y_(n plus 1) eq y_n plus frac(h, 2) dot (k_1 plus k_2)
+$
+
+Имеет второй порядок точности $O(h^2)$.
+
+```python
+def explicit_trapezoid(self, f, t_span, y0, h):
+    t0, t_end = t_span
+    N = int((t_end - t0) / h)
+    t = np.linspace(t0, t_end, N + 1)
+    y = np.zeros(N + 1)
+    y[0] = y0
+
+    for i in range(N):
+        k1 = f(t[i], y[i])
+        k2 = f(t[i] + h, y[i] + h * k1)
+        y[i + 1] = y[i] + (h / 2) * (k1 + k2)
+
+    return t, y
+```
+
+Метод сначала делает предсказание решения в конце интервала (явным Эйлером), затем усредняет производные в начале и конце интервала.
+
+==== Метод Рунге-Кутты 4-го порядка (RK4)
+
+$
+k_1 eq f(t_n, y_n) \
+k_2 eq f(t_n plus frac(h, 2), y_n plus h dot frac(k_1, 2)) \
+k_3 eq f(t_n plus frac(h, 2), y_n plus h dot frac(k_2, 2)) \
+k_4 eq f(t_n plus h, y_n plus h dot k_3) \
+y_(n plus 1) = y_n + frac(h, 6) dot (k_1 plus 2 k_2 plus 2 k_3 plus k_4)
+$
+
+Имеет четвертый порядок точности $O(h^4)$
+
+```python
+def rk4_method(self, f, t_span, y0, h):
+    t0, t_end = t_span
+    N = int((t_end - t0) / h)
+    t = np.linspace(t0, t_end, N + 1)
+    y = np.zeros(N + 1)
+    y[0] = y0
+
+    for i in range(N):
+        k1 = f(t[i], y[i])
+        k2 = f(t[i] + h/2, y[i] + h/2 * k1)
+        k3 = f(t[i] + h/2, y[i] + h/2 * k2)
+        k4 = f(t[i] + h, y[i] + h * k3)
+        y[i + 1] = y[i] + h/6 * (k1 + 2*k2 + 2*k3 + k4)
+
+    return t, y
+```
+
+Метод вычисляет четыре промежуточных значения производной внутри интервала и формирует взвешенное среднее. Это обеспечивает высокую точность.
+
+==== Неявный метод Эйлера
+
+$
+y_(n plus 1) eq y_n plus h dot f(t_(n plus 1), y_(n plus 1))
+$
+
+Имеет первый порядок точности $O(h)$
+
+```python
+def backward_euler(self, f, t_span, y0, h):
+    t0, t_end = t_span
+    N = int((t_end - t0) / h)
+    t = np.linspace(t0, t_end, N + 1)
+    y = np.zeros(N + 1)
+    y[0] = y0
+
+    for i in range(N):
+        t_next = t[i + 1]
+
+        def equation(z):
+            return z - h * f(t_next, z) - y[i]
+
+        y[i + 1] = fsolve(equation, y[i])[0]
+
+    return t, y
+```
+
+Метод использует производную в конце интервала, что требует решения нелинейного уравнения на каждом шаге с помощью `fsolve`. Это обеспечивает безусловную устойчивость для линейных задач.
+
+==== Неявный метод трапеций
+
+$
+y_(n plus 1) eq y_n plus frac(h, 2) dot [f(t_n, y_n) plus f(t_(n plus 1), t_(n + 1)]
+$
+
+Имеет второй порядок точности $O(h^2)$
+
+```python
+def implicit_trapezoid(self, f, t_span, y0, h):
+    t0, t_end = t_span
+    N = int((t_end - t0) / h)
+    t = np.linspace(t0, t_end, N + 1)
+    y = np.zeros(N + 1)
+    y[0] = y0
+
+    for i in range(N):
+        t_n = t[i]
+        t_np1 = t[i + 1]
+        y_n = y[i]
+
+        def equation(z):
+            return z - y_n - (h/2) * (f(t_n, y_n) + f(t_np1, z))
+
+        y[i + 1] = fsolve(equation, y_n)[0]
+
+    return t, y
+```
+
+Метод усредняет производные в начале и конце интервала (как явная трапеция), но использует неизвестное значение $y_(n+1)$, что требует решения нелинейного уравнения. Метод A-устойчив.
+
+=== Задачи
+
+==== 1.4
+
+$
+y' eq y^2 - t, space.quad y(0) eq 1, space.quad t in [0, 1.5], space.quad h eq 0.25
+$
+
+- Нелинейное уравнение типа Риккати
+- Не имеет аналитического решения
+- Нежесткая система - подходит для явных методов
+- Используется для анализа точности и сходимости
+
+```python
+def problem_1_4():
+    def f(t, y):
+        return y**2 - t
+
+    t_span = [0, 1.5]
+    y0 = 1
+    h = 0.25
+
+    solver = ODESolver()
+
+    t_euler, y_euler = solver.euler_method(f, t_span, y0, h)
+    t_trap, y_trap = solver.explicit_trapezoid(f, t_span, y0, h)
+    t_rk4, y_rk4 = solver.rk4_method(f, t_span, y0, h)
+
+    sol_ref = solve_ivp(f, t_span, [y0], method='RK45',
+                       rtol=1e-10, atol=1e-12, dense_output=True)
+```
+
+==== 2.4
+
+$
+y' eq -2000y plus t, space.quad y(0) eq 0, space.quad t in [0, 1], space.quad h eq 1
+$
+
+- Линейное уравнение с большим отрицательным коэффициентом $lambda = -2000$
+- Жесткая система
+- Явный метод Эйлера неустойчив при больших шагах
+- Требует неявных методов для устойчивости
+
+```python
+def problem_2_4():
+    def f(t, y):
+        return -2000 * y + t
+
+    t_span = [0, 1]
+    y0 = 0
+    h = 1
+
+    solver = ODESolver()
+
+    t_euler, y_euler = solver.euler_method(f, t_span, y0, h)
+    t_beuler, y_beuler = solver.backward_euler(f, t_span, y0, h)
+    t_itrap, y_itrap = solver.implicit_trapezoid(f, t_span, y0, h)
+
+    sol_ref = solve_ivp(f, t_span, [y0], method='Radau',
+                       rtol=1e-10, atol=1e-12, dense_output=True)
+```
+
+=== Результаты
+
+==== 1.4
+
+Мы уменьшили интервал с $[0, 1.5]$ до $[0, 1]$ для того, чтобы избежать "взрыва решения", так как уравнение нелинейно.
+
+Результаты программы:
+
+```bash
+=== ЗАДАЧА 1.4: y' = y² - t, y(0) = 1 ===
+Ошибка метода Эйлера: 6.44e+00
+Ошибка явного метода трапеций: 3.76e+00
+Ошибка метода РК4: 6.45e-01
+```
+
+Ошибка получилась порядка $10^(-1)$.
+
+#align(center)[
+  #figure(
+    image("assets/1.png"),
+    supplement: [Рис.],
+    caption: [График решения задачи 1.4]
+  ) <g1>
+]
+
+- RK4 самый точный. Ошибка $tilde 10^(-1)$ в разы меньше остальных. На графике траектория почти совпадает с эталоном.
+
+- Явная трапеция имеет среднюю точность. Ошибка меньше, чем у Эйлера, но заметно хуже RK4.
+
+- Метод эйлера самый неточный. Ошибка огромная, отклонение от эталона быстро растёт.
+
+==== 2.4
+
+Результаты программы:
+
+```bash
+=== ЗАДАЧА 2.4: y' = -2000y + t, y(0) = 0 ===
+Ошибка явного метода Эйлера: 5.00e-04
+Ошибка неявного метода Эйлера: 1.25e-10
+Ошибка неявного метода трапеций: 2.50e-07
+```
+
+
+#align(center)[
+  #figure(
+    image("assets/2.png"),
+    supplement: [Рис.],
+    caption: [График решения задачи 2.4]
+  ) <g2>
+]
+
+- Неявный Эйлер самый точный и устойчивый. Ошибка $tilde 10^(-10)$.
+
+- Неявная трапеция тоже устойчива, но точность хуже. Ошибка $tilde 10^(-7)$.
+
+- Явный Эйлер нестабилен, даёт неправильную форму решения.
+
+=== Структура программы
+
+==== Класс `ODESolver`
+
+Класс инкапсулирует все численные методы:
+
+```python
+class ODESolver:
+    def euler_method(self, f, t_span, y0, h): ...
+    def explicit_trapezoid(self, f, t_span, y0, h): ...
+    def rk4_method(self, f, t_span, y0, h): ...
+    def backward_euler(self, f, t_span, y0, h): ...
+    def implicit_trapezoid(self, f, t_span, y0, h): ...
+```
+
+==== Функция решения задач
+
+Каждая задача решается в отдельной функции:
+
+- `problem_1_4()` - для нежесткой задачи
+- `problem_2_4()` - для жесткой задачи
+
+Функции выполняют:
+
+1. Определение правой части ОДУ
+2. Задание параметров интегрирования
+3. Вызов численных методов
+4. Получение эталонного решения
+5. Вычисление погрешностей
+6. Построение графиков
+
+==== Вычисление эталонного решения
+
+Используются высокоточные методы из библиотеки SciPy:
+
+- `solve_ivp` с методом RK45 для нежестких систем
+- `solve_ivp` с методом Radau для жестких систем
+- Строгие допуски: `rtol=1e-10`, `atol=1e-12`
+
+==== Анализ погрешности
+
+Глобальная погрешность вычисляется как норма разности:
+
+```python
+err = np.linalg.norm(y_numerical - y_reference)
+```
+
+где используется евклидова норма $L_2$.
+
+=== Выводы
+
+О порядке точности:
+
+- Порядок метода существенно влияет на точность при одинаковом шаге
+- Для достижения заданной точности метод более высокого порядка требует меньше вычислений
+- RK4 является оптимальным выбором для нежестких задач
+
+
+О явных и неявных методах:
+
+- Явные методы просты в реализации, но имеют ограничения по устойчивости
+- Неявные методы требуют решения нелинейных уравнений, но обеспечивают устойчивость
+- Для жестких систем неявные методы необходимы
+
+
+О жесткости:
+
+- Жесткие системы характеризуются наличием сильно различающихся масштабов времени
+- Явные методы требуют непрактично малых шагов для жестких систем
+- A-устойчивые неявные методы позволяют использовать большие шаги
+
+=== Приложение
+
+Полный код программы:
+
+```python
+import matplotlib.pyplot as plt
+import math
+import numpy as np
+from scipy.integrate import solve_ivp
+from scipy.optimize import fsolve
+
+
+class ODESolver:
+    def euler_method(self, f, t_span, y0, h):
+        t0, t_end = t_span
+        N = int((t_end - t0) / h)
+        t = np.linspace(t0, t_end, N + 1)
+        y = np.zeros(N + 1)
+        y[0] = y0
+
+        for i in range(N):
+            y[i + 1] = y[i] + h * f(t[i], y[i])
+
+        return t, y
+
+    def explicit_trapezoid(self, f, t_span, y0, h):
+        t0, t_end = t_span
+        N = int((t_end - t0) / h)
+        t = np.linspace(t0, t_end, N + 1)
+        y = np.zeros(N + 1)
+        y[0] = y0
+
+        for i in range(N):
+            k1 = f(t[i], y[i])
+            k2 = f(t[i] + h, y[i] + h * k1)
+            y[i + 1] = y[i] + (h / 2) * (k1 + k2)
+
+        return t, y
+
+    def rk4_method(self, f, t_span, y0, h):
+        t0, t_end = t_span
+        N = int((t_end - t0) / h)
+        t = np.linspace(t0, t_end, N + 1)
+        y = np.zeros(N + 1)
+        y[0] = y0
+
+        for i in range(N):
+            k1 = f(t[i], y[i])
+            k2 = f(t[i] + h/2, y[i] + h/2 * k1)
+            k3 = f(t[i] + h/2, y[i] + h/2 * k2)
+            k4 = f(t[i] + h, y[i] + h * k3)
+            y[i + 1] = y[i] + h/6 * (k1 + 2*k2 + 2*k3 + k4)
+
+        return t, y
+
+    def backward_euler(self, f, t_span, y0, h):
+        t0, t_end = t_span
+        N = int((t_end - t0) / h)
+        t = np.linspace(t0, t_end, N + 1)
+        y = np.zeros(N + 1)
+        y[0] = y0
+
+        for i in range(N):
+            t_next = t[i + 1]
+
+            def equation(z):
+                return z - h * f(t_next, z) - y[i]
+
+            y[i + 1] = fsolve(equation, y[i])[0]
+
+        return t, y
+
+    def implicit_trapezoid(self, f, t_span, y0, h):
+        t0, t_end = t_span
+        N = int((t_end - t0) / h)
+        t = np.linspace(t0, t_end, N + 1)
+        y = np.zeros(N + 1)
+        y[0] = y0
+
+        for i in range(N):
+            t_n = t[i]
+            t_np1 = t[i + 1]
+            y_n = y[i]
+
+            def equation(z):
+                return z - y_n - (h/2) * (f(t_n, y_n) + f(t_np1, z))
+
+            y[i + 1] = fsolve(equation, y_n)[0]
+
+        return t, y
+
+def problem_1_4():
+    print("=== ЗАДАЧА 1.4: y' = y² - t, y(0) = 1 ===")
+
+    def f(t, y):
+        return y**2 - t
+
+    t_span = [0, 1.0]
+    y0 = 1
+    h = 0.25
+
+    solver = ODESolver()
+
+    t_euler, y_euler = solver.euler_method(f, t_span, y0, h)
+    t_trap, y_trap = solver.explicit_trapezoid(f, t_span, y0, h)
+    t_rk4, y_rk4 = solver.rk4_method(f, t_span, y0, h)
+
+    sol_ref = solve_ivp(f, t_span, [y0], method='RK45',
+                       rtol=1e-10, atol=1e-12, dense_output=True)
+    t_ref = np.linspace(t_span[0], t_span[1], 100)
+    y_ref = sol_ref.sol(t_ref)[0]
+
+    y_ref_euler = sol_ref.sol(t_euler)[0]
+    y_ref_trap = sol_ref.sol(t_trap)[0]
+    y_ref_rk4 = sol_ref.sol(t_rk4)[0]
+
+    err_euler = np.linalg.norm(y_euler - y_ref_euler)
+    err_trap = np.linalg.norm(y_trap - y_ref_trap)
+    err_rk4 = np.linalg.norm(y_rk4 - y_ref_rk4)
+
+    print(f"Ошибка метода Эйлера: {err_euler:.2e}")
+    print(f"Ошибка явного метода трапеций: {err_trap:.2e}")
+    print(f"Ошибка метода РК4: {err_rk4:.2e}")
+
+    plt.figure(figsize=(12, 8))
+    plt.plot(t_ref, y_ref, 'k-', linewidth=5, label='Эталонное решение (solve_ivp)')
+    plt.plot(t_euler, y_euler, 'bo-', markersize=10, label='Метод Эйлера')
+    plt.plot(t_trap, y_trap, 'c*-',  markersize=15, label='Явный метод трапеций')
+    plt.plot(t_rk4, y_rk4, 'rs-', markersize=5, label='Метод РК4')
+
+    plt.xlabel('t')
+    plt.ylabel('y(t)')
+    plt.title('Задача 1.4: y\' = y² - t, y(0) = 1')
+    plt.legend()
+    plt.grid(True, alpha=0.3)
+    plt.show()
+
+    return t_euler, y_euler, t_trap, y_trap, t_rk4, y_rk4, t_ref, y_ref
+
+
+if __name__ == "__main__":
+    result1 = problem_1_4()
+
+def problem_2_4():
+    print("\n=== ЗАДАЧА 2.4: y' = -2000y + t, y(0) = 0 ===")
+
+    def f(t, y):
+        return -2000 * y + t
+
+    t_span = [0, 1]
+    y0 = 0
+    h = 1
+
+    solver = ODESolver()
+
+    t_euler, y_euler = solver.euler_method(f, t_span, y0, h)
+    t_beuler, y_beuler = solver.backward_euler(f, t_span, y0, h)
+    t_itrap, y_itrap = solver.implicit_trapezoid(f, t_span, y0, h)
+
+    sol_ref = solve_ivp(f, t_span, [y0], method='Radau',
+                       rtol=1e-10, atol=1e-12, dense_output=True)
+    t_ref = np.linspace(t_span[0], t_span[1], 100)
+    y_ref = sol_ref.sol(t_ref)[0]
+
+    y_ref_euler = sol_ref.sol(t_euler)[0]
+    y_ref_beuler = sol_ref.sol(t_beuler)[0]
+    y_ref_itrap = sol_ref.sol(t_itrap)[0]
+
+    err_euler = np.linalg.norm(y_euler - y_ref_euler)
+    err_beuler = np.linalg.norm(y_beuler - y_ref_beuler)
+    err_itrap = np.linalg.norm(y_itrap - y_ref_itrap)
+
+    print(f"Ошибка явного метода Эйлера: {err_euler:.2e}")
+    print(f"Ошибка неявного метода Эйлера: {err_beuler:.2e}")
+    print(f"Ошибка неявного метода трапеций: {err_itrap:.2e}")
+
+    plt.figure(figsize=(12, 8))
+    plt.plot(t_ref, y_ref, 'k-', linewidth=5, label='Эталонное решение (solve_ivp)')
+    plt.plot(t_euler, y_euler, 'bo-', markersize=4, label='Явный метод Эйлера')
+    plt.plot(t_beuler, y_beuler, 'g^-', markersize=4, label='Неявный метод Эйлера')
+    plt.plot(t_itrap, y_itrap, 'm*-',linewidth=2.2, markersize=4, label='Неявный метод трапеций')
+
+    plt.xlabel('t')
+    plt.ylabel('y(t)')
+    plt.title('Задача 2.4: y\' = -2000y + t, y(0) = 0 (жесткая система)')
+    plt.legend()
+    plt.grid(True, alpha=0.3)
+    plt.show()
+
+    return t_euler, y_euler, t_beuler, y_beuler, t_itrap, y_itrap, t_ref, y_ref
+
+
+if __name__ == "__main__":
+    result2 = problem_2_4()
+```

archive/labs/lab1/scripts/1_4/1_4.m

@@ -0,0 +1,7 @@
+f = @(t,y) y.^2 - t;
+y0 = 1;
+tspan = [0 1.05]; % tspan = [0 1.5];
+opts = odeset('RelTol',1e-12,'AbsTol',1e-14);
+[t1,y1] = ode45(f,tspan,y0,opts);
+plot(t1,y1)
+print -dpng "plot.png"

archive/labs/lab1/scripts/2_4/2_4.m

@@ -0,0 +1,9 @@
+lambda = 2000;
+f = @(t,y) -lambda*y + t;
+y0 = 0;
+tspan = [0 1];
+opts = odeset('RelTol',1e-12,'AbsTol',1e-14);
+[t2,y2] = ode15s(f,tspan,y0,opts);
+plot(t2,y2)
+print -dpng "plot.png"
+

archive/labs/lab2/report.typ

@@ -0,0 +1,527 @@
+#set page(
+  paper: "a4",
+  // margin: (x: 1.8cm, y: 1.5cm),
+)
+#set text(
+  font: "New Computer Modern",
+  size: 14pt
+)
+#set par(
+  /*first-line-indent: (
+    amount: 1.5em,
+    all: true
+  ),*/
+  justify: true,
+  leading: 0.52em,
+)
+#show link: underline
+#set page(footer: context {
+  if counter(page).get().first() > 1 [
+    #align(center)[
+      #counter(page).display("1")
+    ]
+  ]
+  if counter(page).get().first() == 1 [
+    #align(center)[
+      Санкт-Петербург \ 2025
+    ]
+  ]
+})
+
+#set page(header: context {
+  if counter(page).get().first() == 1 [
+    #align(center)[
+      Санкт-Петербургский национальный исследовательский университет информационных технологий, механики и оптики
+    ]
+  ]
+})
+
+
+
+#show raw.where(block: false): box.with(
+  fill: luma(240),
+  inset: (x: 3pt, y: 0pt),
+  outset: (y: 3pt),
+  radius: 2pt,
+)
+
+#show raw.where(block: true): block.with(
+  fill: luma(240),
+  inset: 10pt,
+  radius: 4pt,
+)
+
+// title
+
+#for _ in range(5) { linebreak() }
+
+#align(center)[Расчетно графическая работа №2]
+
+#for _ in range(15) { linebreak() }
+
+
+#align(right)[Выполнили:]
+#align(right)[Левахин Лев]
+#align(right)[Останин Андрей]
+#align(right)[Дощенников Никита]
+#align(right)[Группы: К3221, К3240]
+#align(right)[Проверил:]
+#align(right)[Владимир Владимирович Беспалов]
+
+#pagebreak()
+
+#outline(title: "Содержание")
+
+#align(center)[=== Цель]
+
+Решение неоднородного линейного дифференциального уравнения разными способами. 
+
+Методы будут демонстрироваться на задаче Коши: 
+
+$
+y'' - 4y' + 5y = 4 e^(-2x), space y(0) = 0, space y'(0) = 1
+$
+
+#pagebreak()
+#align(center)[=== Метод специальной правой части]
+
+(Дощенников Никита)
+
+Запишем характеристическое уравнение: 
+
+$
+r^2 - 4r + 5 = 0
+$
+
+Найдем его корни. Посчитаем дискриминант: 
+
+$
+cal(D) = b^2 - 4 a c = 16 - 4 dot 5 = -4 lt 0
+$
+
+Тогда корни: 
+
+$
+r_1 = frac(4 + sqrt(-4), 2) = 2 + i, space.quad r_2 = frac(4 - sqrt(-4), 2) = 2 - i
+$
+
+Запишем общее решение для однородного уравнения. Так как корни комплексные, то общее решение соответствует форме: 
+
+$
+y_"общ" = e^(alpha x) dot (C_1 dot cos beta x + C_2 dot sin beta x)
+$
+
+В нашем случае $alpha = 2, beta = 1$. Тогда: 
+
+$
+y_"общ" = e^(2 x) dot (C_1 dot cos x + C_2 dot sin x)
+$
+
+Специальный вид правой части соответствует виду:
+
+$
+f(x) = P_n (x) dot e^(alpha x)
+$
+
+где $P_0 (x) = 4 space (n = 0), alpha = -2$. 
+
+Так как $alpha = -2$ не является корнем характеристического уравнения, частное решение представим в следующем виде: 
+
+$
+y_"частн" = e^(alpha x) dot Q_0 (x) = e^(-2 x) dot A
+$
+
+Найдем первую и вторую производные для частного решения: 
+
+$
+y'_"частн" eq (A e^(-2x))' = -2 A e^(-2x)
+$
+
+$
+y''_"частн" eq (-2 A e^(-2x))' = 4 A e^(-2x)
+$
+
+Подставим в исходное уравнение: 
+
+$
+4 A e^(-2x) + 8 A e^(-2x) + 5 A e^(-2x) = 4e^(-2x)
+$
+
+Вынеся $e^(-2x)$ за скобки получим: 
+
+$
+e^(-2x) (4A + 8A + 5A) = e^(-2x)(4)
+$
+
+Отсюда видно, что
+
+$
+4A + 8A + 5A = 4 arrow.double 17A = 4 arrow.double A = frac(4, 17)
+$
+
+Тогда 
+
+$
+y_"частн" = A e^(-2x) = frac(4, 17) e^(-2x)
+$
+
+Запишем общее решение: 
+
+$
+y = y_"общ" + y_"частн" = e^(2 x) dot (C_1 dot cos x + C_2 dot sin x) + frac(4, 17) e^(-2x)
+$
+
+Теперь решим задачу Коши. Чтобы определить коэффициенты, подставим $0$ вместо $x$:
+
+$
+y = C_1 + frac(4, 17)
+$
+
+по условию $y(0) = 0$, то есть 
+
+$
+C_1 + frac(4, 17) = 0 arrow.double C_1 = -frac(4, 17).
+$
+
+Найдем производную $y'$:
+
+$
+y' = C_1 (e^(2x) cos x)' + C_2 (e^(2x) sin x)' + frac(4, 17)(e^(-2x))' = \
+= C_1 (2e^(2x) cos x - e^(2x) sin x) + C_2 (2e^(2x) sin x + e^(2x) cos x) - frac(8, 17)e^(-2x) = \
+= e^(2x)((2C_1 + C_2) cos x + (2 C_2 - C_1) sin x) - frac(8, 17) e^(-2x)
+$
+
+По условию $y'(0) = 1$. Подставим $0$ вместо $x$ и получим: 
+
+$
+y'(0) = 2C_1 + C_2 - frac(8, 17) = 1
+$
+
+Подставим $C_1 = -frac(4, 17)$:
+
+$
+-frac(8, 17) + C_2 - frac(8, 17) = 1 arrow.double C_2 = frac(33, 17)
+$
+
+Подставим найденные коэффициенты в итоговое решение: 
+
+$
+y = e^(2x) (-frac(4, 17) cos x + frac(33, 17) sin x) + frac(4, 17) e^(-2x).
+$
+
+
+#pagebreak()
+#align(center)[=== Метод вариации постоянной]
+
+(Дощенников Никита)
+
+Напомним вид общего решения для однородного уравнения: 
+
+$
+y_"общ" = e^(2 x) dot (C_1 dot cos x + C_2 dot sin x)
+$
+
+Раскроем скобки и получим
+
+$
+y_"общ" = C_1 e^(2 x) cos x + C_2 e^(2 x) sin x
+$
+
+Сведем все к системе: 
+
+$
+cases(
+  C'_1 (x) y_1 + C'_2 (x) y_2 = 0,
+  C'_1 (x) y'_1 + C'_2 (x) y'_2 = f(x)
+)
+$
+
+где $y_1 = e^(2x) cos x, space y_2 = e^(2x) sin x, space f(x) = 4 e^(-2x)$.
+
+Найдем производные для $y_1$ и $y_2$:
+
+$
+y_1' = (e^(2x) cos x)' = 2 e^(2x) cos x - e^(2x) sin x
+$
+
+$
+y_2' = (e^(2x) sin x)' = 2 e^(2x) sin x + e^(2x) cos x 
+$
+
+Тогда система примет следующий вид: 
+
+$
+cases(
+  C'_1(x) (e^(2x) cos x) + C'_2(x) (e^(2x) sin x) = 0,
+  C'_1(x) (2 e^(2x) cos x - e^(2x) sin x) + C'_2(x) (2 e^(2x) sin x + e^(2x) cos x) = 4e^(-2x),
+)
+$
+
+Из первого уравнения системы: 
+
+$
+C'_1 = -C'_2 frac(e^(2x) sin x, e^(2x) cos x) = -C'_2 tg x
+$
+
+Подставим во второе уравнение системы и вынесем $C'_2 e^(2x)$: 
+
+$
+C'_2 e^(2x) (-tg x(2 cos x - sin x) + 2 sin x + cos x) = 4 e^(-2x)
+$
+
+Упростим
+
+$
+-tg x(2 cos x - sin x) = -frac(sin x, cos x) (2 cos x - sin x) = -2 sin x + frac(sin^2 x, cos x)
+$
+
+Тогда
+
+$
+(-2 sin x + frac(sin^2 x, cos x)) + 2 sin x + cos x = frac(sin^2 x + cos^2 x, cos x) = frac(1, cos x)
+$
+
+Вернемся в уравнение 2 системы. С учетом упрощения оно примет следующий вид:
+
+$
+C'_2 e^(2x) dot frac(1, cos x) = 4 e^(-2x) 
+$
+
+Отсюда 
+
+$
+C_2' = 4e^(-4x) cos x
+$
+
+Найдем $C'_1$
+
+$
+C'_1 = -C'_2 tg x = -4e^(-4x) cos x frac(sin x, cos x) = -4e^(-4x) sin x
+$
+
+Проинтегрируем и получим
+
+$
+C_1 = -4 integral e^(-4x) sin x d x = frac(4, 17) e^(-4x) (4 sin x + cos x) \
+C_2 = 4 integral e^(-4x) cos x d x = frac(4, 17) e^(-4x) (4 cos x - sin x)
+$
+
+Подставим в формулу частного решения: 
+
+$
+y_"частн" = C_1 (x) y_1 + C_2 (x) y_2 = \
+= frac(4, 17) e^(-4x) (4 sin x + cos x) dot e^(2x) cos x + frac(4, 17) e^(-4x) (4 cos x - sin x) dot e^(2x) sin x = \
+= frac(4, 17) e^(-2x) ((4 sin x + cos x ) cos x + (4 cos x - sin x) sin x) = \
+= frac(4, 17) e^(-2x) (4 sin x cos x + 4 cos x sin x + cos^2 x - sin^2 x) = \
+= frac(4, 17) e^(-2x) (8 sin x cos x + (cos^2 x + sin^2 x)) = \
+= frac(4, 17) e^(-2x) (8 sin x cos x + 1) = \
+= frac(4, 17) e^(-2x) (cos 2 x + 4 sin 2x)
+$
+
+Так как $e^(-2x) cos 2 x$ и $e^(-2x) sin 2x$ -- решения однородного уравнения, то выражение $frac(4, 17) e^(-2x) (cos 2 x + 4 sin 2 x)$ -- это частное решение и часть однородного решения. Тогда окончательно получим
+
+$
+y_"частн" = frac(4, 17) e^(-2x)
+$
+
+Запишем общее решение: 
+
+$
+y = y_"общ" + y_"частн" = e^(2x) (-frac(4, 17) cos x + frac(33, 17) sin x) + frac(4, 17) e^(-2x).
+$
+
+Дальнейшее решение задачи Коши идентично описанному в первом методе.
+
+
+#pagebreak()
+#align(center)[=== Операционный метод] (Левахин Лев)
+
+
+Применим преобразование Лапласа к исходному уравнению:
+$ L[y''] - 4L[y'] + 5L[y] = 4L[e^(-2x)] $
+
+Используем свойства преобразования Лапласа:
+$ (p^2 Y(p) - p y(0) - y'(0)) - 4(p Y(p) - y(0)) + 5Y(p) = 4/(p+2) $
+
+Подставим начальные условия $y(0)=0$, $y'(0)=1$:
+$ (p^2 Y(p) - 1) - 4p Y(p) + 5Y(p) = 4/(p+2) $
+
+Сгруппируем члены с $Y(p)$:
+$ (p^2 - 4p + 5)Y(p) - 1 = 4/(p+2) $
+$ (p^2 - 4p + 5)Y(p) = 1 + 4/(p+2) = (p+6)/(p+2) $
+
+Выразим $Y(p)$:
+$ Y(p) = (p+6)/((p+2)(p^2 - 4p + 5)) $
+
+Разложим на простейшие дроби. Для этого сначала представим:
+$ Y(p) = A/(p+2) + (B*p + C)/(p^2 - 4p + 5) $
+
+Умножим обе части на $(p+2)(p^2 - 4p + 5)$:
+$ p+6 = A(p^2 - 4p + 5) + (B*p + C)(p+2) $
+$ p+6 = A p^2 - 4A p + 5A + B p^2 + 2B p + C p + 2C $
+$ p+6 = (A+B)p^2 + (-4A + 2B + C)p + (5A + 2C) $
+
+Получим систему уравнений:
+$
+cases(
+  A + B = 0,
+  -4A + 2B + C = 1,
+  5A + 2C = 6
+)
+$
+A + B = 0, при p^2
+-4A + 2B + C = 1, при p^1
+5A + 2C = 6, при p^0
+
+Из первого уравнения: $B = -A$. Подставим в остальные:
+$
+cases(
+  -4A - 2A + C = 1,
+  5A + 2C = 6
+)
+$
+$
+cases(
+  -6A + C = 1,
+  5A + 2C = 6
+)$
+
+Решим систему: из первого $C = 1 + 6A$, подставим во второе:
+$ 5A + 2(1 + 6A) = 6 $
+$ 5A + 2 + 12A = 6 $
+$ 17A = 4 $
+$ A = 4/17 $
+
+Тогда:
+$ B = -4/17 $
+$ C = 1 + 6*(4/17) = 1 + 24/17 = 41/17 $
+
+Итак:
+$ Y(p) = (4/17)/(p+2) + ((-4/17)p + 41/17)/(p^2 - 4p + 5) $
+
+Преобразуем второе слагаемое. Заметим, что $p^2 - 4p + 5 = (p-2)^2 + 1$:
+$ Y(p) = (4/17)/(p+2) + (-4/17 * p + 41/17)/((p-2)^2 + 1) $
+
+Выделим в числителе второй дроби слагаемое $p-2$:
+$ -4/17 * p + 41/17 = -4/17 * (p-2) - 8/17 + 41/17 = -4/17 * (p-2) + 33/17 $
+
+Тогда:
+$ Y(p) = (4/17)/(p+2) + (-4/17 * (p-2) + 33/17)/((p-2)^2 + 1) $
+$ Y(p) = (4/17)/(p+2) - (4/17) * (p-2)/((p-2)^2 + 1) + (33/17)/((p-2)^2 + 1) $
+
+Применим обратное преобразование Лапласа, используя свойства:
+$ L^(-1)[1/(p+2)] = e^(-2t) $
+$ L^(-1)[(p-2)/((p-2)^2 + 1)] = e^(2t) cos t $
+$ L^(-1)[1/((p-2)^2 + 1)] = e^(2t) sin t $
+
+Получаем решение:
+$ y(x) = (4/17)e^(-2x) - (4/17)e^(2x) cos x + (33/17)e^(2x) sin x $
+
+Или в более компактной форме:
+$ y(x) = e^(2x)(-frac(4, 17) cos x + frac(33, 17) sin x) + frac(4, 17)e^(-2x) $
+
+Полученное решение полностью совпадает с результатами предыдущих методов.
+
+
+#pagebreak()
+#align(center)[=== Метод разложения в ряд] (Останин Андрей)
+
+Ищем решение в виде степенного ряда:
+$
+y(x) = sum_(n=0)^infinity a_(n) x^n
+$
+
+Тогда производные имеют вид:
+$
+y'(x) = sum_(n=1)^infinity n a_(n) x^(n-1), quad
+y''(x) = sum_(n=2)^infinity n(n-1) a_(n) x^(n-2)
+$
+
+Правая часть уравнения раскладывается в ряд Тейлора:
+$
+4 e^(-2x) = sum_(n=0)^infinity b_(n) x^n, quad
+b_(n) = frac(4 (-2)^n, n!)
+$
+
+Приведём ряды для производных к одинаковым степеням $x^n$:
+$
+y'' = sum_(n=0)^infinity (n+2)(n+1) a_(n+2) x^n
+$
+
+$
+y' = sum_(n=0)^infinity (n+1) a_(n+1) x^n
+$
+
+Подставляя полученные выражения в исходное дифференциальное уравнение, получаем:
+$
+sum_(n=0)^infinity [
+(n+2)(n+1) a_(n+2)
+- 4 (n+1) a_(n+1)
++ 5 a_(n)
+] x^n
+=
+sum_(n=0)^infinity b_(n) x^n
+$
+
+Отсюда следует рекуррентное соотношение для коэффициентов:
+$
+(n+2)(n+1) a_(n+2)
+- 4 (n+1) a_(n+1)
++ 5 a_(n)
+=
+b_(n)
+$
+
+Начальные условия задачи Коши дают:
+$
+a_(0) = y(0) = 0, quad a_(1) = y'(0) = 1
+$
+
+Найдём первые коэффициенты ряда.
+
+При $n = 0$:
+$
+2 a_(2) - 4 a_(1) + 5 a_(0) = 4
+$
+$
+2 a_(2) - 4 = 4 arrow.double a_(2) = 4
+$
+
+При $n = 1$:
+$
+6 a_(3) - 8 a_(2) + 5 a_(1) = -8
+$
+$
+6 a_(3) - 32 + 5 = -8 arrow.double a_(3) = frac(19, 6)
+$
+
+При $n = 2$:
+$
+12 a_(4) - 12 a_(3) + 5 a_(2) = 8
+$
+$
+12 a_(4) - 38 + 20 = 8 arrow.double a_(4) = frac(13, 6)
+$
+
+При $n = 3$:
+$
+20 a_(5) - 16 a_(4) + 5 a_(3) = - frac(16, 3)
+$
+$
+20 a_(5) - frac(208, 6) + frac(95, 6) = - frac(16, 3)
+arrow.double a_(5) = frac(27, 40)
+$
+
+Таким образом, решение в виде степенного ряда имеет вид:
+$
+y(x) =
+x
++ 4 x^2
++ frac(19, 6) x^3
++ frac(13, 6) x^4
++ frac(27, 40) x^5
++ dots
+$
+
+Полученное разложение совпадает с рядом Тейлора точного решения, найденного ранее другими методами.
+

archive/notes/notes.typ

@@ -0,0 +1,26 @@
+#set text(size: 1.3em)
+
+#align(center)[=== Однородные]
+
+1. Если $k_1 eq.not k_2$ - действительные и различные, то #align(center)[$y eq C_1 dot e^(k_1 x) + C_2 dot e^(k_2 x)$]
+
+2. Если $k_1 eq k_2$ - действительные и совпавшие, то #align(center)[$y eq e^(k_1 x) dot (C_1 plus C_2 x)$]
+
+3. Если $k_(1, 2) eq alpha plus.minus beta i$ - комплексные корни, то #align(center)[$y eq e^(alpha x) dot (C_1 dot cos beta x plus C_2 dot sin beta x).$]
+
+#align(center)[=== Неоднородные]
+
+1. Специальный вид правой части: #align(center)[$f(x) eq P_n (x) dot e^(alpha x).$]
+    - Если $alpha$ не является корнем характеристического уравнения, то #align(center)[$y_"частное" eq e^(alpha x) dot Q_n (x).$]
+    - Если $alpha$ - корень характеристического уравнения кратности $S$ ($S in {1, 2}$), то #align(center)[$y_"частное" eq x^S dot e^(alpha x) dot Q_n (x).$]
+
+2. Специальный вид правой части: #align(center)[$f(x) eq e^(alpha x) dot (P_n (x) dot cos beta x + Q_m (x) dot sin beta x), space.quad N eq max(n, m).$]
+
+    - Если $alpha plus.minus beta i$ не являются корнями характеристического уравнения, то #align(center)[$y_"частное" eq e^(alpha x) dot (P_N (x) dot cos beta x + Q_N (x) dot sin beta x).$]
+    - Если $alpha plus.minus beta i$ - корни характеристического уравнения, то #align(center)[$y_"частное" eq x dot e^(alpha x) dot (P_N (x) dot cos beta x + Q_N (x) dot sin beta x).$]
+
+=== Метод Лагранжа (Метод вариации постоянных)
+
+$
+cases(C'_1 (x) y_1 plus C'_2 (x) y_2 eq 0, C'_1 (x) y'_1 plus C'_2 (x) y'_2 eq f(x))
+$

archive/practice/lab1/main.py

@@ -0,0 +1,82 @@
+import sys
+import matplotlib.pyplot as plt
+
+def f(x: float, y: float) -> float:
+    # f(x, y) = x^2 - 2y
+    return x**2 - 2 * y
+
+def euler(h: float, x_min: float, x_max: float, y_0: float):
+    # y_(n + 1) = y_n + h * f(x_n, y_n)
+    x_i: list[float] = [x_min]
+    y_i: list[float] = [y_0]
+
+    f_i = f(x_i[-1], y_i[-1])
+    h_f_i = h * f_i
+
+    for i in range(int((x_max - x_min) / h)):
+        y_i.append(y_i[-1] + h_f_i)
+        x_i.append(x_i[-1] + h)
+
+        f_i = f(x_i[-1], y_i[-1])
+        h_f_i = h * f_i
+    
+    return x_i, y_i
+
+def trapezoid(h: float, x_min: float, x_max: float, y_0: float):
+    # y_(n + 1) = y_n + h * (f(x_n, y_n) + f(x_(n + 1), y~_(n + 1)))/2
+    # y~_(n + 1) = y_n + h * f(x_n, y_n)
+    
+    x_i: list[float] = [x_min]
+    y_i: list[float] = [y_0]
+
+    y_tilde = y_i[-1] + h * f(x_i[-1], y_i[-1])
+
+    for i in range(int((x_max - x_min) / h)):
+        y = y_i[-1] + h * (f(x_i[-1], y_i[-1]) + f(x_i[-1] + h, y_tilde)) / 2
+        y_i.append(y)
+        x_i.append(x_i[-1] + h)
+
+        y_tilde = y_i[-1] + h * f(x_i[-1], y_i[-1])
+    
+    return x_i, y_i
+
+def rk4(h: float, x_min: float, x_max: float, y_0: float):
+    # k_1 = f(x_n, y_n)
+    # k_2 = f(x_n + h/2, y_n + h/2 k_1)
+    # k_3 = f(x_n + h/2, y_n + h/2 k_2)
+    # k_4 = f(x_n + h, y_n + hk_3)
+    # y_(n + 1) = y_n + h/6 (k_1 + 2k_2 + 2k_3 + k_4)
+
+    x_i: list[float] = [x_min]
+    y_i: list[float] = [y_0]
+
+    for i in range(int((x_max - x_min) / h)):
+        k_1 = f(x_i[-1], y_i[-1])
+        k_2 = f(x_i[-1] + h/2, y_i[-1] + h/2 * k_1)
+        k_3 = f(x_i[-1] + h/2, y_i[-1] + h/2 * k_2)
+        k_4 = f(x_i[-1] + h, y_i[-1] + h * k_3)
+        y = y_i[-1] + h/6 * (k_1 + 2 * k_2 + 2 * k_3 + k_4)
+
+        y_i.append(y)
+        x_i.append(x_i[-1] + h)
+
+    return x_i, y_i
+
+def draw_graph(x: list, y: list):
+    plt.plot(x, y)
+    plt.show()
+
+def main() -> None:
+    # y' = t^2 - 2y, y(0) = 1, t in [0, 1]
+
+    methods = {"euler": euler,
+               "trapezoid": trapezoid,
+               "rk4": rk4}
+
+    if sys.argv[1] in methods.keys():
+        x, y = methods[sys.argv[1]](0.01, 0, 1, 1)
+
+    draw_graph(x, y)
+
+if __name__ == "__main__":
+    main()