Дифференциальные уравнения с запаздывающим аргументом

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

Здесь предполагается, что вы ознакомились с руководством об обыкновенных дифференциальных уравнениях.

Дифференциальные уравнения с запаздыванием — это уравнения, имеющие запаздывающий аргумент. Чтобы задавать запаздывающий аргумент, определение функции для дифференциального уравнения с запаздыванием расширяется и включает в себя функцию истории h(p, t), которая использует интерполяции во всей истории решения для формирования непрерывного продолжения прошлого решателя и зависит от параметров p и времени t. Сигнатура функции для дифференциального уравнения с запаздыванием имеет вид f(u, h, p, t) для вычислений не на месте и f(du, u, h, p, t) для вычислений на месте.

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

Для приведенной задачи заметим, что является константой, и поэтому можно применять метод, использующий это поведение. Сначала напишем уравнение, используя соответствующую сигнатуру функции. Написание уравнений в основном выполняется одинаково, хотя мы используем функцию истории, сначала интерполируя, а затем выбирая компоненты. Таким образом, i-й компонент в момент времени t-tau задается с помощью h(p, t-tau)[i]. Компоненты без задержек записываются как в уравнении ODE.

Поэтому функция для данной модели имеет следующий вид:

using DifferentialEquations
function bc_model(du, u, h, p, t)
    p0, q0, v0, d0, p1, q1, v1, d1, d2, beta0, beta1, tau = p
    hist3 = h(p, t - tau)[3]
    du[1] = (v0 / (1 + beta0 * (hist3^2))) * (p0 - q0) * u[1] - d0 * u[1]
    du[2] = (v0 / (1 + beta0 * (hist3^2))) * (1 - p0 + q0) * u[1] +
            (v1 / (1 + beta1 * (hist3^2))) * (p1 - q1) * u[2] - d1 * u[2]
    du[3] = (v1 / (1 + beta1 * (hist3^2))) * (1 - p1 + q1) * u[2] - d2 * u[3]
end

bc_model (generic function with 1 method)

Теперь мы создаем DDEProblem. Сигнатура

prob = DDEProblem(f, u0, h, tspan, p = SciMLBase.NullParameters();
                  constant_lags = [], dependent_lags = [], kwargs...)

очень похожа на ODE, где теперь нужно задать запаздывания и функцию h. h — это функция истории, которая объявляет, какими были значения до момента запуска модели. Здесь мы будем считать, что за все время до t0 значения были равны 1, и определим h как функцию не на месте:

h(p, t) = ones(3)

h (generic function with 1 method)

Для использования модели с постоянным запаздыванием необходимо объявить запаздывания. Здесь мы будем использовать tau=1.

tau = 1
lags = [tau]

1-element Vector{Int64}:
 1

Далее выбираем временной диапазон (0.0,10.0) для решения и создаем тип задачи:

p0 = 0.2;
q0 = 0.3;
v0 = 1;
d0 = 5;
p1 = 0.2;
q1 = 0.3;
v1 = 1;
d1 = 1;
d2 = 1;
beta0 = 1;
beta1 = 1;
p = (p0, q0, v0, d0, p1, q1, v1, d1, d2, beta0, beta1, tau)
tspan = (0.0, 10.0)
u0 = [1.0, 1.0, 1.0]

prob = DDEProblem(bc_model, u0, h, tspan, p; constant_lags = lags)

DDEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 10.0)
u0: 3-element Vector{Float64}:
 1.0
 1.0
 1.0

Эффективным способом решения этой задачи (с учетом постоянных запаздываний) является использование решателя MethodOfSteps. Благодаря особенностям Julia метод решателя ODE OrdinaryDiffEq.jl преобразуется в метод дифференциальных уравнений с запаздыванием, который отличается высокой эффективностью благодаря преимуществам компилятора. Хорошим вариантом является метод пятого порядка Tsit5():

alg = MethodOfSteps(Tsit5())

