Интерфейс интегратора

Интерфейс интегратора позволяет интерактивно пошагово выполнять численное решение дифференциального уравнения. С помощью этого интерфейса можно легко отслеживать результаты, изменять задачу во время выполнения и динамически продолжать решение по своему усмотрению.

Инициализация и пошаговый переход

Для инициализации интегратора используется следующий синтаксис:

integrator = init(prob, alg; kwargs...)

Принимаемые именованные аргументы совпадают с параметрами решателя, используемыми в solve, а возвращаемым значением является integrator, удовлетворяющее typeof(integrator)<:DEIntegrator. Можно вручную выбрать шаг с помощью команды step!:

step!(integrator)

что позволит сделать один успешный шаг. Кроме того:

step!(integrator, dt, stop_at_tdt = false)

при передаче dt интегратор продолжит выполнять шаги до integrator.t+dt, а при задании stop_at_tdt=true будет добавлен tstop, чтобы интегратор сделал шаг до integrator.t+dt.

Для проверки успешности выполнения шага интеграции можно вызвать функцию check_error(integrator), которая возвращает один из кодов возврата.

Этот тип также реализует интерфейс итератора, поэтому с помощью итератора take можно выполнять шаг n раз (или до последнего tstop):

for i in take(integrator, n)
end

Дойти до конца можно с помощью solve!(integrator) или используя интерфейс итератора:

for i in integrator
end

Кроме того, доступны некоторые вспомогательные итераторы, помогающие следить за ходом решения. Например, итератор tuples позволяет просматривать значения:

for (u, t) in tuples(integrator)
    @show u, t
end

а итератор intervals позволяет просматривать полный интервал:

for (uprev, tprev, u, t) in intervals(integrator)
    @show tprev, t
end

Кроме того, можно заставить итератор возвращать конкретные моменты времени с помощью TimeChoiceIterator:

ts = range(0, stop = 1, length = 11)
for (u, t) in TimeChoiceIterator(integrator, ts)
    @show u, t
end

Наконец, можно динамически управлять «конечной точкой». Инициализация просто делает prob.tspan[2] последним значением для tstop, а многие итераторы останавливаются на конечном значении tstop. Однако step! всегда будет делать шаг, и вы можете динамически добавлять новые значения tstops, изменяя переменную в поле параметров: add_tstop!(integrator,new_t).

Наконец, чтобы решить последний tstop, вызовите solve!(integrator). Выполнение init, а затем solve! эквивалентно solve.

# CommonSolve.step! — Function

step!(integ::DEIntegrator [, dt [, stop_at_tdt]])

Perform one (successful) step on the integrator.

Alternative, if a dt is given, then step! the integrator until there is a temporal difference ≥ dt in integ.t. When true is passed to the optional third argument, the integrator advances exactly dt.

# SciMLBase.check_error — Function

check_error(integrator)

Проверяет состояние integrator и возвращает один из кодов возврата.

# SciMLBase.check_error! — Function

check_error!(integrator)

То же, что и check_error, но также задает код возврата решения (integrator.sol.retcode) и выполняет postamble!.

Обработка интеграторов

Тип integrator<:DEIntegrator содержит всю информацию для промежуточного решения дифференциального уравнения. Полезными являются следующие поля:

t — время предлагаемого шага
u — значение на предлагаемом шаге
p — данные, предоставленные пользователем
opts — общие параметры решателя
alg — алгоритм, связанный с решением
f — решаемая функция
sol — текущее состояние решения
tprev — последняя точка времени
uprev — значение в последней точке времени
tdir — знак направления времени

Функция f обычно является оболочкой функции, предоставляемой при создании конкретной задачи. Например, при решении ODEProblem f будет функцией ODEFunction. Для доступа к функции справа, предоставленной пользователем при создании ODEProblem, используйте SciMLBase.unwrapped_f(integrator.f.f).

p — это данные (параметра), которые предоставляются пользователем в виде именованного аргумента в init. opts содержит все общие параметры решателя и может быть модифицирован для изменения характеристик решателя. Например, чтобы изменить абсолютный допуск для будущих временных интервалов, можно сделать следующее.

integrator.opts.abstol = 1e-9

