Описания методов использования временных шагов
Общая настройка
Все методы начинаются с вычисления масштабной расчетной погрешности для каждой скалярной компоненты :
По этой масштабной расчетной погрешности мы вычисляем норму. Обычно такой нормой является норма Хайрера:
Эта норма эффективна, так как не меняется при добавлении новых частей в дифференциальное уравнение: она масштабирует погрешность по количеству уравнений, так что шаги для независимых уравнений будут выполняться так же, как и при решении одного.
Во всех случаях шаг отклоняется, если , так как это означает, что погрешность больше допуска. Шаг принимается, если 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 снижается ниже теоретического оптимального значения.
Поскольку пропорциональное управление является неровным, т. е. может вызывать большие изменения между шагами, оно может повлиять на устойчивость явных методов. Таким образом, по умолчанию оно применяется только к неявным решателям низкого порядка.
#
OrdinaryDiffEq.IController — Type
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
qbeta1 — это коэффициент усиления пропорциональной части, а beta2 — коэффициент усиления исторической части. qoldinit — это инициализированное значение для истории коэффициента усиления.
#
OrdinaryDiffEq.PIController — Type
PIController(beta1, beta2)Пропорционально-интегральный (PI) контроллер является широко распространенным контроллером временного шага с улучшенными свойствами стабильности по сравнению с IController. Этот контроллер используется по умолчанию для большинства алгоритмов в OrdinaryDiffEq.jl.
Построим PI-контроллер размера шага, адаптирующий временной шаг, на основании формулы
Δtₙ₊₁ = εₙ₊₁^β₁ * εₙ^β₂ * Δtₙ
где εᵢ — это обратные величины оценок погрешностей, дифференцированные по допуску (Hairer, Nørsett, Wanner, 2010, Section IV.2). Коэффициент размера шага умножается на коэффициент запаса gamma и обрезается до интервала [qmin, qmax]. Шаг будет принят, если расчетная погрешность integrator.EEst меньше или равна единице. В противном случае шаг отклоняется, и выполняется повторная попытка с прогнозируемым размером шага.
|
Коэффициенты |
Справочные материалы
-
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
Пропорционально-интегрально-дифференциальный контроллер (ПИД-контроллер)
#
OrdinaryDiffEq.PIDController — Type
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 |
|---|---|---|---|
базовый |
|
|
|
PI42 |
|
|
|
PI33 |
|
|
|
PI34 |
|
|
|
H211PI |
|
|
|
H312PID |
|
|
|
|
В отличие от |
|
В отличие от других контроллеров, |
Справочные материалы
-
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#
OrdinaryDiffEq.PredictiveController — Type
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