Документация Engee

Описания методов использования временных шагов

Общая настройка

Все методы начинаются с вычисления масштабной расчетной погрешности для каждой скалярной компоненты :

По этой масштабной расчетной погрешности мы вычисляем норму. Обычно такой нормой является норма Хайрера:

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

Во всех случаях шаг отклоняется, если , так как это означает, что погрешность больше допуска. Шаг принимается, если err^{scaled}<1.

Интегральный контроллер (стандартный контроллер)

Алгоритм пропорционального управления является стандартным алгоритмом для использования адаптивных временных шагов. Обратите внимание, что в DifferentialEquations.jl он не является заданным по умолчанию, так как это, как правило, плохо сказывается на производительности. Однако он объясняется первым, потому что является самым распространенным алгоритмом, на действии которого основываются другие алгоритмы.

Управление просто изменяет dt пропорционально погрешности. Существует экспонента, основанная на порядке алгоритма, которая восходит к результату Чечино для оптимального размера шага в целях уменьшения погрешности. Алгоритм имеет следующий вид:

qtmp = integrator.EEst^(1 / (alg_adaptive_order(integrator.alg) + 1)) /
       integrator.opts.gamma
@fastmath q = max(inv(integrator.opts.qmax), min(inv(integrator.opts.qmin), qtmp))
integrator.dtnew = integrator.dt / q

Таким образом, q — это масштабирующий множитель для dt, и он должен находиться в диапазоне от qmin до qmax. gamma — коэффициент запаса 0.9: на сколько dt снижается ниже теоретического оптимального значения.

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

IController()

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

Построим интегральный (I) контроллер размера шага, адаптирующий временной шаг, на основании формулы

Δtₙ₊₁ = εₙ₊₁^(1/k) * Δtₙ

где k = get_current_adaptive_order(alg, integrator.cache) + 1, а εᵢ — обратная величина оценки погрешности integrator.EEst, дифференцированная по допуску (Hairer, Nørsett, Wanner, 2008, Section II.4). Коэффициент размера шага умножается на коэффициент запаса gamma и обрезается до интервала [qmin, qmax]. Шаг будет принят, если расчетная погрешность integrator.EEst меньше или равна единице. В противном случае шаг отклоняется, и выполняется повторная попытка с прогнозируемым размером шага.

Справочные материалы

  • Hairer, Nørsett, Wanner (2008) Решение обыкновенных дифференциальных уравнений I Нежесткие задачи DOI:https://doi.org/10.1007/978-3-540-78862-1[]10.1007/978-3-540-78862-1

Пропорционально-интегральный контроллер (ПИ-контроллер)

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

EEst, beta1, q11, qold, beta2 = integrator.EEst, integrator.opts.beta1, integrator.q11,
                                integrator.qold, integrator.opts.beta2
@fastmath q11 = EEst^beta1
@fastmath q = q11 / (qold^beta2)
integrator.q11 = q11
@fastmath q = max(inv(integrator.opts.qmax),
                  min(inv(integrator.opts.qmin), q / integrator.opts.gamma))
if q <= integrator.opts.qsteady_max && q >= integrator.opts.qsteady_min
    q = one(q)
end
q

beta1 — это коэффициент усиления пропорциональной части, а beta2 — коэффициент усиления исторической части. qoldinit — это инициализированное значение для истории коэффициента усиления.

PIController(beta1, beta2)

Пропорционально-интегральный (PI) контроллер является широко распространенным контроллером временного шага с улучшенными свойствами стабильности по сравнению с IController. Этот контроллер используется по умолчанию для большинства алгоритмов в OrdinaryDiffEq.jl.

Построим PI-контроллер размера шага, адаптирующий временной шаг, на основании формулы

Δtₙ₊₁ = εₙ₊₁^β₁ * εₙ^β₂ * Δtₙ