Поле sol содержит текущее решение. Это текущее решение включает в себя функцию интерполяции, если она доступна, и, таким образом, integrator.sol(t) позволяет эффективно интерполировать все текущее решение. Кроме того, для типа integrator предусмотрена «функция интерполяции текущего интервала» с помощью integrator(t,deriv::Type=Val{0};idxs=nothing,continuity=:left). При этом для вычисления интерполяции используется только информация решателя из интервала [tprev,t], и допускается экстраполяция за пределы этого интервала.

Примечание об изменяемости

Будьте осторожны: не следует напрямую изменять поля t и u интегратора. В противном случае это приведет к снижению точности интерполятора и может нанести вред некоторым алгоритмам. Если необходимо ввести прерывающиеся изменения, следует использовать обратные вызовы. Изменения внутри обратного вызова affect!, окруженного saves, обеспечивают безошибочную обработку прерывания.

В качестве низкоуровневой альтернативы обратным вызовам можно использовать set_t!, set_u! и set_ut! для изменения состояний интегратора. Учтите, что у некоторых интеграторов могут отсутствовать эффективные способы изменения u и t. В этом случае set_*! столь же неэффективны, как и reinit!.

# SciMLBase.set_t! — Function

set_t!(integrator::DEIntegrator, t)

Устанавливает текущую временную точку integrator в значение t.

# SciMLBase.set_u! — Function

set_u!(integrator::DEIntegrator, u)
set_u!(integrator::DEIntegrator, sym, val)

Устанавливает текущее состояние integrator в значение u. Кроме того, может устанавливать состояние переменной sym в значение val.

# SciMLBase.set_ut! — Function

set_ut!(integrator::DEIntegrator, u, t)

Устанавливает текущее состояние integrator в значения u и t.

Интегратор и решение

Интегратор и решение имеют совершенно разные действия, поскольку у них совершенно разные значения. Тип typeof(sol) <: DESolution — это тип с историей: он хранит все (запрашиваемые) временные точки и интерполирует/действует, используя наиболее близкие по времени значения. В свою очередь, тип typeof(integrator)<:DEIntegrator является локальным объектом. Ему известны только время интервала, который он в данный момент охватывает, текущие кэши и значения, а также текущее состояние решателя (текущие параметры, допуски и т. д.). Они предназначены для совершенно разных целей:

Интерполяция интегратора (integrator) может экстраполировать как вперед, так и назад во времени. Это используется для оценки событий и внутренним образом применяется для прогнозирования.
Интегратор (integrator) является полностью изменяемым при итерации. Это означает, что каждый раз, когда используется аффект итератора, он будет выполнять временные шаги от текущего времени. Это означает, что first(integrator)!=first(integrator), поскольку интегратор (integrator) сделает один шаг для вычисления левой части и затем еще один шаг (не возвращаясь назад). Так итератор может продолжать динамически выполнять пошаговый переход, хотя следует отметить, что это может нарушить некоторые предположения о неизменяемости, обычно относящиеся к итераторам.

Если нужен объект решения, его можно найти в integrator.sol.

Интерфейс функции

В дополнение к интерфейсу типов предоставляется интерфейс функций, который позволяет безопасно изменять тип интегратора и обеспечивает единообразное использование во всей экосистеме (для пакетов/алгоритмов, реализующих функции). В состав интерфейса входят следующие функции:

Элементы управления сохранением

# SciMLBase.savevalues! — Function

savevalues!(integrator::DEIntegrator,
  force_save=false) -> Tuple{Bool, Bool}

Пытается сохранить переменные состояния и времени в текущей временной точке или точке saveat, при необходимости используя интерполяцию. Возвращает кортеж (saved, savedexactly). Если функция savevalues! сохранила значение, saved имеет значение true, а если функция savevalues! выполнила сохранение в текущей временной точке, savedexactly имеет значение true.

Действует следующий приоритет (порядок) сохранения.

save_on
- saveat
- force_save
- save_everystep

Кэши

# SciMLBase.get_tmp_cache — Function

get_tmp_cache(i::DEIntegrator)

Возвращает кортеж векторов внутреннего кэша, которые можно безопасно использовать как временные массивы. Следует применять для интерфейса интегратора и обратных вызовов, которым требуются массивы для записи, чтобы не выделять память. Длина кортежа зависит от метода.

# SciMLBase.full_cache — Function

full_cache(i::DEIntegrator)

Возвращается итератор по массивам кэша для метода. Можно использовать для изменения внутренних значений нужным образом.

Элементы управления пошаговым переходом

# SciMLBase.u_modified! — Function

u_modified!(i::DEIntegrator,bool)

