Адаптивный размер шага

В математике и численном анализе адаптивный размер шага используется в некоторых методах численного решения обыкновенных дифференциальных уравнений (включая частный случай численного интегрирования ) с целью контроля ошибок метода и обеспечения свойств устойчивости , таких как A- стабильность . Использование адаптивного размера шага имеет особое значение, когда размер производной сильно варьируется. Например, при моделировании движения спутника вокруг Земли по стандартной орбите Кеплера метода фиксированного шага по времени, такого как метод Эйлера может быть достаточно . Однако все сложнее, если кто-то хочет смоделировать движение космического корабля, принимая во внимание как Землю, так и Луну, как в задаче трех тел . Там возникают сценарии, в которых можно делать большие шаги по времени, когда космический корабль находится далеко от Земли и Луны, но если космический корабль приближается к столкновению с одним из планетарных тел, тогда необходимы малые шаги по времени. метод Ромберга и Рунге-Кутта-Фельберг являются примерами методов численного интегрирования, в которых используется адаптивный размер шага.

Пример

Для простоты в следующем примере используется простейший метод интегрирования — метод Эйлера ; на практике Рунге-Кутты, предпочтительны методы более высокого порядка, такие как методы из-за их превосходных свойств сходимости и стабильности.

Рассмотрим задачу начального значения

$y'(t)=f(t,y(t)),\qquad y(a)=y_{a}$

где y и f могут обозначать векторы (в этом случае это уравнение представляет собой систему связанных ОДУ с несколькими переменными).

