Численные методы: интегрирование и решение ОДУ

Формулы прямоугольников, трапеций и Симпсона, оценка погрешности и правило Рунге, методы Эйлера и Рунге-Кутта для дифференциальных уравнений, устойчивость схем.

Численные методы нужны там, где аналитического решения нет: интеграл не берётся, дифференциальное уравнение нелинейно. Практическая часть курсовой обычно состоит в реализации нескольких схем и сравнении их точности.

Численное интегрирование

Отрезок [a, b] делят на n частей шагом h = (b − a)/n

Левых прямоугольников:  I ≈ h·Σ f(xᵢ),        i = 0…n−1
Правых:                 I ≈ h·Σ f(xᵢ),        i = 1…n
Средних:                I ≈ h·Σ f(xᵢ + h/2)

Трапеций:
  I ≈ h·[ f(a)/2 + f(x₁) + … + f(x_{n−1}) + f(b)/2 ]

Симпсона (n обязательно чётное):
  I ≈ (h/3)·[ f(a) + 4·(нечётные) + 2·(чётные) + f(b) ]

Правило Рунге для оценки погрешности:

  R ≈ (I_h/2 − I_h) / (2^p − 1)

p — порядок метода (2 для трапеций, 4 для Симпсона)

Уточнённое значение (экстраполяция Ричардсона):
  I ≈ I_h/2 + R

Алгоритм с автоматическим выбором шага:
считаем с шагом h и h/2, если |R| > ε — дробим шаг дальше.

def trapezoid(f, a, b, n):
    h = (b - a) / n
    s = (f(a) + f(b)) / 2
    for i in range(1, n):
        s += f(a + i * h)
    return s * h

def simpson(f, a, b, n):
    if n % 2:
        n += 1                       # число отрезков должно быть чётным
    h = (b - a) / n
    s = f(a) + f(b)
    for i in range(1, n):
        s += f(a + i * h) * (4 if i % 2 else 2)
    return s * h / 3

def integrate_auto(f, a, b, eps=1e-8, p=4):
    """Автоматический выбор шага по правилу Рунге."""
    n = 4
    prev = simpson(f, a, b, n)
    while n < 10**7:
        n *= 2
        cur = simpson(f, a, b, n)
        r = (cur - prev) / (2**p - 1)
        if abs(r) < eps:
            return cur + r, n
        prev = cur
    raise RuntimeError('не сошлось')

import math
val, n = integrate_auto(lambda x: math.exp(-x*x), 0, 1)
print(f'∫e^(-x²)dx на [0,1] = {val:.10f} при n = {n}')
# 0.7468241328, точное значение 0.746824132812...

Задача Коши

Дано:  y' = f(x, y),  y(x₀) = y₀
Найти: y(x) на отрезке

Метод Эйлера (первого порядка):
  y_{i+1} = yᵢ + h·f(xᵢ, yᵢ)

Модифицированный Эйлера, предиктор-корректор (второго порядка):
  ỹ = yᵢ + h·f(xᵢ, yᵢ)
  y_{i+1} = yᵢ + (h/2)·[ f(xᵢ, yᵢ) + f(x_{i+1}, ỹ) ]

Рунге-Кутта 4-го порядка:
  k₁ = h·f(xᵢ, yᵢ)
  k₂ = h·f(xᵢ + h/2, yᵢ + k₁/2)
  k₃ = h·f(xᵢ + h/2, yᵢ + k₂/2)
  k₄ = h·f(xᵢ + h, yᵢ + k₃)
  y_{i+1} = yᵢ + (k₁ + 2k₂ + 2k₃ + k₄)/6

def euler(f, x0, y0, x_end, h):
    xs, ys = [x0], [y0]
    x, y = x0, y0
    while x < x_end - 1e-12:
        y += h * f(x, y)
        x += h
        xs.append(x); ys.append(y)
    return xs, ys

def runge_kutta4(f, x0, y0, x_end, h):
    xs, ys = [x0], [y0]
    x, y = x0, y0
    while x < x_end - 1e-12:
        k1 = h * f(x, y)
        k2 = h * f(x + h/2, y + k1/2)
        k3 = h * f(x + h/2, y + k2/2)
        k4 = h * f(x + h,   y + k3)
        y += (k1 + 2*k2 + 2*k3 + k4) / 6
        x += h
        xs.append(x); ys.append(y)
    return xs, ys

# Проверка на y' = y, y(0) = 1 → точное решение e^x
import math
for h in (0.1, 0.05):
    _, ye = euler(lambda x, y: y, 0, 1, 1, h)
    _, yr = runge_kutta4(lambda x, y: y, 0, 1, 1, h)
    print(f'h={h}: Эйлер {ye[-1]:.6f} (ошибка {abs(ye[-1]-math.e):.2e}), '
          f'РК4 {yr[-1]:.8f} (ошибка {abs(yr[-1]-math.e):.2e})')

Системы и уравнения высших порядков

Уравнение n-го порядка сводят к системе n уравнений первого:

  y'' + 2y' + 5y = 0

  Замена: u = y, v = y'

  u' = v
  v' = −2v − 5u

Далее применяют тот же Рунге-Кутта, но к вектору (u, v).
В коде это означает, что y становится списком, а f возвращает список.

Устойчивость

Явные методы устойчивы только при достаточно малом шаге. Для жёстких систем — где процессы идут с сильно разными скоростями — ограничение на шаг становится непрактичным, и применяют неявные схемы: неявный Эйлер, методы Гира. Они требуют решения уравнения на каждом шаге, зато устойчивы при любом шаге.

Признак жёсткости: собственные числа матрицы системы
отличаются на порядки.

Пример: y' = −1000·y + 3000 − 2000·e^(−x)

Явный Эйлера требует h < 0,002 — иначе решение расходится,
хотя само решение гладкое и медленно меняется.
Неявный Эйлера считает это уравнение при h = 0,1 без проблем.

Что показать в отчёте

1. Реализовать не менее двух методов разного порядка.
2. Сравнить с точным решением там, где оно известно, и построить график ошибки.
3. Проверить порядок точности численно: при уменьшении шага вдвое ошибка должна упасть в 2^p раз.
4. Применить правило Рунге для автоматического выбора шага.
5. Построить таблицу: шаг, число вычислений функции, погрешность, время работы.
6. Сделать вывод о применимости методов к вашей задаче.