Задает значение типа bool, которое указывает, имело ли место изменение переменной u, что позволяет решателю обрабатывать точку разрыва. По умолчанию имеет значение true, если используется обратный вызов. Приводит к повторному вычислению производной в точке t+dt. В этом нет необходимости, если применяется алгоритм FSAL и u не имеет разрыва непрерывности в конце интервала. Таким образом, если переменная u не изменяется в обратном вызове, единственный вызов для вычисления производной можно устранить посредством u_modified!(integrator,false).

# SciMLBase.get_proposed_dt — Function

get_proposed_dt(i::DEIntegrator)

Возвращает предлагаемое значение dt для следующего временного шага.

# SciMLBase.set_proposed_dt! — Function

set_proposed_dt(i::DEIntegrator,dt)
set_proposed_dt(i::DEIntegrator,i2::DEIntegrator)

Задает предлагаемое значение dt для следующего временного шага. Если вторым аргументом является DEIntegrator, задает для первого аргумента такие же временные шаги, как и для второго. Обратите внимание, что вследствие пропорционально-интегрального управления и ускорения шагов в большинстве случаев это не ограничивается сопоставлением множителей.

# SciMLBase.terminate! — Function

terminate!(i::DEIntegrator[, retcode = :Terminated])

Завершает выполнение интегратора, очищая tstops. Может использоваться в событиях и обратных вызовах для немедленного завершения процесса нахождения решения. Кроме того, можно указать retcode (см. раздел Коды возврата (RetCodes)).

# SciMLBase.change_t_via_interpolation! — Function

change_t_via_interpolation!(integrator::DEIntegrator,t,modify_save_endpoint=Val{false})

Изменяет текущее значение t и все соответствующие значения с использованием локальной интерполяции. Если текущее решение уже сохранено, можно указать необязательное значение modify_save_endpoint, чтобы также изменить конечную точку sol аналогичным образом.

# SciMLBase.add_tstop! — Function

add_tstop!(i::DEIntegrator,t)

Добавляет tstop в момент времени t.

# SciMLBase.has_tstop — Function

has_tstop(i::DEIntegrator)

Проверяет, определены ли в интеграторе моменты остановки.

# SciMLBase.first_tstop — Function

first_tstop(i::DEIntegrator)

Возвращает первый момент остановки интегратора.

# SciMLBase.pop_tstop! — Function

pop_tstop!(i::DEIntegrator)

Удаляет последний момент остановки из интегратора и возвращает его.

# SciMLBase.add_saveat! — Function

add_saveat!(i::DEIntegrator,t)

Добавляет временную точку saveat в t.

Изменение размера

# Base.resize! — Function

resize!(integrator::DEIntegrator,k::Int)

Изменяет размер ДУ на k. Конец массива отсекается, либо в его конец добавляются пустые значения в зависимости от условия k > length(integrator.u).

# Base.deleteat! — Function

deleteat!(integrator::DEIntegrator,idxs)

Сжимает ОДУ, удаляя компоненты idxs.

# SciMLBase.addat! — Function

addat!(integrator::DEIntegrator,idxs,val)

Расширяет ОДУ, добавляя компоненты idxs. Индексы должны быть непрерывными.

# SciMLBase.resize_non_user_cache! — Function

resize_non_user_cache!(integrator::DEIntegrator,k::Int)

Изменяет размер внутренних кэшей (недоступных пользователю) так, чтобы они были совместимы с ДУ размера k. В том числе изменяется размер кэшей якобианов.

Зачастую функция resize! просто изменяет размер переменных full_cache, а затем вызывает данную функцию. Для некоторых операций с AbstractArray требуется более точный контроль.

# SciMLBase.deleteat_non_user_cache! — Function

deleteat_non_user_cache!(integrator::DEIntegrator,idxs)

Удаляет внутренние кэши в позициях с индексами idxs с помощью функции deleteat!. В том числе изменяется размер кэшей якобианов.

Зачастую функция deleteat! просто удаляет (deleteat!) переменные full_cache, а затем вызывает данную функцию. Для некоторых операций с AbstractArray требуется более точный контроль.

# SciMLBase.addat_non_user_cache! — Function

addat_non_user_cache!(i::DEIntegrator,idxs)

Добавляет внутренние кэши в позициях с индексами idxs с помощью функции addat!. В том числе изменяется размер кэшей якобианов.