MethodOfSteps{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, NLFunctional{Rational{Int64}, Rational{Int64}}, false}(Tsit5(; stage_limiter! = trivial_limiter!, step_limiter! = trivial_limiter!, thread = static(false),), NLFunctional{Rational{Int64}, Rational{Int64}}(1//100, 1//5, 10))

Для решения задач с малыми допусками можно с успехом использовать алгоритм BS3() (эта комбинация аналогична MATLAB dde23, но более эффективна), а для больших допусков алгоритм Vern6() выдаст решение шестого порядка.

Чтобы решить задачу с помощью этого алгоритма, мы делаем то же самое, что и с другими методами общего интерфейса:

sol = solve(prob, alg)

retcode: Success
Interpolation: specialized 4th order "free" interpolation
t: 40-element Vector{Float64}:
  0.0
  0.05836370738671701
  0.12916108569539597
  0.21568702493539388
  0.31476012203954973
  0.42647876001237145
  0.5486009950019299
  0.6801498148290146
  0.8195261166349427
  0.9657507075079192
  ⋮
  5.637340314503988
  6.0
  6.465321265034228
  6.9807287937908855
  7.578753298968584
  8.241987232937237
  8.98383276855844
  9.794183131169206
 10.0
u: 40-element Vector{Vector{Float64}}:
 [1.0, 1.0, 1.0]
 [0.7447276972023367, 0.9674847415959038, 0.9739922644849844]
 [0.520865690149123, 0.9216176741982295, 0.9429345997713768]
 [0.33647961749786465, 0.860710765286255, 0.9053954304799672]
 [0.20402080446565815, 0.7893156306081849, 0.8627636946755286]
 [0.11605349292143427, 0.710940127223295, 0.8151154583589272]
 [0.06263661020799872, 0.6308046470869514, 0.7637408792490028]
 [0.03223510074296827, 0.5524928215361826, 0.7096033440842167]
 [0.01594663432245627, 0.47891242866857714, 0.654067407727854]
 [0.007620742290428411, 0.41158157168839415, 0.5982822197402151]
 ⋮
 [1.9479598719208864e-12, 0.0025766300373552546, 0.0197374598472333]
 [3.3296189152665293e-13, 0.001729122637428431, 0.014435113877197158]
 [5.5064416255269584e-14, 0.0010364504244003874, 0.009606890211128613]
 [1.1838460929002524e-14, 0.0005879461631044232, 0.006079732552649997]
 [5.109127121132652e-15, 0.0003045520635702376, 0.003549648850830764]
 [4.075064652637717e-15, 0.0001468398691910801, 0.0019394744851741882]
 [6.623872359855281e-15, 6.493848650511485e-5, 0.0009786579729539678]
 [1.912280041308691e-14, 2.6637039533047875e-5, 0.0004599558275824244]
 [6.69983491663282e-15, 2.1240382256782126e-5, 0.0003792544387315653]

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

using Plots
plot(sol)

Ускорение интерполяций с помощью Idxs

Ускорить решение предыдущей задачи можно двумя различными способами. Прежде всего, если необходимо интерполировать несколько значений из предыдущего времени, можно использовать форму на месте для функции истории h(out, p, t), которая записывает вывод в out. В этом случае необходимо также задать начальные условия истории как и на месте. Для предыдущего примера это просто

h(out, p, t) = (out .= 1.0)

h (generic function with 2 methods)

и затем DDE имеет следующий вид:

const out = zeros(3) # Определяет переменную кэша
function bc_model(du, u, h, p, t)
    h(out, p, t - tau) # обновляется до корректной функции истории
    du[1] = (v0 / (1 + beta0 * (out[3]^2))) * (p0 - q0) * u[1] - d0 * u[1]
    du[2] = (v0 / (1 + beta0 * (out[3]^2))) * (1 - p0 + q0) * u[1] +
            (v1 / (1 + beta1 * (out[3]^2))) * (p1 - q1) * u[2] - d1 * u[2]
    du[3] = (v1 / (1 + beta1 * (out[3]^2))) * (1 - p1 + q1) * u[2] - d2 * u[3]
end

bc_model (generic function with 1 method)

Однако в большинстве случаев мы можем поступить еще хитрее. Нам потребовалось интерполировать прошлые значения только в индексе 3. Вместо того чтобы генерировать множество массивов, мы можем запросить конкретно это значение, передав ключевое слово idxs = 3. Теперь функция DDE bc_model имеет следующий вид:

function bc_model(du, u, h, p, t)
    u3_past_sq = h(p, t - tau; idxs = 3)^2
    du[1] = (v0 / (1 + beta0 * (u3_past_sq))) * (p0 - q0) * u[1] - d0 * u[1]
    du[2] = (v0 / (1 + beta0 * (u3_past_sq))) * (1 - p0 + q0) * u[1] +
            (v1 / (1 + beta1 * (u3_past_sq))) * (p1 - q1) * u[2] - d1 * u[2]
    du[3] = (v1 / (1 + beta1 * (u3_past_sq))) * (1 - p1 + q1) * u[2] - d2 * u[3]
end

bc_model (generic function with 1 method)

Обратите внимание, что для этого необходимо определить исторические значения

h(p, t; idxs = nothing) = typeof(idxs) <: Number ? 1.0 : ones(3)

h (generic function with 2 methods)

где idxs может быть целым числом, для которого нужно вычислить переменную в истории, и здесь для любого числа idxs мы возвращаем 1.0. Заметим, что если бы нужно было использовать прошлые значения i-й производной, в функции DDE мы бы вызвали функцию истории h(p, t, Val{i}) и должны были бы определить диспетчеризацию следующего вида:

h(p, t, ::Type{Val{1}}) = zeros(3)

h (generic function with 3 methods)

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

Обсуждение функциональных форм для функции истории также можно найти на странице DDEProblem.

Необъявленные и зависимые от состояния задержки посредством остаточного контроля

Возможно, вы заметили, что DifferentialEquations.jl позволяет решать задачи с необъявленными запаздываниями, поскольку интерполировать h можно при любом значении. Эту возможность следует использовать с осторожностью. Необъявленные запаздывания могут увеличить погрешность решения. При наличии необъявленных запаздываний рекомендуется применять метод с остаточным управлением, например MethodOfSteps(RK4()). С его помощью можно использовать интерполированные производные, решать функционально-дифференциальные уравнения с помощью квадратур в интерполянте и т. д. Однако следует иметь в виду, что остаточное управление выполняет решение с низким уровнем точности, поэтому допуски следует делать очень малыми, а доверять решению более чем на 2—3 знака после запятой не следует.

Примечание. Метод MethodOfSteps(RK4()) с необъявленными запаздываниями аналогичен ddesd MATLAB. Так, например, следующее аналогично решению примера выше из предыдущего раздела с остаточным управлением:

prob = DDEProblem(bc_model, u0, h, tspan)
alg = MethodOfSteps(RK4())
sol = solve(prob, alg)

retcode: Success
Interpolation: 3rd order Hermite
t: 56-element Vector{Float64}:
  0.0
  0.028686567733943524
  0.05496370797935472
  0.09262168203274358
  0.1328816015643946
  0.1804089117440632
  0.23156096928227454
  0.28745950203007636
  0.34649352838806813
  0.40865997317696995
  ⋮
  6.709304549678902
  7.077771808058965
  7.457378129664066
  7.849229500378771
  8.253885262007339
  8.671507191596193
  9.102489930712723
  9.548109990519311
 10.0
u: 56-element Vector{Vector{Float64}}:
 [1.0, 1.0, 1.0]
 [0.8651377233256287, 0.9847919107281841, 0.9871577813957303]
 [0.7576257549925667, 0.9695375729958094, 0.9754952861529256]
 [0.6264188253434527, 0.9459514608421642, 0.9589125735972589]
 [0.511173996718998, 0.9190777328442522, 0.9413127662016809]
 [0.4021017755857873, 0.8859139628064613, 0.9206598661690842]
 [0.31056686797828226, 0.8492819001921053, 0.8985424920191215]
 [0.23418948608094928, 0.8089369051538083, 0.8744786633816477]
 [0.17382236795861314, 0.766681887677747, 0.8491794697919689]
 [0.1269922932046665, 0.7231641434076941, 0.8226801184713568]
 ⋮
 [7.633678240249587e-14, 0.0007930624927795012, 0.007743143315919018]
 [2.290654525439246e-14, 0.00052886086662893, 0.005575060518185558]
 [7.192689065388156e-15, 0.000348386861669049, 0.00396242735079288]
 [2.3938297589461707e-15, 0.00022643491366957353, 0.0027774610777958984]
 [8.537747837490129e-16, 0.00014511825271079198, 0.001919152324510346]
 [3.290737503207264e-16, 9.169024318664183e-5, 0.0013070717660449717]
 [1.3830954331558912e-16, 5.7090180277528864e-5, 0.0008771566981009424]
 [6.431419123095559e-17, 3.498102566610629e-5, 0.0005793516695673948]
 [3.127939839039937e-17, 2.128729410362161e-5, 0.00037958433259645345]

Отметим, что данный метод позволяет решать задачи с запаздываниями, зависящими от состояния.

Зависимое от состояния отслеживание непрерывности задержек

Запаздывания, зависящие от состояния, представляют собой задачи, в которых запаздывание может быть функцией текущего состояния. Для их эффективного решения используется отслеживание прерываний. Для этого в DifferentialEquations.jl необходимо передать функции запаздывания g(u,p,t) в качестве ключевого слова dependent_lags определению DDEProblem. В остальном все одинаково, и эта задача решается с использованием общего интерфейса.

Мы можем решить описанную выше задачу с зависимым отслеживанием запаздывания, объявив зависимые запаздывания и решив ее с помощью алгоритма MethodOfSteps:

prob = DDEProblem(bc_model, u0, h, tspan; dependent_lags = ((u, p, t) -> tau,))
alg = MethodOfSteps(Tsit5())
sol = solve(prob, alg)

retcode: Success
Interpolation: specialized 4th order "free" interpolation
t: 31-element Vector{Float64}:
  0.0
  0.05836370738671701
  0.12916108569539597
  0.21568702493539388
  0.31476012203954973
  0.42647876001237145
  0.5486009950019299
  0.6801498148290146
  0.8195261166349427
  0.9657507075079192
  ⋮
  4.70131900781237
  5.319727380784704
  5.975944439730393
  6.672902478790312
  7.412222703550066
  8.203867606081308
  9.055215255977002
  9.965347765311543
 10.0
u: 31-element Vector{Vector{Float64}}:
 [1.0, 1.0, 1.0]
 [0.7447276972023367, 0.9674847415959038, 0.9739922644849844]
 [0.520865690149123, 0.9216176741982295, 0.9429345997713768]
 [0.33647961749786465, 0.860710765286255, 0.9053954304799672]
 [0.20402080446565815, 0.7893156306081849, 0.8627636946755286]
 [0.11605349292143427, 0.710940127223295, 0.8151154583589272]
 [0.06263661020799872, 0.6308046470869514, 0.7637408792490028]
 [0.03223510074296827, 0.5524928215361826, 0.7096033440842167]
 [0.01594663432245627, 0.47891242866857714, 0.654067407727854]
 [0.007620742290428411, 0.41158157168839415, 0.5982822197402151]
 ⋮
 [1.637789748146323e-9, 0.007211518284286898, 0.043270644077427596]
 [8.555301578118273e-10, 0.0036540032665950783, 0.025862653910404754]
 [6.387919765560543e-10, 0.0017756654863386229, 0.014739430025735469]
 [6.957421403133639e-10, 0.000825017766449486, 0.007996169982628906]
 [1.1062342089489479e-9, 0.0003658747495832752, 0.004126394925649491]
 [2.7421619004585237e-9, 0.0001531915519910343, 0.0020084764604955713]
 [1.0932031990743894e-8, 6.0070746916231975e-5, 0.0009160350421131426]
 [6.748471780530628e-8, 2.2069323013970547e-5, 0.00039185666579902566]
 [5.655280359438763e-8, 2.1246239705310403e-5, 0.00037932174234340285]

Здесь единичная задержка t-tau рассматривалась как запаздывание, зависящее от состояния. Разумеется, вы можете заменить этот кортеж функций любыми функциями, соответствующими вашим запаздываниям.