Реализация адаптивной числовой интеграции
Я написал пакет Python 3, который может интегрировать системы обыкновенных дифференциальных уравнений. При написании этого я столкнулся с проблемой, о которой я не могу найти никакой информации.
Проблема:
По сути, у меня есть адаптивный алгоритм, который регулирует размер шага в соответствии с локальной ошибкой (фактический метод - это реализация алгоритма RKF с коэффициентами Кэша-Карпа) и после определенного момента, если размер шага достаточно велик, метод будет превышать время окончания.
Мое решение:
Я установил его так, что, если текущий размер шага приведет к тому, что интеграция превысит время окончания, тогда я установлю размер шага равным половине значения оставшегося времени. В форме уравнения: dt = (t_end - t_current)/2
Но это, при определенных требованиях относительной ошибки, используемой в адаптивном алгоритме, приведет к тому, что шаг по времени станет крошечным. Приблизительно порядка 1e-15, что приводит к численной нестабильности, поскольку t_current никогда не увеличивается из-за ограничений точности станка.
Мой вопрос:
Мне было интересно, существует ли более эффективный способ реализации алгоритма интегрирования, чтобы время окончания было точно ~15 десятичных разрядов относительно ожидаемого времени окончания без создания таких числовых нестабильностей.