Зачастую функция addat! просто добавляет (addat!) переменные full_cache, а затем вызывает данную функцию. Для некоторых операций с AbstractArray требуется более точный контроль.

Повторная инициализация

# SciMLBase.reinit! — Function

reinit!(integrator::DEIntegrator,args...; kwargs...)

Функция reinit позволяет перезапустить интегрирование с нового значения.

Аргументы

u0: начальное значение u. Значение по умолчанию — integrator.sol.prob.u0.

Именованные аргументы

t0: начальная точка времени. Значение по умолчанию — integrator.sol.prob.tspan[1].
tf: конечная точка времени. Значение по умолчанию — integrator.sol.prob.tspan[2].
erase_sol=true: следует ли начать процесс без других значений в решении или оставить предыдущее решение.
tstops, d_discontinuities и saveat: кэш для хранения этих значений. По умолчанию используется исходный кэш.
reset_dt: указывает, следует ли сбросить текущее значение dt с помощью алгоритма автоматического определения dt. Значение по умолчанию — (integrator.dtcache == zero(integrator.dt)) && integrator.opts.adaptive
reinit_callbacks: указывает, следует ли инициализировать обратные вызовы повторно (для этого служит initialize_save). Значение по умолчанию — true.
reinit_cache: указывает, следует ли повторно выполнить функцию инициализации кэша (то есть сбросить FSAL без выделения памяти под векторы). Для получения правильного результата обычно требуется значение true. Значение по умолчанию — true.

Кроме того, можно выполнить функцию auto_dt_reset!, которая запускает алгоритм автоматической инициализации dt.

# SciMLBase.auto_dt_reset! — Function

auto_dt_reset!(integrator::DEIntegrator)

Запускает алгоритм автоматической инициализации dt.

Прочее

# SciMLBase.get_du — Function

get_du(i::DEIntegrator)

Возвращает производную в точке t.

# SciMLBase.get_du! — Function

get_du!(out,i::DEIntegrator)

Записывает текущую производную в точке t в out.

Обратите внимание, что не все эти функции будут реализованы для каждого алгоритма. Некоторые имеют жесткие ограничения. Например, Sundials.jl не может изменять размер задач. Если у функций нет ограничений, возникнет ошибка.

Дополнительные параметры

Для дальнейшего управления интегратором можно дополнительно указать следующие параметры в init (или изменить в opts).

advance_to_tstop: step! продолжит переход к следующему значению в tstop.
stop_at_next_tstop: итераторы будут останавливаться на следующем значении tstop.

Например, если требуется выполнять итерацию, но останавливаться только на определенных значениях, то можно выбрать следующее:

integrator = init(prob, Tsit5(); dt = 1 // 2^(4), tstops = [0.5], advance_to_tstop = true)
for (u, t) in tuples(integrator)
    @test t ∈ [0.5, 1.0]
end

когда вход в тело цикла будет выполняться только на значениях в tstops (здесь prob.tspan[2]==1.0 и, следовательно, есть два значения tstops для попадания). Кроме того, можно решить (solve!) только до 0.5 следующим образом:

integrator = init(prob, Tsit5(); dt = 1 // 2^(4), tstops = [0.5])
integrator.opts.stop_at_next_tstop = true
solve!(integrator)

Шаблон графика

Как и для типа DESolution, для типа DEIntegrator существует шаблон графика. Поскольку тип DEIntegrator является типом локального состояния в текущем интервале, plot(integrator) возвращает решение в текущем интервале. Для шаблона графика доступны те же параметры, что и для sol, то есть можно выбрать переменные с помощью именованного аргумента idxs, изменить plotdensity или включить/выключить denseplot.

Кроме того, поскольку integrator является итератором, его можно использовать в команде Plots.jl animate для итеративного построения анимации решения при решении дифференциального уравнения.

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

using DifferentialEquations, DiffEqProblemLibrary, Plots

# Линейное уравнение ODE, которое начинается в точке 0.5 и решается от t=0.0 до t=1.0
prob = ODEProblem((u, p, t) -> 1.01u, 0.5, (0.0, 1.0))

using Plots
integrator = init(prob, Tsit5(); dt = 1 // 2^(4), tstops = [0.5])
pyplot(show = true)
plot(integrator)
for i in integrator
    display(plot!(integrator, idxs = (0, 1), legend = false))
end
step!(integrator);
plot!(integrator, idxs = (0, 1), legend = false);
savefig("iteratorplot.png")