где εᵢ — это обратные величины оценок погрешностей, дифференцированные по допуску (Hairer, Nørsett, Wanner, 2010, Section IV.2). Коэффициент размера шага умножается на коэффициент запаса gamma и обрезается до интервала [qmin, qmax]. Шаг будет принят, если расчетная погрешность integrator.EEst меньше или равна единице. В противном случае шаг отклоняется, и выполняется повторная попытка с прогнозируемым размером шага.

Коэффициенты beta1, beta2 не дифференцируются по порядку метода, в отличие от PIDController. Для PIController такое дифференцирование по порядку должно выполняться при построении контроллера.

Справочные материалы

  • Hairer, Nørsett, Wanner (2010) Решение обыкновенных дифференциальных уравнений II Жесткие и дифференциально-алгебраические задачи DOI:https://doi.org/10.1007/978-3-642-05221-7[]10.1007/978-3-642-05221-7

  • Hairer, Nørsett, Wanner (2008) Решение обыкновенных дифференциальных уравнений I Нежесткие задачи DOI:https://doi.org/10.1007/978-3-540-78862-1[]10.1007/978-3-540-78862-1

Пропорционально-интегрально-дифференциальный контроллер (ПИД-контроллер)

PIDController(beta1, beta2, beta3=zero(beta1);
              limiter=default_dt_factor_limiter,
              accept_safety=0.81)

Пропорционально-интегральный дифференциальный (PID) контроллер является обобщением PIController и может обладать улучшенными свойствами устойчивости и эффективности.

Построим PID-контроллер размера шага, адаптирующий временной шаг, на основании формулы

Δtₙ₊₁ = εₙ₊₁^(β₁/k) * εₙ^(β₂/k) * εₙ₋₁^(β₃/ k) * Δtₙ

где k = min(alg_order, alg_adaptive_order) + 1, и εᵢ являются обратными величинами оценок погрешностей, дифференцированными по допуску (Söderlind, 2003). Коэффициент размера шага ограничивается с помощью limiter со значением по умолчанию,

limiter(x) = one(x) + atan(x - one(x))

в соответствии с предложением Сёдерлинда (Söderlind) и Ванга (Wang) (2006). Шаг будет принят, если изменение прогнозируемого размера шага превышает accept_safety. В противном случае шаг отклоняется, и выполняется повторная попытка с прогнозируемым размером шага.

Ниже приводятся некоторые стандартные параметры контроллера, предлагаемые в документации.

Контроллер beta1 beta2 beta3

базовый

1.00

0.00

0

PI42

0.60

-0.20

0

PI33

2//3

-1//3

0

PI34

0.70

-0.40

0

H211PI

1//6

1//6

0

H312PID

1//18

1//9

1//18

В отличие от PIController, коэффициенты beta1, beta2, beta3 дифференцируются по порядку метода. Таким образом, стандартные контроллеры, например PI42, могут использовать одни и те же коэффициенты beta1, beta2, beta3 для разных алгоритмов.

В отличие от других контроллеров, PIDController не использует именованные аргументы qmin, qmax для ограничения изменения временного шага или коэффициента запаса gamma. Эти общие именованные аргументы заменяются с помощью limiter и accept_safety для гарантии корректного поведения (Сёдерлинд и Ванг, 2006). В связи с этим работа PIDController отличается от работы PIController, даже если beta1, beta2 адаптированы соответствующим образом и iszero(beta3).

Справочные материалы

  • Söderlind (2003) Digital Filters in Adaptive Time-Stepping DOI:https://doi.org/10.1145/641876.641877[]10.1145/641876.641877

  • Söderlind, Wang (2006) Adaptive time-stepping and computational stability DOI:https://doi.org/10.1016/j.cam.2005.03.008[]10.1016/j.cam.2005.03.008

  • Ranocha, Dalcin, Parsani, Ketcheson (2021) Optimized Runge-Kutta Methods with Automatic Step Size Control for Compressible Computational Fluid Dynamics arXiv:2104.06836

