Зависящие от времени события в ODE

Question

Я недавно начал с Джулией и хотел реализовать одну из моих обычных проблем - реализовать зависящие от времени события.

На данный момент у меня есть:

# Packages
using Plots
using DifferentialEquations

# Parameters
k21 = 0.14*24
k12 = 0.06*24
ke = 1.14*24
α = 0.5
β = 0.05
η = 0.477
μ = 0.218
k1 = 0.5
V1 = 6

# Time
maxtime = 10
tspan = (0.0, maxtime)

# Dose
stim = 100

# Initial conditions
x0 = [0 0 2e11 8e11]

# Model equations
function system(dy, y, p, t)
  dy[1] = k21*y[2] - (k12 + ke)*y[1]
  dy[2] = k12*y[1] - k21*y[2]
  dy[3] = (α - μ - η)*y[3] + β*y[4] - k1/V1*y[1]*y[3]
  dy[4] = μ*y[3] - β*y[4]
end

# Events
eventtimes = [2, 5]
function condition(y, t, integrator)
    t - eventtimes
end
function affect!(integrator)
    x0[1] = stim
end
cb = ContinuousCallback(condition, affect!)

# Solve
prob = ODEProblem(system, x0, tspan)
sol = solve(prob, Rodas4(), callback = cb)

# Plotting
plot(sol, layout = (2, 2))

Но вывод, который дают, не является правильным. Более конкретно, события не принимаются во внимание, и начальное условие, по-видимому, не 0 для y1, а stim.

Любая помощь будет принята с благодарностью.

Answer 1

t - eventtimes не работает, потому что один скаляр, а другой вектор. Но для этого случая гораздо проще просто использовать DiscreteCallback. Когда вы делаете его DiscreteCallback, вы должны заранее установить время остановки так, чтобы оно совпадало с 2 и 5 для обратного вызова. Вот пример:

# Packages
using Plots
using DifferentialEquations

# Parameters
k21 = 0.14*24
k12 = 0.06*24
ke = 1.14*24
α = 0.5
β = 0.05
η = 0.477
μ = 0.218
k1 = 0.5
V1 = 6

# Time
maxtime = 10
tspan = (0.0, maxtime)

# Dose
stim = 100

# Initial conditions
x0 = [0 0 2e11 8e11]

# Model equations
function system(dy, y, p, t)
  dy[1] = k21*y[2] - (k12 + ke)*y[1]
  dy[2] = k12*y[1] - k21*y[2]
  dy[3] = (α - μ - η)*y[3] + β*y[4] - k1/V1*y[1]*y[3]
  dy[4] = μ*y[3] - β*y[4]
end

# Events
eventtimes = [2.0, 5.0]
function condition(y, t, integrator)
    t ∈ eventtimes
end
function affect!(integrator)
    integrator.u[1] = stim
end
cb = DiscreteCallback(condition, affect!)

# Solve
prob = ODEProblem(system, x0, tspan)
sol = solve(prob, Rodas4(), callback = cb, tstops = eventtimes)

# Plotting
plot(sol, layout = (2, 2))

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

В любом случае, обратите внимание, что affect был изменен на

function affect!(integrator)
    integrator.u[1] = stim
end

Необходимо изменить текущее значение u, иначе он ничего не будет делать.