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

Обработка событий и функции обратных вызовов

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

Типы обратных вызовов

Типы обратных вызовов определяются следующим образом. Существует три примитивных типа обратного вызова: ContinuousCallback, DiscreteCallback и VectorContinuousCallback:

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

  • DiscreteCallback применяется, когда его функция условия имеет значение true, но условие определяется только в конце каждого шага интеграции.

  • VectorContinuousCallback работает как вектор ContinuousCallbacks и позволяет пользователю задавать вектор непрерывных обратных вызовов, каждый из которых содержит уравнения одновременного поиска корней. Применяется эффект, соответствующий первому (самому раннему) выполненному условию. VectorContinuousCallback эффективнее набора CallbackSet обратных вызовов ContinuousCallbackпри увеличении числа обратных вызовов. Таким образом, это несколько более сложное определение, которое обеспечивает лучшее масштабирование.

ContinuousCallback

ContinuousCallback(condition,affect!,affect_neg!;
                   initialize = INITIALIZE_DEFAULT,
                   finalize = FINALIZE_DEFAULT,
                   idxs = nothing,
                   rootfind=LeftRootFind,
                   save_positions=(true,true),
                   interp_points=10,
                   abstol=10eps(),reltol=0,repeat_nudge=1//100)
function ContinuousCallback(condition,affect!;
                   initialize = INITIALIZE_DEFAULT,
                   finalize = FINALIZE_DEFAULT,
                   idxs = nothing,
                   rootfind=LeftRootFind,
                   save_positions=(true,true),
                   affect_neg! = affect!,
                   interp_points=10,
                   abstol=10eps(),reltol=0,repeat_nudge=1//100)

Содержит один обратный вызов, condition которого является непрерывной функцией. Обратный вызов срабатывает, когда результатом вычисления этой функции является 0.

Аргументы

  • condition. Это функция condition(u,t,integrator), используемая для объявления того, когда следует использовать обратный вызов. Обратный вызов инициируется, если условие достигает 0 в пределах временного интервала. Сведения об интеграторе (integrator) см. в документе Интерфейс интегратора.

  • affect!. Это функция affect!(integrator), позволяющая изменять текущее состояние интегратора. Если функция affect_neg! не передана, она вызывается, когда условие (condition) оказывается равным 0 (в корне), а пересечение является либо восходящим (от отрицательного к положительному), либо нисходящим (от положительного к отрицательному). Необходимо явно передать nothing в качестве аргумента affect_neg!, если он должен вызываться только на пересечениях, например ContinuousCallback(condition, affect!, nothing). Более подробная информация о том, что можно сделать, приведена на странице руководства Интерфейс интегратора. Изменения u в этой функции безопасны.

  • affect_neg!=affect!. Это функция affect_neg!(integrator), позволяющая изменять текущее состояние интегратора. Она вызывается, когда условие (condition) оказывается равным 0 (в корне), а пересечение является нисходящим (от положительного к отрицательному). Более подробная информация о том, что можно сделать, приведена на странице руководства Интерфейс интегратора. Изменения u в этой функции безопасны.

  • rootfind=LeftRootFind. Это флаг, указывающий тип поиска корней для поиска расположения события. Если задано значение LeftRootfind, решение будет возвращено в точку, где condition==0, и если решение не является точным, используется левый предел корня. Если задано значение RightRootFind, для решения будет задан правый предел корня. В противном случае системы и affect! произойдут в момент t+dt. Обратите внимание, что эти перечисления не экспортируются, и поэтому их нужно указывать как SciMLBase.LeftRootFind, SciMLBase.RightRootFind или SciMLBase.NoRootFind.

  • interp_points=10. Количество интерполированных точек для проверки условия. Поиск условия осуществляется путем проверки наличия другого знака у какой-либо точки интерполяции или конечной точки. Если interp_points=0, условия будут замечены только в том случае, если знак условия (condition) в момент t будет иным, чем в момент t+dt. Такое поведение не является надежным при осциллирующем характере решения, поэтому рекомендуется использовать некоторые точки интерполяции (их вычисление дешево). 0 в пределах временного интервала.

  • save_positions=(true,true). Логический кортеж для указания, выполнять ли сохранение до и после affect!. Это сохранение будет происходить непосредственно перед событием и после него, только во время события, и не зависит от таких параметров, как saveat, save_everystep и т. д. (т. е. если saveat=[1.0,2.0,3.0], точка сохранения все равно может быть добавлена в 2.1, если true). Для корректной (безошибочной) обработки прерывающихся изменений, таких как изменение u, необходимо задать save_positions=(true,true).

  • idxs=nothing. Компоненты, которые будут интерполированы в условие. Значение по умолчанию — nothing, что означает, что u будут все компоненты.

  • initialize. Это функция (c,u,t,integrator), которую можно использовать для инициализации состояния обратного вызова c. Она должна изменить аргумент c, а возврат игнорируется.

  • finalize. Это функция (c,u,t,integrator), которую можно использовать для финализации состояния обратного вызова c. Она может изменить аргумент c, а возврат игнорируется.

  • abstol=1e-14 и reltol=0. С их помощью задается отклонение от нуля для вычислителя корней: если начальное состояние меньше допуска от нуля, корень не будет обнаружен. Это необходимо, чтобы предотвратить возникновение повторных событий сразу после события поиска корней.

  • repeat_nudge = 1//100. Используется для установки следующей точки тестирования после ранее найденного нуля. По умолчанию установлено значение 1//100, т. е. после обратного вызова следующая проверка знака будет происходить не в момент t, а в момент t + dt*1//100, чтобы избежать повторов.

DiscreteCallback

DiscreteCallback(condition,affect!;
                 initialize = INITIALIZE_DEFAULT,
                 finalize = FINALIZE_DEFAULT,
                 save_positions=(true,true))

Аргументы

  • condition. Это функция condition(u,t,integrator), используемая для объявления того, когда следует использовать обратный вызов. Обратный вызов инициируется, если условие определяется как имеющее значение true. Сведения об интеграторе (integrator) см. в документе Интерфейс интегратора.

    • affect!: This is the function affect!(integrator) where one is allowed to

affect!. Это функция affect!(integrator), позволяющая изменять текущее состояние интегратора. Более подробная информация о том, что можно сделать, приведена на странице руководства [Интерфейс интегратора](../basics/integrator.md#integrator).

  • save_positions. Логический кортеж для указания, выполнять ли сохранение до и после affect!. Это сохранение будет происходить непосредственно перед событием и после него, только во время события, и не зависит от таких параметров, как saveat, save_everystep и т. д. (т. е. если saveat=[1.0,2.0,3.0], точка сохранения все равно может быть добавлена в 2.1, если true). Для корректной (безошибочной) обработки прерывающихся изменений, таких как изменение u, необходимо задать save_positions=(true,true).

  • initialize. Это функция (c,u,t,integrator), которую можно использовать для инициализации состояния обратного вызова c. Она должна изменить аргумент c, а возврат игнорируется.

  • finalize. Это функция (c,u,t,integrator), которую можно использовать для финализации состояния обратного вызова c. Она должна изменить аргумент c, а возврат игнорируется.

CallbackSet

struct CallbackSet{T1<:Tuple, T2<:Tuple} <: SciMLBase.DECallback

Несколько обратных вызовов могут быть объединены в набор CallbackSet. CallbackSet строится путем передачи конструктора ContinuousCallback, DiscreteCallback, VectorContinuousCallback или других экземпляров CallbackSet:

CallbackSet(cb1,cb2,cb3)

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

  • ContinuousCallback и VectorContinuousCallback применяются до DiscreteCallback (это связано с тем, что часто реализуется поиск событий, который будет выполнять поиск временного шага с возвратом к меньшему, чем dt).

  • Для ContinuousCallback и VectorContinuousCallback поиск времени событий осуществляется методом поиска корней, и применяется только первое воздействие ContinuousCallback или VectorContinuousCallback.

  • После этого DiscreteCallback применяются по порядку. Обратите внимание, что порядок имеет значение только для условий: если предыдущий обратный вызов изменяет u так, что следующий обратный вызов больше не определяет условие как имеющее значение true, его воздействие (affect) не будет применено.

VectorContinuousCallback

VectorContinuousCallback(condition,affect!,affect_neg!,len;
                         initialize = INITIALIZE_DEFAULT,
                         finalize = FINALIZE_DEFAULT,
                         idxs = nothing,
                         rootfind=LeftRootFind,
                         save_positions=(true,true),
                         interp_points=10,
                         abstol=10eps(),reltol=0,repeat_nudge = 1//100)
VectorContinuousCallback(condition,affect!,len;
                   initialize = INITIALIZE_DEFAULT,
                   finalize = FINALIZE_DEFAULT,
                   idxs = nothing,
                   rootfind=LeftRootFind,
                   save_positions=(true,true),
                   affect_neg! = affect!,
                   interp_points=10,
                   abstol=10eps(),reltol=0,repeat_nudge=1//100)

От также является подтипом AbstractContinuousCallback. Использовать CallbackSet нецелесообразно, если у вас много обратных вызовов, поскольку он плохо масштабируется. Для этого и существует VectorContinuousCallback — он позволяет иметь один обратный вызов для нескольких событий.

Аргументы

  • condition. Это функция condition(out, u, t, integrator), которая должна сохранять значение условия в массиве out по нужному индексу. Максимальный индекс out должен быть указан в свойстве len обратного вызова. Таким образом, можно получить цепочку len событий, которая при out[i] = 0 приведет к активации i-го события.

  • affect!. Это функция affect!(integrator, event_index), которая позволяет изменять интегратор (integrator) и сообщает, какое событие произошло, используя event_idx, т. е. выдает индекс i, для которого out[i] оказался равен нулю.

  • len. Количество объединенных в цепочку обратных вызовов. Это значение обязательно должно быть указано.

Остальные аргументы имеют то же значение, что и в ContinuousCallback.

Использование обратных вызовов

Тип обратного вызова затем передается решателю (или интегратору) через именованный аргумент callback:

sol = solve(prob, alg, callback = cb)

Можно указать nothing, один DiscreteCallback или ContinuousCallback, или VectorContinuousCallback, либо CallbackSet.

Примечание о сохранении

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

Изменение пошагового перехода в обратном вызове

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

Примеры вызовов DiscreteCallback

Пример 1. Вмешательства в заданные моменты времени

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

using DifferentialEquations
function f(du, u, p, t)
    du[1] = -u[1]
end
u0 = [10.0]
const V = 1
prob = ODEProblem(f, u0, (0.0, 10.0))
sol = solve(prob, Tsit5())
using Plots;
plot(sol);

Теперь предположим, что пациенту нужно дать дозу, равную 10, в момент времени t==4. Для этого можно использовать DiscreteCallback, который будет истинным только при t==4:

condition(u, t, integrator) = t == 4
affect!(integrator) = integrator.u[1] += 10
cb = DiscreteCallback(condition, affect!)
DiscreteCallback{typeof(Main.condition), typeof(Main.affect!), typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(SciMLBase.FINALIZE_DEFAULT)}(Main.condition, Main.affect!, SciMLBase.INITIALIZE_DEFAULT, SciMLBase.FINALIZE_DEFAULT, Bool[1, 1])

Если затем решить задачу с включенным этим обратным вызовом, никаких изменений не произойдет:

sol = solve(prob, Tsit5(), callback = cb)
plot(sol)

Причина отсутствия изменений заключается в том, что DiscreteCallback применяется только в определенное время, а интегратор никогда не достигает этого момента. Таким образом, чтобы применить условие, нам нужно сделать так, чтобы решатель ODE сделал шаг именно в момент t=4. Для этого можно воспользоваться аргументом tstops:

sol = solve(prob, Tsit5(), callback = cb, tstops = [4.0])
plot(sol)

и тем самым будет достигнут желаемый результат.

Тогда для применения многократных доз требуется лишь достигнуть нескольких точек. Например, для выдачи дозы в моменты времени t=4 и t=8 можно сделать следующее.

dosetimes = [4.0, 8.0]
condition(u, t, integrator) = t ∈ dosetimes
affect!(integrator) = integrator.u[1] += 10
cb = DiscreteCallback(condition, affect!)
sol = solve(prob, Tsit5(), callback = cb, tstops = dosetimes)
plot(sol)

Затем с помощью этого механизма можно произвольно усложнить модель. Например, допустим, сейчас существует 3 момента выдачи дозы, но доза выдается только в том случае, если текущая концентрация ниже 1,0. Кроме того, теперь доза составляет 10t, а не просто 10. Эта модель реализуется просто:

dosetimes = [4.0, 6.0, 8.0]
condition(u, t, integrator) = t ∈ dosetimes && (u[1] < 1.0)
affect!(integrator) = integrator.u[1] += 10integrator.t
cb = DiscreteCallback(condition, affect!)
sol = solve(prob, Tsit5(), callback = cb, tstops = dosetimes)
plot(sol)

PresetTimeCallback

Поскольку события в заданное время --это очень частое явление, в DifferentialEquations.jl предусмотрены предварительно встроенные обратные вызовы в библиотеке обратных вызовов. Функция PresetTimeCallback(tstops,affect!) принимает массив значений времени и применяемую функцию affect!. Таким образом, чтобы выполнить простой пример с двумя дозами с помощью этого обратного вызова, можно сделать следующее.

dosetimes = [4.0, 8.0]
affect!(integrator) = integrator.u[1] += 10
cb = PresetTimeCallback(dosetimes, affect!)
sol = solve(prob, Tsit5(), callback = cb)
plot(sol)

Обратите внимание, что в этой версии tstops задается автоматически.

Пример 2. Задача управления

Другим примером DiscreteCallback является переключение параметров задачи управления. Параметризованная система ODE имеет следующий вид:

Функция ODE будет использовать это поле следующим образом:

function f(du, u, p, t)
    du[1] = -0.5 * u[1] + p
    du[2] = -0.5 * u[2]
end
f (generic function with 1 method)

Теперь настроим механизм управления. Это будет простая настройка, использующая заданные моменты времени, в которые мы будем изменять p. В момент t=5.0 нужно увеличить значение p, а в момент t=8.0 — уменьшить значение p. Используя интерфейс DiscreteCallback, закодируем эти условия следующим образом.

const tstop1 = [5.0]
const tstop2 = [8.0]

function condition(u, t, integrator)
    t in tstop1
end

function condition2(u, t, integrator)
    t in tstop2
end
condition2 (generic function with 1 method)

Теперь нам нужно применить эффект, когда эти условия будут достигнуты. При достижении условия (condition) (при t=5.0) мы увеличим p до 1,5. При достижении условия (condition2) мы уменьшим p до -1.5. Это выполняется с помощью функций:

function affect!(integrator)
    integrator.p = 1.5
end

function affect2!(integrator)
    integrator.p = -1.5
end
affect2! (generic function with 1 method)

Используя эти функции, можно построить обратные вызовы:

using DifferentialEquations
save_positions = (true, true)

cb = DiscreteCallback(condition, affect!, save_positions = save_positions)

save_positions = (false, true)

cb2 = DiscreteCallback(condition2, affect2!, save_positions = save_positions)

cbs = CallbackSet(cb, cb2)
CallbackSet{Tuple{}, Tuple{DiscreteCallback{typeof(Main.condition), typeof(Main.affect!), typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(SciMLBase.FINALIZE_DEFAULT)}, DiscreteCallback{typeof(Main.condition2), typeof(Main.affect2!), typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(SciMLBase.FINALIZE_DEFAULT)}}}((), (DiscreteCallback{typeof(Main.condition), typeof(Main.affect!), typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(SciMLBase.FINALIZE_DEFAULT)}(Main.condition, Main.affect!, SciMLBase.INITIALIZE_DEFAULT, SciMLBase.FINALIZE_DEFAULT, Bool[1, 1]), DiscreteCallback{typeof(Main.condition2), typeof(Main.affect2!), typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(SciMLBase.FINALIZE_DEFAULT)}(Main.condition2, Main.affect2!, SciMLBase.INITIALIZE_DEFAULT, SciMLBase.FINALIZE_DEFAULT, Bool[0, 1])))

Теперь определим начальное условие. Начнем при [10.0;10.0] с p=0.0.

u0 = [10.0; 10.0]
p = 0.0
prob = ODEProblem(f, u0, (0.0, 10.0), p)
ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 10.0)
u0: 2-element Vector{Float64}:
 10.0
 10.0

Наконец, решаем задачу. Обратите внимание, что необходимо передать значения tstop, равные 5.0 и 8.0, чтобы решатель точно достиг этих временных точек:

const tstop = [5.0; 8.0]
sol = solve(prob, Tsit5(), callback = cbs, tstops = tstop)
using Plots;
plot(sol);

Из графика ясно, как механизм управления повлиял на результат.

Пример 3. AutoAbstol

В Simulink от MATLAB имеется возможность автоматического абсолютного допуска. В данном примере мы реализуем обратный вызов, который добавит это поведение в любой решатель JuliaDiffEq, реализующий интерфейс интегратора (integrator) и обратного вызова.

Алгоритм имеет следующий вид. Значение по умолчанию задано для запуска в 1e-6, хотя пользователь получит возможность выбора. Затем по мере выполнения моделирования на каждом шаге абсолютный допуск устанавливается равным максимальному значению, которое было достигнуто к этому моменту, умноженному на относительный допуск. Именно такое поведение мы реализуем в affect!.

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

condition = function (u, t, integrator)
    true
end
#1 (generic function with 1 method)

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

mutable struct AutoAbstolAffect{T}
    curmax::T
end
# Создайте `affect!` для этого:
function (p::AutoAbstolAffect)(integrator)
    p.curmax = max(p.curmax, integrator.u)
    integrator.opts.abstol = p.curmax * integrator.opts.reltol
    u_modified!(integrator, false)
end

affect!(integrator) будет использовать внутреннее изменяемое значение curmax для обновления абсолютного допуска интегратора согласно алгоритму.

Наконец, мы можем заключить его в небольшой конструктор:

function AutoAbstol(save = true; init_curmax = 1e-6)
    affect! = AutoAbstolAffect(init_curmax)
    condition = (u, t, integrator) -> true
    save_positions = (save, false)
    DiscreteCallback(condition, affect!, save_positions = save_positions)
end
AutoAbstol (generic function with 2 methods)

При этом будет создан DiscreteCallback из реализованных нами функций affect! и condition. Теперь

using DifferentialEquations
cb = AutoAbstol(true; init_curmax = 1e-6)
DiscreteCallback{Main.var"#4#5", Main.AutoAbstolAffect{Float64}, typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(SciMLBase.FINALIZE_DEFAULT)}(Main.var"#4#5"(), Main.AutoAbstolAffect{Float64}(1.0e-6), SciMLBase.INITIALIZE_DEFAULT, SciMLBase.FINALIZE_DEFAULT, Bool[1, 0])

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

function g(u, p, t)
    -u[1]
end
u0 = 10.0
const V = 1
prob = ODEProblem(g, u0, (0.0, 10.0))
integrator = init(prob, BS3(), callback = cb)
at1 = integrator.opts.abstol
step!(integrator)
at2 = integrator.opts.abstol
at1 <= at2
true
step!(integrator)
at3 = integrator.opts.abstol
at2 <= at3
true

Отметим, что данный пример содержится в библиотеке обратных вызовов, представляющей собой библиотеку полезных обратных вызовов для решателей JuliaDiffEq.

Примеры вызовов ContinuousCallback

Пример 1. Прыгающий мяч

Рассмотрим случай с прыгающим мячом. Пусть первая переменная y — это высота, изменяющаяся в зависимости от скорости v, где скорость всегда изменяется на -g, которая является гравитационной постоянной. Уравнение имеет следующий вид:

function f(du, u, p, t)
    du[1] = u[2]
    du[2] = -p
end
f (generic function with 1 method)

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

function condition(u, t, integrator) # Событие, когда condition(u,t,integrator) == 0
    u[1]
end
condition (generic function with 1 method)

Обратите внимание, что здесь мы использовали значения u, а не значения из integrator. Это связано с тем, что значения u,t будут соответствующим образом изменены в точках интерполяции, что позволит реализовать поведение поиска корней.

Теперь посмотрим, что нужно делать при наступлении события. В этом случае мы просто инверсируем скорость (вторую переменную).

function affect!(integrator)
    integrator.u[2] = -integrator.u[2]
end
affect! (generic function with 1 method)

Таким образом, обратный вызов задается следующим образом:

using DifferentialEquations
cb = ContinuousCallback(condition, affect!)
ContinuousCallback{typeof(Main.condition), typeof(Main.affect!), typeof(Main.affect!), typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(SciMLBase.FINALIZE_DEFAULT), Float64, Int64, Rational{Int64}, Nothing, Int64}(Main.condition, Main.affect!, Main.affect!, SciMLBase.INITIALIZE_DEFAULT, SciMLBase.FINALIZE_DEFAULT, nothing, SciMLBase.LeftRootFind, 10, Bool[1, 1], 1, 2.220446049250313e-15, 0, 1//100)

Затем можно выполнить решение и построить график:

u0 = [50.0, 0.0]
tspan = (0.0, 15.0)
p = 9.8
prob = ODEProblem(f, u0, tspan, p)
sol = solve(prob, Tsit5(), callback = cb)
using Plots;
plot(sol);

Как видно из полученного изображения, DifferentialEquations.jl использует интерполяцию для уточнения времени события и применяет событие в нужный момент. Таким образом, можно не беспокоиться о том, что адаптивная дискретизация по времени «проскочит» событие, так как это будет сделано за вас. Обратите внимание, что макрос события сохранит значение (значения) на разрыве.

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

u0 = [50.0, 0.0]
tspan = (0.0, 100.0)
prob = ODEProblem(f, u0, tspan, p)
sol = solve(prob, Tsit5(), callback = cb)
plot(sol)

Обработка изменяющейся динамики и точности

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

function dynamics!(du, u, p, t)
    du[1] = u[2]
    du[2] = p[1] * -9.8
end
function floor_aff!(int)
    int.p[1] = 0
    int.u[2] = 0
    @show int.u[1], int.t
end
floor_event = ContinuousCallback(condition, floor_aff!)
u0 = [1.0, 0.0]
p = [1.0]
prob = ODEProblem{true}(dynamics!, u0, (0.0, 1.75), p)
sol = solve(prob, Tsit5(), callback = floor_event)
plot(sol)

Обратите внимание, что в конце мяч находится не в точке 0.0, как можно было бы предположить из условия, а в точке 4.329177480185359e-16. Из результата внутри функции affect видно, что это значение, которое она имела в момент события t=0.4517539514526232. Почему после обработки событий значение не стало равно нулю? Если бы было выполнено моделирование по nextfloat(0.4517539514526232) = 0.45175395145262326, то мы бы увидели, что значение u[1] = -1.2647055847076505e-15. В этом можно убедиться, изменив аргумент rootfind обратного вызова:

floor_event = ContinuousCallback(condition, floor_aff!, rootfind = SciMLBase.RightRootFind)
u0 = [1.0, 0.0]
p = [1.0]
prob = ODEProblem{true}(dynamics!, u0, (0.0, 1.75), p)
sol = solve(prob, Tsit5(), callback = floor_event)
sol[end] # [-1.2647055847076505e-15, 0.0]
2-element Vector{Float64}:
 -5.9048760274925596e-15
  0.0

Это означает, что не существует 64-разрядного числа с плавающей запятой t, при котором условие будет равно нулю. По умолчанию, если не существует момента t, когда condition=0, вычислитель корней по умолчанию выбирает число с плавающей запятой точно перед событием (LeftRootFind). Таким образом, многочисленные ограничения сохраняются по умолчанию (т. е. мяч никогда не опускается ниже пола). Однако если требуется, чтобы условие выполнялось именно после наступления события, необходимо добавить такое изменение в функцию affect!. Например, исправление ошибки в данном случае заключается в добавлении int.u[1] = 0 к affect!, т. е:

function floor_aff!(int)
    int.p[1] = 0
    int.u[1] = 0
    int.u[2] = 0
    @show int.u[1], int.t
end
floor_event = ContinuousCallback(condition, floor_aff!)
u0 = [1.0, 0.0]
p = [1.0]
prob = ODEProblem{true}(dynamics!, u0, (0.0, 1.75), p)
sol = solve(prob, Tsit5(), callback = floor_event)
sol[end] # [0.0,0.0]
2-element Vector{Float64}:
 0.0
 0.0

и теперь поведение «приклеивания» идеально подходит для плавающей запятой.

Обработка точек накопления

Теперь рассмотрим прыгающий мяч в ситуации с трением. После отскока мы укажем скорость как -v/2. Поскольку скорость с каждым разом уменьшается вдвое, мы должны наблюдать поведение в стиле Зенона и видеть точку накопления отскоков. Мы будем использовать некоторые дополнительные параметры для подсчета количества отскоков (до бесконечности) и нахождения точки накопления. Давайте посмотрим.

function dynamics!(du, u, p, t)
    du[1] = u[2]
    du[2] = -9.8
end
function floor_aff!(int)
    int.u[2] *= -0.5
    int.p[1] += 1
    int.p[2] = int.t
end
floor_event = ContinuousCallback(condition, floor_aff!)
u0 = [1.0, 0.0]
p = [0.0, 0.0]
prob = ODEProblem{true}(dynamics!, u0, (0.0, 2.0), p)
sol = solve(prob, Tsit5(), callback = floor_event)
plot(sol)

Показатели говорят, что мяч отскочил всего восемь раз, прежде чем он опустился ниже уровня пола. Что произошло? Возникла ошибка с плавающей запятой. Поскольку нет уверенности, что существуют числа с плавающей запятой, позволяющие сделать condition=0, используется эвристика, гарантирующая, что нуль не будет случайно обнаружен в nextfloat(t) после перезапуска моделирования (в противном случае оно будет многократно находить одно и то же событие). Однако рано или поздно способность определять минутные различия с плавающей запятой даст сбой, и в результате, который должен быть бесконечно большим количеством отскоков, в итоге будет пропущен один отскок.

В связи с этим возникает два вопроса.

  1. Как можно повысить точность вычисления накопления?

  2. Как сделать так, чтобы оно продолжалось корректно?

Для (1) заметим, что точность вычислений с плавающей запятой зависит от текущего значения dt. Если известно, что приближается точка накопления, можно использовать set_proposed_dt! для уменьшения значения dt и поиска следующей точки отскока. С помощью t - tprev можно узнать длину предыдущего интервала для данного расчета. Для этого примера можно задать dt значение (t - tprev)/10, чтобы обеспечить постоянно возрастающую точность проверки.

Однако в какой-то момент мы попадем в машинный эпсилон, где t + eps(t) == t, поэтому не сможем измерить бесконечное множество отскоков и будем ограничены точностью представления чисел с плавающей запятой. Использование альтернативных числовых типов, таких как ArbFloats.jl, позволяет делать это с очень высокой точностью, но все же не бесконечно. Таким образом, необходимо определить допуск, после которого мы будем считать, что накопление достигнуто, и определить поведение завершения. В данном случае когда dt<1e-12, мы находимся почти на границе точности Float64 (eps(1.0) = 2.220446049250313e-16), поэтому изменим положение и скорость ровно на нуль.

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

function dynamics!(du, u, p, t)
    du[1] = u[2]
    du[2] = p[1] * -9.8
end
function floor_aff!(int)
    int.u[2] *= -0.5
    if int.dt > 1e-12
        set_proposed_dt!(int, (int.t - int.tprev) / 100)
    else
        int.u[1] = 0
        int.u[2] = 0
        int.p[1] = 0
    end
    int.p[2] += 1
    int.p[3] = int.t
end
floor_event = ContinuousCallback(condition, floor_aff!)
u0 = [1.0, 0.0]
p = [1.0, 0.0, 0.0]
prob = ODEProblem{true}(dynamics!, u0, (0.0, 2.0), p)
sol = solve(prob, Tsit5(), callback = floor_event)
plot(sol)

В этом исправленном варианте видно, что после 41 отскока точка накопления достигается при t = 1.355261854357056. Чтобы действительно увидеть накопление, рассмотрим этот момент более детально.

p1 = plot(sol, idxs = 1, tspan = (1.25, 1.40))
p2 = plot(sol, idxs = 1, tspan = (1.35, 1.36))
p3 = plot(sol, idxs = 1, tspan = (1.354, 1.35526))
p4 = plot(sol, idxs = 1, tspan = (1.35526, 1.35526185))
plot(p1, p2, p3, p4)

Думаю, Зенон гордился бы нашим решением.

Пример 2. Завершение интеграции

Часто требуется завершить интеграцию при выполнении некоторого условия. Для этого в обратном вызове используется terminate!(integrator) в качестве affect!.

В этом примере мы будем решать дифференциальное уравнение:

using DifferentialEquations
u0 = [1.0, 0.0]
function fun2(du, u, p, t)
    du[2] = -u[1]
    du[1] = u[2]
end
tspan = (0.0, 10.0)
prob = ODEProblem(fun2, u0, tspan)
ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 10.0)
u0: 2-element Vector{Float64}:
 1.0
 0.0

решениями которой являются косинус и -синус, соответственно. Мы будем решать задачу до тех пор, пока часть синуса u[2] не станет положительной. Мы можем искать две вещи.

Значение DiscreteCallback приведет к остановке на первом шаге, так как условие выполнено. Например, можно использовать следующее.

condition(u, t, integrator) = u[2] > 0
affect!(integrator) = terminate!(integrator)
cb = DiscreteCallback(condition, affect!)
sol = solve(prob, Tsit5(), callback = cb)
retcode: Terminated
Interpolation: specialized 4th order "free" interpolation
t: 13-element Vector{Float64}:
 0.0
 0.0009990009990009992
 0.010989010989010992
 0.07985922249873038
 0.2403882280626971
 0.48125583199780264
 0.7872724750204353
 1.1770558023021032
 1.630266513673432
 2.1040018658103854
 2.6513321708517577
 3.2715754793582668
 3.2715754793582668
u: 13-element Vector{Vector{Float64}}:
 [1.0, 0.0]
 [0.9999995009985435, -0.000999000832833342]
 [0.9999396214263468, -0.010988789821184267]
 [0.9968129466114029, -0.07977436592568171]
 [0.9712456180254766, -0.23807970964822986]
 [0.8864142948207253, -0.4628927349710369]
 [0.7057801319307301, -0.7084309070362929]
 [0.383645008268814, -0.9234805642321081]
 [-0.05943639522181044, -0.9982319800633571]
 [-0.5082985573738588, -0.86118077710699]
 [-0.882213099527795, -0.47085082387073085]
 [-0.9915642392403551, 0.12963031216471343]
 [-0.9915642392403551, 0.12963031216471343]

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

condition(u, t, integrator) = u[2]
affect!(integrator) = terminate!(integrator)
cb = ContinuousCallback(condition, affect!)
sol = solve(prob, Tsit5(), callback = cb)
using Plots;
plot(sol);

Заметим, что здесь используется поиск корней для приближенного определения «точного» момента пересечения. Аналитически мы знаем, что значением является pi, и интеграция завершается в точке

sol.t[end] # 3,1415902502224307
3.1415902497670927

Использование более точной интеграции повышает точность этого прогноза:

sol = solve(prob, Vern8(), callback = cb, reltol = 1e-12, abstol = 1e-12)
#π = 3,141592653589703...
sol.t[end] # 3,1415926535896035
3.141592653577113

Теперь предположим, что нужно найти момент окончания первого периода, т. е. будем игнорировать восходящее пересечение и останавливаемся только на нисходящем пересечении. Для этого мы игнорируем affect! и передаем affect! только для второго:

condition(u, t, integrator) = u[2]
affect!(integrator) = terminate!(integrator)
cb = ContinuousCallback(condition, nothing, affect!)
sol = solve(prob, Tsit5(), callback = cb)
plot(sol)

Обратите внимание, что передача только одной функции affect! аналогична передаче ContinuousCallback(condition,affect!,affect!), т. е. и восходящее, и нисходящие пересечения активируют событие. Таким образом, использование ContinuousCallback(condition,affect!,nothing) будет таким же, как и выше, поскольку первым событием является восходящее пересечение.

Пример 3. Рост популяции клеток

Другой интересный вопрос связан с моделями изменяющихся размеров. Возможность обработки таких событий является уникальной особенностью DifferentialEquations.jl. Здесь мы рассмотрим популяцию клеток. Мы начинаем с 1 клетки с белком X, которая линейно увеличивается со временем с параметром скорости α. Поскольку мы будем менять размер популяции, запишем модель в общем виде.

const α = 0.3
function f(du, u, p, t)
    for i in 1:length(u)
        du[i] = α * u[i]
    end
end
f (generic function with 1 method)

Модель состоит в том, что всякий раз, когда концентрация белка X достигает значения 1, запускается деление клетки. Поэтому мы проверяем, не достигла ли какая-либо концентрация единицы:

function condition(u, t, integrator) # Событие, когда condition(u,t,integrator) == 0
    1 - maximum(u)
end
condition (generic function with 1 method)

Опять же, вспомним, что эта функция находит события, когда condition==0, поэтому 1-maximum(u) положительна до тех пор, пока концентрация X в клетке не достигнет единицы, после чего и активируется событие. В момент события клетка делится на две клетки, выдавая каждой из новых клеток случайное количество белка. Это можно сделать, изменив размер кэша (прибавив 1 к длине всех кэшей) и задав значения этих двух клеток на момент события:

function affect!(integrator)
    u = integrator.u
    resize!(integrator, length(u) + 1)
    maxidx = findmax(u)[2]
    Θ = rand()
    u[maxidx] = Θ
    u[end] = 1 - Θ
    nothing
end
affect! (generic function with 1 method)

Как отмечается в документе Интерфейс интегратора, resize!(integrator,length(integrator.u)+1) используется для изменения длины всех внутренних кэшей (к которым относится и u) на их текущую длину + 1, что приводит к росту системы ODE. Затем следующий код задает новые концентрации белка. Теперь можно выполнить решение:

using DifferentialEquations
callback = ContinuousCallback(condition, affect!)
u0 = [0.2]
tspan = (0.0, 10.0)
prob = ODEProblem(f, u0, tspan)
sol = solve(prob, callback = callback)
retcode: Success
Interpolation: specialized 4th order "free" interpolation, specialized 2nd order "free" stiffness-aware interpolation
t: 17-element Vector{Float64}:
  0.0
  0.12735293592636002
  0.7397598129840375
  1.7483495595956313
  2.9466710993053042
  4.439321082686529
  5.36480663246384
  5.36480663246384
  6.029810268947888
  6.029810268947888
  7.487472001951411
  7.487472001951411
  8.109955749760646
  8.109955749760646
  9.489295456025763
  9.489295456025763
 10.0
u: 17-element Vector{Vector{Float64}}:
 [0.2]
 [0.20778902193773555]
 [0.24969628283119377]
 [0.33792440742897883]
 [0.4841131394025877]
 [0.757568052770951]
 [0.9999999999999998]
 [0.819139933567359, 0.18086006643264096]
 [0.9999999999999999, 0.2207926375228641]
 [0.354219472299632, 0.2207926375228641, 0.645780527700368]
 [0.5485137087688463, 0.3419004259993857, 0.9999999999999988]
 [0.5485137087688463, 0.3419004259993857, 0.1703446117757983, 0.8296553882242017]
 [0.6611343897167807, 0.41209932563831214, 0.2053197197217085, 0.9999999999999999]
 [0.6611343897167807, 0.41209932563831214, 0.2053197197217085, 0.5262502129944281, 0.4737497870055719]
 [0.9999999999999998, 0.6233215697868154, 0.31055670815984043, 0.795980697993739, 0.716571085053552]
 [0.17875277410805313, 0.6233215697868154, 0.31055670815984043, 0.795980697993739, 0.716571085053552, 0.8212472258919469]
 [0.20834910543035007, 0.7265257398580072, 0.3619759897621752, 0.9277722664408693, 0.8352146997053893, 0.9572222065113961]

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

using Plots
plot(sol.t, map((x) -> length(x), sol[:]), lw = 3,
     ylabel = "Number of Cells", xlabel = "Time")

Теперь проверим клетку. Мы все еще можем использовать интерполяцию для получения красивого графика концентрации клетки 1 с течением времени. Это делается с помощью команды:

ts = range(0, stop = 10, length = 100)
plot(ts, map((x) -> x[1], sol.(ts)), lw = 3,
     ylabel = "Amount of X in Cell 1", xlabel = "Time")

Обратите внимание, что при каждом достижении единицы клетка делится, предоставляя клетке 1 случайное количество X, которое затем растет до следующего деления.

Обратите внимание на один макрос, который не был показан в данном примере, — deleteat! для кэшей. Например, для удаления второй ячейки можно использовать следующее.

deleteat!(integrator, 2)

Это позволяет строить сложные модели популяций с рождаемостью и смертностью.

Пример вызова VectorContinuousCallback

Пример 1. Прыгающий мяч с несколькими стенами

Это похоже на приведенный выше пример с прыгающим мячом, но теперь у нас есть еще две вертикальные стены — x = 0 и x = 10.0. Функция ODEF имеет следующий вид:

function f(du, u, p, t)
    du[1] = u[2]
    du[2] = -p
    du[3] = u[4]
    du[4] = 0.0
end
f (generic function with 1 method)

где u[1] обозначает y-координату, u[2] — скорость в y-направлении, u[3] — x-координату и u[4] — скорость в x-направлении. Сделаем VectorContinuousCallback длиной 2: 1 для столкновения с осью x, 1 для стен, параллельных оси y.

function condition(out, u, t, integrator) # Событие, когда condition(out,u,t,integrator) == 0
    out[1] = u[1]
    out[2] = (u[3] - 10.0)u[3]
end

function affect!(integrator, idx)
    if idx == 1
        integrator.u[2] = -0.9integrator.u[2]
    elseif idx == 2
        integrator.u[4] = -0.9integrator.u[4]
    end
end
using DifferentialEquations
cb = VectorContinuousCallback(condition, affect!, 2)
VectorContinuousCallback{typeof(Main.condition), typeof(Main.affect!), typeof(Main.affect!), typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(SciMLBase.FINALIZE_DEFAULT), Float64, Int64, Rational{Int64}, Nothing, Int64}(Main.condition, Main.affect!, Main.affect!, 2, SciMLBase.INITIALIZE_DEFAULT, SciMLBase.FINALIZE_DEFAULT, nothing, SciMLBase.LeftRootFind, 10, Bool[1, 1], 1, 2.220446049250313e-15, 0, 1//100)

Очевидно, что out[2] будет равно нулю, если u[3] (координата x) равна либо 0.0, либо 10.0. И когда это происходит, мы инверсируем скорость с некоторым коэффициентом возмещения (0.9).

Завершение оставшейся части кода

u0 = [50.0, 0.0, 0.0, 2.0]
tspan = (0.0, 15.0)
p = 9.8
prob = ODEProblem(f, u0, tspan, p)
sol = solve(prob, Tsit5(), callback = cb, dt = 1e-3, adaptive = false)
using Plots;
plot(sol, idxs = (1, 3));