Нам даны функция f ( t , y ) и начальные условия ( a , ya ₎ , и мы заинтересованы в поиске решения при t = b . Пусть y ( b ) обозначает точное решение в точке b , и пусть y _b обозначает решение, которое мы вычисляем. Мы пишем $y_{b}+\varepsilon =y(b)$, где $\varepsilon$ – ошибка численного решения.

Для последовательности ( t _n ) значений t с t _n = a + nh метод Эйлера дает приближения к соответствующим значениям y ( t _n ) как

$y_{n+1}^{(0)}=y_{n}+hf(t_{n},y_{n})$

Локальная ошибка усечения этого приближения определяется выражением

$\tau _{n+1}^{(0)}=y(t_{n+1})-y_{n+1}^{(0)}$

и по теореме Тейлора можно показать, что (при условии, что f достаточно гладкое) локальная ошибка усечения пропорциональна квадрату размера шага:

$\tau _{n+1}^{(0)}=ch^{2}$

где c — некоторая константа пропорциональности.

Мы отметили это решение и его ошибку значком $(0)$.

Значение c нам неизвестно. Давайте теперь снова применим метод Эйлера с другим размером шага, чтобы сгенерировать второе приближение к y ( t _{n +1} ). Мы получаем второе решение, которое обозначим значком $(1)$. Возьмите новый размер шага равным половине исходного размера шага и примените два шага метода Эйлера. Это второе решение, по-видимому, более точное. Поскольку нам приходится применять метод Эйлера дважды, локальная ошибка (в худшем случае) вдвое превышает исходную ошибку.

$y_{n+{\frac {1}{2}}}=y_{n}+{\frac {h}{2}}f(t_{n},y_{n})$

$y_{n+1}^{(1)}=y_{n+{\frac {1}{2}}}+{\frac {h}{2}}f(t_{n+{\frac {1}{2}}},y_{n+{\frac {1}{2}}})$

$\tau _{n+1}^{(1)}=c\left({\frac {h}{2}}\right)^{2}+c\left({\frac {h}{2}}\right)^{2}=2c\left({\frac {h}{2}}\right)^{2}={\frac {1}{2}}ch^{2}={\frac {1}{2}}\tau _{n+1}^{(0)}$

$y_{n+1}^{(1)}+\tau _{n+1}^{(1)}=y(t+h)$

Здесь мы предполагаем коэффициент ошибки $c$ постоянен на интервале $[t,t+h]$. В действительности скорость его изменения пропорциональна $y^{(3)}(t)$. Вычитание решений дает оценку ошибки:

$y_{n+1}^{(1)}-y_{n+1}^{(0)}=\tau _{n+1}^{(1)}$

Эта локальная оценка ошибки имеет точность третьего порядка.

Оценку локальной ошибки можно использовать для определения размера шага. $h$ должны быть изменены для достижения желаемой точности. Например, если местная толерантность ${\text{tol}}$ разрешено, мы могли бы позволить h развиваться следующим образом:

$h\rightarrow 0.9\times h\times \min \left(\max \left(\left({\frac {\text{tol}}{2\left|\tau _{n+1}^{(1)}\right|}}\right)^{1/2},0.3\right),2\right)$

The $0.9$ является фактором безопасности, гарантирующим успех при следующей попытке. Минимум и максимум предназначены для предотвращения резких изменений по сравнению с предыдущим размером шага. В принципе, это должно давать погрешность около $0.9\times {\text{tol}}$ в следующей попытке. Если $|\tau _{n+1}^{(1)}|<{\text{tol}}$, считаем шаг успешным, а оценку ошибки используем для улучшения решения:

$y_{n+1}^{(2)}=y_{n+1}^{(1)}+\tau _{n+1}^{(1)}$

На самом деле это решение имеет третий порядок точности в локальной области (второй порядок в глобальной области), но поскольку для него нет оценки ошибки, это не помогает уменьшить количество шагов. Этот метод называется экстраполяцией Ричардсона .

Начиная с начального размера шага $h=b-a$, эта теория облегчает нашу контролируемую интеграцию ОДУ с точки $a$ к $b$, используя оптимальное количество шагов с учетом локальной толерантности к ошибкам. Недостаток состоит в том, что размер шага может стать непомерно малым, особенно при использовании метода Эйлера низкого порядка .

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

Встроенные оценки ошибок

Адаптивные методы размера шага, которые используют так называемую «встроенную» оценку ошибки, включают методы Богацкого-Шампина , Рунге-Кутты-Фельберга , Кэша-Карпа и Дормана-Принса . Эти методы считаются более эффективными в вычислительном отношении, но имеют меньшую точность оценок ошибок.

Чтобы проиллюстрировать идеи встроенного метода, рассмотрим следующую схему, которая обновляет $y_{n}$:

$y_{n+1}=y_{n}+h_{n}\psi (t_{n},y_{n},h_{n})$

$t_{n+1}=t_{n}+h_{n}$

Следующий шаг $h_{n}$ прогнозируется на основе предыдущей информации $h_{n}=g(t_{n},y_{n},h_{n-1})$.

Для встроенного метода РК вычисление $\psi$ включает метод RK низшего порядка ${\tilde {\psi }}$. Ошибкатогда можно просто записать как

${\textrm {err}}_{n}(h)={\tilde {y}}_{n+1}-y_{n+1}=h({\tilde {\psi }}(t_{n},y_{n},h_{n})-\psi (t_{n},y_{n},h_{n}))$

${\textrm {err}}_{n}$ это ненормированная ошибка. Чтобы его нормализовать, мы сравниваем его с заданным пользователем допуском, которыйсостоит из абсолютной толерантности и относительной толерантности:

${\textrm {tol}}_{n}={\textrm {Atol}}+{\textrm {Rtol}}\cdot \max(|y_{n}|,|y_{n-1}|)$

$E_{n}={\textrm {norm}}({\textrm {err}}_{n}/{\textrm {tol}}_{n})$

Затем сравниваем нормализованную ошибку $E_{n}$ против 1, чтобы получить предсказанное $h_{n}$:

$h_{n}=h_{n-1}(1/E_{n})^{1/(q+1)}$

Параметр q — порядок, соответствующий методу РК. ${\tilde {\psi }}$, который имеет меньший порядок.Приведенная выше формула прогнозирования правдоподобна в том смысле, что она увеличивает шаг, если предполагаемая локальная ошибка меньше, чемтолерантность, и в противном случае он сокращает шаг.

Приведенное выше описание представляет собой упрощенные процедуры, используемые при управлении размером шага для явных решателей RK. Более детальное лечениеможно найти в учебнике Хайрера. Решатель ОДУ во многих языках программирования использует эту процедуру в качестве стратегии по умолчанию для адаптивного управления размером шага, которая добавляет другие инженерные параметры, чтобы сделать систему более стабильной.

См. также

• Адаптивная квадратура
• Адаптивное числовое дифференцирование

Дальнейшее чтение

• Уильям Х. Пресс, Сол А. Теукольский, Уильям Т. Веттерлинг, Брайан П. Фланнери, Численные рецепты на языке C , второе издание, CAMBRIDGE UNIVERSITY PRESS, 1992. ISBN   0-521-43108-5
• Кендалл Э. Аткинсон, Численный анализ , второе издание, John Wiley & Sons, 1989. ISBN   0-471-62489-6