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

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

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

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

Шаг 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⁻⁸ на одинаковой сетке выдачи.

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

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