Ускорение Густавсона

Алгоритм ускорения Густафссона ускоряет изменения, чтобы алгоритмы могли быстрее меняться для обработки быстрых переходных процессов. Таким образом, этот алгоритм хорошо подходит для жестких решателей, где можно ожидать подобные ситуации, и является стандартным для алгоритмов, подобных методам (E)SDIRK.

gamma = integrator.opts.gamma
niters = integrator.cache.newton_iters
fac = min(gamma,
          (1 + 2 * integrator.alg.max_newton_iter) * gamma /
          (niters + 2 * integrator.alg.max_newton_iter))
expo = 1 / (alg_order(integrator.alg) + 1)
qtmp = (integrator.EEst^expo) / fac
@fastmath q = max(inv(integrator.opts.qmax), min(inv(integrator.opts.qmin), qtmp))
if q <= integrator.opts.qsteady_max && q >= integrator.opts.qsteady_min
    q = one(q)
end
integrator.qold = q
q

В данном случае niters — это количество итераций Ньютона, которое потребовалось на последнем шаге алгоритма. Обратите внимание, что эти значения используются по-разному в зависимости от принятия и отклонения. Если шаг принят, применяется следующая логика:

if integrator.success_iter > 0
    expo = 1 / (alg_adaptive_order(integrator.alg) + 1)
    qgus = (integrator.dtacc / integrator.dt) *
           (((integrator.EEst^2) / integrator.erracc)^expo)
    qgus = max(inv(integrator.opts.qmax),
               min(inv(integrator.opts.qmin), qgus / integrator.opts.gamma))
    qacc = max(q, qgus)
else
    qacc = q
end
integrator.dtacc = integrator.dt
integrator.erracc = max(1e-2, integrator.EEst)
integrator.dt / qacc

При отклонении происходит то же самое, что и при пропорциональном управлении:

if integrator.success_iter == 0
    integrator.dt *= 0.1
else
    integrator.dt = integrator.dt / integrator.qold
end
PredictiveController()

Алгоритм ускорения Густафссона ускоряет изменения, чтобы алгоритмы могли быстрее меняться для обработки быстрых переходных процессов. Таким образом, этот алгоритм хорошо подходит для жестких решателей, где можно ожидать подобные ситуации, и является стандартным для алгоритмов, подобных методам (E)SDIRK.

gamma = integrator.opts.gamma
niters = integrator.cache.newton_iters
fac = min(gamma,(1+2*integrator.alg.max_newton_iter)*gamma/(niters+2*integrator.alg.max_newton_iter))
expo = 1/(alg_order(integrator.alg)+1)
qtmp = (integrator.EEst^expo)/fac
@fastmath q = max(inv(integrator.opts.qmax),min(inv(integrator.opts.qmin),qtmp))
if q <= integrator.opts.qsteady_max && q >= integrator.opts.qsteady_min
  q = one(q)
end
integrator.qold = q
q

В данном случае niters — это количество итераций Ньютона, которое потребовалось на последнем шаге алгоритма. Обратите внимание, что эти значения используются по-разному в зависимости от принятия и отклонения. Если шаг принят, применяется следующая логика.

if integrator.success_iter > 0
  expo = 1/(alg_adaptive_order(integrator.alg)+1)
  qgus=(integrator.dtacc/integrator.dt)*(((integrator.EEst^2)/integrator.erracc)^expo)
  qgus = max(inv(integrator.opts.qmax),min(inv(integrator.opts.qmin),qgus/integrator.opts.gamma))
  qacc=max(q,qgus)
else
  qacc = q
end
integrator.dtacc = integrator.dt
integrator.erracc = max(1e-2,integrator.EEst)
integrator.dt/qacc

При отклонении происходит то же самое, что и при IController:

if integrator.success_iter == 0
  integrator.dt *= 0.1
else
  integrator.dt = integrator.dt/integrator.qold
end