Численные методы без магии

Метод Эйлера, адаптивный решатель и проверка точности результата.

От уравнения к последовательности точек

Компьютер не хранит непрерывную функцию целиком. Он вычисляет состояния в конечном наборе точек времени. Простейший метод использует текущую производную для прогноза следующего состояния. Это явный метод Эйлера:

Шаг h измеряется в секундах. Метод предполагает, что на маленьком интервале наклон почти не меняется. Для электромеханической системы обе новые переменные нужно вычислять из одного старого состояния. Если обновить ток и сразу использовать новый ток для скорости, получится уже другой численный метод.

i_old, w_old = i, w
di = (U - R*i_old - k*w_old)/L
dw = (k*i_old - b*w_old - load)/J
i = i_old + h*di
w = w_old + h*dw

Устойчивость не равна точности

Для тестового уравнения x′ = −x/τ метод Эйлера даёт xn+1 = (1 − h/τ)xn. Модуль множителя должен быть меньше единицы, иначе ошибка будет расти. Для этого примера условие устойчивости: 0 < h < 2τ. Но шаг 1,5τ, хотя формально устойчив, даёт чередование знаков, которого нет в точном решении.

Для хорошей точности шаг обычно выбирают существенно меньше самой быстрой характерной постоянной времени, а затем проверяют сходимость. Универсального числа точек для всех задач нет. Чем резче воздействие и чем быстрее динамика, тем внимательнее нужно исследовать сетку.

Адаптивный решатель

В SciPy функция solve_ivp решает задачу y′ = f(t,y) с заданным начальным состоянием. Метод RK45 автоматически меняет внутренний шаг по оценке локальной ошибки. Параметр t_eval задаёт точки выдачи результата, а не фиксированный шаг интегрирования. Тысяча красивых точек на графике сама по себе не гарантирует точного расчёта.

sol = solve_ivp(motor, (0, 1.5), [0, 0],
                t_eval=np.linspace(0, 1.5, 1501),
                rtol=1e-8, atol=[1e-10, 1e-8],
                max_step=0.005)
assert sol.success, sol.message
current, speed = sol.y

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

Скачки и переключения

Если нагрузка скачком меняется в момент 0,5 с, разумно закончить первый расчёт точно в этой точке, взять его последнее состояние и начать второй интервал. Состояние сохраняется, но правая часть меняется. Это особенно полезно, если короткий импульс мог бы оказаться между внутренними шагами решателя.

Для событий, время которых заранее неизвестно, например пересечения скоростью заданного порога, существуют функции событий. Важно учитывать, что поиск события тоже связан с шагами и может пропустить несколько пересечений внутри одного шага.

Три независимые проверки

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

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

Когда задача становится жёсткой

Если электрическая динамика намного быстрее механической, явный метод может потребовать очень много шагов, даже когда интересует медленный разгон. Это один из признаков жёсткости. Для таких задач полезны неявные методы Radau или BDF. Однако менять метод следует после проверки масштаба параметров и единиц, а не использовать его для маскировки ошибки в уравнениях.

Практическая работа

Решите RL-задачу методом Эйлера с шагами 2; 1; 0,5 и 0,25 мс. Для каждого найдите максимальную абсолютную ошибку относительно аналитического решения. Постройте зависимость ошибки от шага. Затем сравните решение двигателя при rtol = 10⁻⁵ и 10⁻⁸ на одинаковой сетке выдачи.

Область устойчивости явного Эйлера

Для тестового уравнения x′ = λx шаг Эйлера даёт x_n+1 = (1 + hλ)x_n. Чтобы затухающее точное решение не превращалось в растущую последовательность, нужно |1 + hλ| < 1. Для отрицательного вещественного λ это означает 0 < h < 2/|λ|.

Условие проверяют для всех собственных значений матрицы A. Для базового двигателя быстрый полюс −97,407 с⁻¹ ограничивает h значением около 0,02053 с. Но шаг 0,02 с, хотя формально устойчив, воспроизводит быстрый режим плохо: множитель ≈ −0,948 создаёт чередование знаков. Для монотонного затухания этой моды нужно h ≤ 1/97,407 ≈ 0,01027 с, а для хорошей точности обычно ещё меньше.

Методы Рунге — Кутты и порядок ошибки

Явный Эйлер использует наклон в начале шага и имеет глобальную ошибку порядка h на фиксированном интервале для достаточно гладкой задачи. Классический RK4 использует четыре оценки наклона:

Глобальная ошибка RK4 имеет порядок h⁴ при достаточной гладкости. Скачок входа внутри шага нарушает предпосылки: разделяйте расчёт в известные моменты коммутации. Больший порядок не устраняет ограничение устойчивости явного метода на жёстких задачах.

Как измерить сходимость

Выполните расчёты с шагами h, h/2 и h/4 и сравните значения в одинаковые моменты времени. Для одной величины задайте d₁ = |x_h − x_h/2| и d₂ = |x_h/2 − x_h/4|. Если ведущая ошибка пропорциональна h^p, то:

Оценка работает в асимптотической области, когда погрешность округления и разрывы не доминируют. Для вектора состояния нельзя напрямую складывать амперы и рад/с: нормируйте компоненты на выбранные масштабы или анализируйте отдельно.

В адаптивном решателе rtol и atol задают масштабы контроля локальной ошибки: atol + rtol·|x|. Они не являются гарантией такой же глобальной ошибки. Параметр t_eval задаёт точки выдачи результата, а не внутренний шаг интегрирования. Для жёсткой системы сравните RK45 с Radau или BDF, проверяя точность и число вычислений.

Дополнительное чтение: SciPy: параметры solve_ivp и выбор решателя. Обозначения в источнике могут отличаться от принятых в уроке.

Задачи для самостоятельного решения

1. Решите x′ = −10x, x₀ = 1 до t = 0,1 с методом Эйлера с h = 0,05; 0,025; 0,0125 с. Сравните с e⁻¹.

2. Для λ = −2 + 20j с⁻¹ найдите верхнюю границу шага Эйлера. Почему оценка только по действительной части неверна?

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

Что запомнить

Сначала проверяют постановку, затем устойчивость и точность численного метода. Не путайте сетку графика с внутренними шагами решателя. Основной справочник по параметрам — официальная документация solve_ivp.