Улучшение расчетной модели. Интегрирование с адаптивным шагом.

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

Сейчас численное интегрирование идет с постоянным временным шагом. Не важно при этом, являются взаимодействия интенсивными (к примеру, пролет на "гравипраще" спутника рядом с массивным объектом (Юпитер, например), так что за считанные минуты скорость спутника увеличивается на километры в секунду) или медленным (вековой дрейф обледеневшего булыжника в облаке Оорта).

Желая получить высокую точность, мы можем уменьшать шаг интегрирования до минут и даже секунд, но будем платить за это падением производительности. И наоборот - можно ускорять расчеты, жертвуя точностью.

К счастью, существуют способы оценки точности численного интегрирования, благодаря которым можно делать вывод о правильности выбора шага. Один из простых и явных методов оценки - правило Рунге (стр.10, формула 68) или (стр.15, формула 1.12.19). Внимание, внутри много теоретического матана, стоящего "под капотом" формул численных методов.

Принцип достаточно прост - проводится расчет по выбранной величине шага, затем выполняется такой же расчет в две ступени с вдвое меньшим шагом. Результаты сравниваются с учетом порядка метода интегрирования, и если ошибка по абсолютной величине больше заранее заданной точности, то нам нужно уменьшить шаг в два раза.

И наоборот, если при сравнении двух обычных шагов и одного - с удвоенным шагом ошибка остается в допустимых пределах, то мы можем удвоить величину шага.

С 54 по 59 строку - проверка на необходимость уменьшить шаг, с 60 по 65 - на увеличение шага

Для отработки нового приема выберем функцию с как можно более неоднородным поведением (чтобы и производная вела себя так же). К примеру:

Выглядит малость запутанно. Это наша функция

Выглядит откровенно гадко. Это производная, с ней и будем работать. А предыдущая функция для нее - первообразная

Запустим интегрирование с точки X0 = -3,5 и до X1 = -0,75.

Начальный шаг - 0,01.

График нашей функции выведем черным, он будет эталоном. Точки численного решения выведем поверх красным. Сравниваем и рассуждаем.

Итоги численного расчета достаточно плотно ложатся на график. В экстремумах, где функция круто возрастает или падает, точки лежат плотнее

Проверим расчет по формуле Ньютона-Лейбница, а результат выведем в консоль:

К формуле Ньютона-Лейбница мы обращаемся при выводе в консоль на строке 78.

Результат вывод в консоль. Абсолютная погрешность чуть меньше двух тысячных. Относительная - 1.2%. Неплохо

Метод управления шагом опробован. Пока что в тепличных условиях для функции одной переменной (пусть и достаточно страшной функции). Теперь надо попытаться прикрутить его к нашей модели орбитального движения. Спойлер. Это будет долго и тяжко.

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