From cd7f71e6529b801e08390861fb60f60711b384dc Mon Sep 17 00:00:00 2001 From: David Widmann Date: Mon, 4 Sep 2017 02:50:05 +0200 Subject: [PATCH] Improve and add DDE problems --- src/DiffEqProblemLibrary.jl | 12 +- src/dde_premade_problems.jl | 571 ++++++++++++++++++++++++++++-------- 2 files changed, 461 insertions(+), 122 deletions(-) diff --git a/src/DiffEqProblemLibrary.jl b/src/DiffEqProblemLibrary.jl index f6da6b3..d9db11c 100644 --- a/src/DiffEqProblemLibrary.jl +++ b/src/DiffEqProblemLibrary.jl @@ -26,12 +26,20 @@ export prob_ode_linear, prob_ode_bigfloatlinear, prob_ode_2Dlinear, #DAE Example Problems export prob_dae_resrob -#DDE Example Problems +# DDE Example Problems +# examples with constant delays export prob_dde_1delay, prob_dde_1delay_notinplace, prob_dde_1delay_scalar_notinplace, prob_dde_2delays, prob_dde_2delays_notinplace, prob_dde_2delays_scalar_notinplace, prob_dde_1delay_long, prob_dde_1delay_long_notinplace, prob_dde_1delay_long_scalar_notinplace, prob_dde_2delays_long, - prob_dde_2delays_long_notinplace, prob_dde_2delays_long_scalar_notinplace + prob_dde_2delays_long_notinplace, prob_dde_2delays_long_scalar_notinplace, + prob_dde_mackey, prob_dde_wheldon, prob_dde_qs, +# examples with vanishing time dependent delays + prob_ddde_neves1, prob_dde_neves_thompson, +# examples with state dependent delays + prob_dde_paul1, prob_dde_paul2, prob_dde_mahaffy1, prob_dde_mahaffy2, +# examples with vanishing state dependent delays + prob_neves2, prob_dde_gatica #FEM Example Problems export prob_femheat_moving, prob_femheat_pure, prob_femheat_diffuse, diff --git a/src/dde_premade_problems.jl b/src/dde_premade_problems.jl index 7de1cc2..7cf47c0 100644 --- a/src/dde_premade_problems.jl +++ b/src/dde_premade_problems.jl @@ -4,47 +4,56 @@ ### In-place function -f_1delay = function (t,u,h,du) - du[1] = - h(t-1)[1] +function f_1delay(t, u, h, du) + du[1] = - h(t - oneunit(t))[1] / oneunit(t) end -function (f::typeof(f_1delay))(::Type{Val{:analytic}}, t, u₀) - if t < 0 - return 0 - elseif t < 1 +function f_1delay(::Type{Val{:analytic}}, t, u₀) + z = t/oneunit(t) + + if z < 0 + return zero(u₀) + elseif z < 1 return u₀ - elseif t < 2 - return u₀ * (2 - t) - elseif t < 3 - return u₀ * (8 - 6t + t^2) / 2 - elseif t < 4 - return u₀ * (51 - 45t + 12t^2 - t^3) / 6 - elseif t < 5 - return u₀ * (460 - 436t + 144t^2 - 20t^3 + t^4) / 24 - elseif t < 6 - return u₀ * (5_425 - 5_305t + 1_970t^2 - 350t^3 + 30t^4 - t^5) / 120 - elseif t < 7 - return u₀ * (79_206 - 78_486t + 31_260t^2 - 6_420t^3 + 720t^4 - 42t^5 + t^6) / 720 - elseif t < 8 - return u₀ * (1_377_985 - 1_372_945t + 571_767t^2 - 128_975t^3 + 17_045t^4 - - 1_323t^5 + 56t^6 - t^7) / 5040 - elseif t < 9 - return u₀ * (27_801_096 - 27_760_776t + 11_914_168t^2 - 2_866_808t^3 + 423_080t^4 - - 39_256t^5 + 2_240t^6 - 72t^7 + t^8) / 40_320 - elseif t ≤ 10 - return u₀ * (637_630_353 - 637_267_473t + 279_414_396t^2 - 70_442_316t^3 + - 11_247_894t^4 - 1_179_990t^5 + 81_396t^6 - 3_564t^7 + 90t^8 - t^9) / - 362_880 + elseif z ≤ 10 + if z < 2 + c = @evalpoly(z, 2, -1) + elseif z < 3 + c = @evalpoly(z, 4, -3, 1//2) + elseif z < 4 + c = @evalpoly(z, 17//2, -15//2, 2, -1//6) + elseif z < 5 + c = @evalpoly(z, 115//6, -109//6, 6, -5//6, 1//24) + elseif z < 6 + c = @evalpoly(z, 1085//24, -1061//24, 197//12, -35//12, 1//4, -1//120) + elseif z < 7 + c = @evalpoly(z, 13201//120, -13081//120, 521//12, -107//12, 1, -7//120, 1//720) + elseif z < 8 + c = @evalpoly(z, 39371//144, -39227//144, 27227//240, -3685//144, 487//144, + -21//80, 1//90, -1//5040) + elseif z < 9 + c = @evalpoly(z, 1158379//1680, -1156699//1680, 212753//720, -51193//720, + 1511//144, -701//720, 1//18, -1//560, 1//40320) + else + c = @evalpoly(z, 23615939//13440, -23602499//13440, 7761511//10080, + -279533//1440, 89269//2880, -1873//576, 323//1440, -11//1120, + 1//4032, -1//362880) + end + + return c .* u₀ else error("This analytical solution is only valid on (-∞,10]") end end +build_prob_dde_1delay(u₀, ::T=u₀) where {T} = + DDEProblem(f_1delay, t->[zero(u₀)], [u₀], (zero(T), T(10)), [oneunit(T)]) + """ - prob_dde_1delay(u₀) + prob_dde_1delay -Model problem of finding a solution ``u(t)`` in the time span ``t \\in [0,10]`` to the -delay differential equation +Model problem of finding a solution ``u(t)`` in the time span ``t \\in [0,10]`` to the delay +differential equation ```math \\frac{du}{dt} = -u(t-1) @@ -59,86 +68,99 @@ u₀ & \\text{if } t = 0, \\end{cases} ``` -for ``t \\leq 0``. Hence the problem is discontinuous at ``t = 0`` for all ``u₀ \\neq 0``. +for ``t \\leq 0`` with ``u₀ = 1``. Note that the problem is discontinuous at ``t = 0`` for +all ``u₀ \\neq 0``. An analytical solution of this problem is provided for ``t \\in (-\\infty,10]``. """ -prob_dde_1delay(u₀) = DDEProblem(f_1delay, t->[0.0], [u₀], (0.0, 10.0), [1]) +prob_dde_1delay = build_prob_dde_1delay(1.0) ### Not in-place function -f_1delay_notinplace = function (t,u,h) - - h(t-1) +function f_1delay_notinplace(t, u, h) + - h(t - oneunit(t)) ./ oneunit(t) end -(f::typeof(f_1delay_notinplace))(::Type{Val{:analytic}}, t, u0) = f_1delay(Val{:analytic}, - t, u0) +f_1delay_notinplace(::Type{Val{:analytic}}, t, u0) = f_1delay(Val{:analytic}, t, u0) #### Vectorized history function +build_prob_dde_1delay_notinplace(u₀, ::T=u₀) where {T} = + DDEProblem(f_1delay_notinplace, t->[zero(u₀)], [u₀], (zero(T), T(10)), [oneunit(T)]) + """ - prob_dde_1delay_notinplace(u₀) + prob_dde_1delay_notinplace Same as [`prob_dde_1delay`](@ref), but purposefully implemented with a not in-place function. """ -prob_dde_1delay_notinplace(u₀) = - DDEProblem(f_1delay_notinplace, t->[0.0], [u₀], (0.0, 10.0), [1]) +prob_dde_1delay_notinplace = build_prob_dde_1delay_notinplace(1.0) #### Scalar history function +build_prob_dde_1delay_scalar_notinplace(u₀, ::T=u₀) where {T} = + DDEProblem(f_1delay_notinplace, t -> zero(u₀), u₀, (zero(T), T(10)), [oneunit(T)]) + """ - prob_dde_1delay_scalar_notinplace(u₀) + prob_dde_1delay_scalar_notinplace Same as [`prob_dde_1delay_notinplace`](@ref), but purposefully implemented with a scalar history function. """ -prob_dde_1delay_scalar_notinplace(u₀) = - DDEProblem(f_1delay_notinplace, t->0.0, u₀, (0.0, 10.0), [1]) +prob_dde_1delay_scalar_notinplace = build_prob_dde_1delay_scalar_notinplace(1.0) ## Two constant delays ### In-place function -f_2delays = function (t,u,h,du) - du[1] = - h(t-1/3)[1] - h(t-1/5)[1] +function f_2delays(t::T, u, h, du) where {T<:Number} + du[1] = (- h(t - T(1//3))[1] - h(t - T(1//5))[1]) / oneunit(t) end -function (f::typeof(f_2delays))(::Type{Val{:analytic}}, t, u₀) - if t < 0 - return 0 - elseif t < 1/5 +function f_2delays(::Type{Val{:analytic}}, t, u₀) + z = t/oneunit(t) + + if z < 0 + return zero(u₀) + elseif z < 1/5 return u₀ - elseif t < 1/3 - return u₀ * (6 - 5t) / 5 - elseif t < 2/5 - return u₀ * (23 - 30t) / 15 - elseif t < 8/15 - return u₀ * (242 - 360t + 75t^2) / 150 - elseif t < 3/5 - return u₀ * (854 - 1_560t + 675t^2) / 450 - elseif t < 2/3 - return u₀ * (4_351 - 8_205t + 4_050t^2 - 375t^3) / 2_250 - elseif t < 11/15 - return u₀ * (1_617 - 3_235t + 1_725t^2 - 125t^3) / 750 - elseif t < 4/5 - return u₀ * (7_942 - 17_280t + 11_475t^2 - 2_250t^3) / 3_375 - elseif t < 13/15 - return u₀ * (319_984 - 702_720t + 480_600t^2 - 108_000t^3 + 5_625t^4) / 135_000 - elseif t < 14/15 - return u₀ * (40_436 - 94_980t + 72_900t^2 - 19_500t^3 + 625t^4) / 15_000 - elseif t ≤ 1 - return u₀ * (685_796 - 1_670_388t + 1_392_660t^2 - 467_100t^3 + 50_625t^4) / 243_000 + elseif z ≤ 1 + if z < 1/3 + c = @evalpoly(z, 6//5, -1) + elseif z < 2/5 + c = @evalpoly(z, 23//15, -2) + elseif z < 8/15 + c = @evalpoly(z, 121//75, -12//5, 1//2) + elseif z < 3/5 + c = @evalpoly(z, 427//225, -52//15, 3//2) + elseif z < 2/3 + c = @evalpoly(z, 4351//2250, -547//150, 9//5, -1//6) + elseif z < 11/15 + c = @evalpoly(z, 539//250, -647//150, 23//10, -1//6) + elseif z < 4/5 + c = @evalpoly(z, 7942//3375, -128//25, 17//5, -2//3) + elseif z < 13/15 + c = @evalpoly(z, 39998//16875, -1952//375, 89//25, -4//5, 1//24) + elseif z < 14/15 + c = @evalpoly(z, 10109//3750, -1583//250, 243//50, -13//10, 1//24) + else + c = @evalpoly(z, 171449//60750, -139199//20250, 2579//450, -173//90, 5//24) + end + + return c .* u₀ else error("This analytical solution is only valid on (-∞,1]") end end +build_prob_dde_2delays(u₀, ::T=u₀) where {T} = + DDEProblem(f_2delays, t -> [zero(u₀)], [u₀], (zero(T), oneunit(T)), [T(1//3), T(1//5)]) + """ - prob_dde_2delays(u₀) + prob_dde_2delays -Model problem of finding a solution ``u(t)`` in the time span ``t \\in [0,1]`` to the delay -differential equation +Model problem of finding a solution ``u(t)`` in the time span ``t \\in [0,1]``, where +``t`` is of type `T`, to the delay differential equation ```math \\frac{du}{dt} = -u(t-1/3)-u(t-1/5) @@ -148,47 +170,54 @@ with history function ```math u(t) = \\begin{cases} -0 & \text{if } t < 0 \\ -u₀ & \text{if } t = 0 -\end{cases} +0 & \\text{if } t < 0 \\ +u₀ & \\text{if } t = 0 +\\end{cases} ``` -for ``t \\leq 0``. Hence the problem is discontinuous at ``t = 0`` for all ``u₀ \\neq 0``. +for ``t \\leq 0`` with ``u₀ = 1``. Note that the problem is discontinuous at ``t = 0`` for +all ``u₀ \\neq 0``. An analytical solution of this problem is provided for ``t \\in (-\\infty,1]``. """ -prob_dde_2delays(u₀) = DDEProblem(f_2delays, t->[0.0], [u₀], (0.0, 1.0), [1//3, 1//5]) +prob_dde_2delays = build_prob_dde_2delays(1.0) ### Not in-place function -f_2delays_notinplace = function (t,u,h) - - h(t-1/3) - h(t-1/5) +function f_2delays_notinplace(t::T, u, h) where {T} + (- h(t - T(1//3)) .- h(t - T(1//5))) ./ oneunit(t) end -(f::typeof(f_2delays_notinplace))(::Type{Val{:analytic}}, t, u0) = +f_2delays_notinplace(::Type{Val{:analytic}}, t, u0) = f_2delays(Val{:analytic}, t, u0) #### Vectorized history function +build_prob_dde_2delays_notinplace(u₀, ::T=u₀) where {T} = + DDEProblem(f_2delays_notinplace, t -> [zero(u₀)], [u₀], (zero(T), oneunit(T)), + [T(1//3), T(1//5)]) + """ - prob_dde_2delays_notinplace(u₀) + prob_dde_2delays_notinplace Same as [`prob_dde_2delays`](@ref), but purposefully implemented with a not in-place function. """ -prob_dde_2delays_notinplace(u₀) = - DDEProblem(f_2delays_notinplace, t->[0.0], [u₀], (0.0, 1.0), [1//3, 1//5]) +prob_dde_2delays_notinplace = build_prob_dde_2delays_notinplace(1.0) #### Scalar history function +build_prob_dde_2delays_scalar_notinplace(u₀, ::T=u₀) where {T} = + DDEProblem(f_2delays_notinplace, t -> zero(u₀), u₀, (zero(T), oneunit(T)), + [T(1//3), T(1//5)]) + """ - prob_dde_2delays_scalar_notinplace(u₀) + prob_dde_2delays_scalar_notinplace -Same as [`prob_dde_2delays_notinplace`](@ref), but purposefully implemented with a scalar -history function. +Same as [`prob_dde_2delays_notinplace`](@ref), but purposefully implemented with a +scalar history function. """ -prob_dde_2delays_scalar_notinplace(u₀) = - DDEProblem(f_2delays_notinplace, t->0.0, u₀, (0.0, 1.0), [1//3, 1//5]) +prob_dde_2delays_scalar_notinplace = build_prob_dde_2delays_scalar_notinplace(1.0) # DDE examples without analytical solution @@ -196,11 +225,16 @@ prob_dde_2delays_scalar_notinplace(u₀) = ### In-place function -f_1delay_long = function (t,u,h,du) - du[1] = - h(t-0.2)[1] + u[1] +function f_1delay_long(t::T, u, h, du) where {T} + du[1] = (- h(t - T(1//5))[1] + u[1]) / oneunit(t) end +build_prob_dde_1delay_long(u₀, ::T=u₀) where {T} = + DDEProblem(f_1delay_long, t -> [zero(u₀)], [u₀], (zero(T), T(100)), [T(1//5)]) + """ + prob_dde_1delay_long + Model problem of finding a solution ``u(t)`` in the time span ``t \\in [0,100]`` to the delay differential equation @@ -213,44 +247,58 @@ with history function ```math u(t) = \\begin{cases} 0 & \\text{if } t < 0,\\ -1 & \\text{if } t = 0, +u₀ & \\text{if } t = 0, \\end{cases} ``` -for ``t \\leq 0``. Hence the problem is discontinuous at ``t = 0``. +for ``t \\leq 0`` with ``u₀ = 1``. Note that the problem is discontinuous at ``t = 0`` for +all ``u₀ \\neq 0``. """ -prob_dde_1delay_long = DDEProblem(f_1delay_long, t->[0.0], [1.0], - (0.0, 100.0),[0.2]) +prob_dde_1delay_long = build_prob_dde_1delay_long(1.0) ### Not in-place function -f_1delay_long_notinplace = function (t,u,h) - - h(t-0.2) + u +function f_1delay_long_notinplace(t::T, u, h) where {T} + (- h(t - T(1//5)) .+ u ) ./ oneunit(t) end +build_prob_dde_1delay_long_notinplace(u₀, ::T=u₀) where {T} = + DDEProblem(f_1delay_long_notinplace, t -> [zero(u₀)], [u₀], (zero(T), T(100)), + [T(1//5)]) + """ -Same as [`prob_dde_1delay_long`](@ref), but purposefully implemented with a not in-place -function. + prob_dde_1delay_long_notinplace + +Same as [`prob_dde_1delay_long`](@ref), but purposefully implemented with a not +in-place function. """ -prob_dde_1delay_long_notinplace = - DDEProblem(f_1delay_long_notinplace, t->[0.0], [1.0], (0.0, 100.0), [0.2]) +prob_dde_1delay_long_notinplace = build_prob_dde_1delay_long_notinplace(1.0) + +build_prob_dde_1delay_long_scalar_notinplace(u₀, ::T=u₀) where {T} = + DDEProblem(f_1delay_long_notinplace, t -> zero(u₀), u₀, (zero(T), T(100)), [T(1//5)]) """ + prob_dde_1delay_long_scalar_notinplace + Same as [`prob_dde_1delay_long_notinplace`](@ref), but purposefully implemented with a scalar history function. """ -prob_dde_1delay_long_scalar_notinplace = - DDEProblem(f_1delay_long_notinplace, t->0.0, 1.0, (0.0, 100.0), [0.2]) +prob_dde_1delay_long_scalar_notinplace = build_prob_dde_1delay_long_scalar_notinplace(1.0) ## Two constant delays ### In-place function -f_2delays_long = function (t,u,h,du) - du[1] = - h(t-1/3)[1] - h(t-1/5)[1] +function f_2delays_long(t::T, u, h, du) where {T} + du[1] = (- h(t - T(1//3))[1] - h(t - T(1//5))[1]) / oneunit(t) end +build_prob_dde_2delays_long(u₀, ::T=u₀) where {T} = + DDEProblem(f_2delays_long, t -> [zero(u₀)], [u₀], (zero(T), T(100)), [T(1//3), T(1//5)]) + """ + prob_dde_2delays_long + Model problem of finding a solution ``u(t)`` in the time span ``t \\in [0,100]`` to the delay differential equation @@ -263,37 +311,320 @@ with history function ```math u(t) = \\begin{cases} 0 & \\text{if } t < 0,\\ -1 & \\text{if } t = 0, +u₀ & \\text{if } t = 0, \\end{cases} ``` -for ``t < 0``. Hence the problem is discontinuous at ``t = 0``. +for ``t < 0`` with ``u₀ = 1``. Note that the problem is discontinuous at ``t = 0`` for all +``u₀ \\neq 0``. """ -prob_dde_2delays_long = DDEProblem(f_2delays_long, t->[0.0], [1.0], - (0.0, 100.0), [1//3, 1//5]) +prob_dde_2delays_long = build_prob_dde_2delays_long(1.0) ### Not in-place function -f_2delays_long_notinplace = function (t,u,h) - - h(t-1/3) - h(t-1/5) +function f_2delays_long_notinplace(t::T, u, h) where {T} + (- h(t - T(1//3)) .- h(t - T(1//5))) ./ oneunit(t) end #### Vectorized history function +build_prob_dde_2delays_long_notinplace(u₀, ::T=u₀) where {T} = + DDEProblem(f_2delays_long_notinplace, t -> [zero(u₀)], [u₀], (zero(T), T(100)), + [T(1//3), T(1//5)]) + """ -Same as [`prob_dde_2delays_long`](@ref), but purposefully implemented with a not in-place -function. + prob_dde_2delays_long_notinplace + +Same as [`prob_dde_2delays_long`](@ref), but purposefully implemented with a not +in-place function. """ -prob_dde_2delays_long_notinplace = - DDEProblem(f_2delays_long_notinplace, t->[0.0], [1.0], - (0.0, 100.0), [1//3, 1//5]) +prob_dde_2delays_long_notinplace = build_prob_dde_2delays_long_notinplace(1.0) #### Scalar history function +build_prob_dde_2delays_long_scalar_notinplace(u₀, ::T=u₀) where {T} = + DDEProblem(f_2delays_long_notinplace, t -> zero(u₀), u₀, (zero(T), T(100)), + [T(1//3), T(1//5)]) + """ -Same as [`prob_dde_2delays_long_notinplace`](@ref), but purposefully implemented with a scalar -history function. + prob_dde_2delays_long_scalar_notinplace + +Same as [`prob_dde_2delays_long_notinplace`](@ref), but purposefully implemented with a +scalar history function. +""" +prob_dde_2delays_long_scalar_notinplace = build_prob_dde_2delays_long_scalar_notinplace(1.0) + +#= +The following examples are taken from: +W.H. Enright and H. Hayashi, The evaluation of numerical software for delay +differential equations. +=# + +function f_dde_mackey(t, u, h, du) + du[1] = 0.2*h(t-14)[1]/(1 + h(t-14)[1]^10) - 0.1*u[1] +end + +""" + prob_dde_mackey + +Model of blood production with constant delay (M. C. Mackey and L. Glass, Oscillation and +chaos in physiological control systems, 1977). +""" +prob_dde_mackey = DDEProblem(f_dde_mackey, t -> [0.5], [0.5], (0.0, 500.0), [14]) + +function f_dde_wheldon(t, u, h, du) + du[1] = 1.1/(1 + sqrt(10)*(h(t-20)[1])^(5/4)) - 10*u[1]/(1 + 40*u[2]) + du[2] = 100*u[1]/(1 + 40*u[2]) - 2.43*u[2] +end +u0_wheldon = [1.05767027/3; 1.030713491/3] + +""" + prob_dde_wheldon + +Model of chronic granulocytic leukemia with constant delay (T. Wheldon, J. Kirk and +H. Finlay, Cyclical granulopoiesis in chronic granulocytic leukemia: A simulation study, 1974). +""" +prob_dde_wheldon = DDEProblem(f_dde_wheldon, t -> u0_wheldon, u0_wheldon, (0., 100.), [20]) + +function f_dde_neves1(t, u, h, du) + du[1] = 1 - h(exp(1-1/t))[1] +end +# only valid for specific history function +function f_dde_neves1(::Type{Val{:analytic}}, t, u₀) + 0 < t ≤ 10 && u₀ == [log(0.1)] && return [log(t)] + error("This analytical solution is only valid on (0, 10] and for history function ϕ(t) = ln(t) for 0 < t ≤ 0.1") +end + +""" + prob_dde_neves1 + +DDE with vanishing time dependent delay at ``t = 1`` (K. W. Neves, Automatic integratorion +of functional differential equations: An approach, 1975). +""" +prob_dde_neves_1 = DDEProblem(f_dde_neves1, t -> [log(t)], [log(0.1)], (0.1, 10.), [], + [(t, u) -> t - exp(1 - 1/t)]) + +function f_dde_neves_thompson(t, u, h, du) + if h(t/2)[1] < 0 + du[1] = 1 - u[1] + else + du[1] = -1 - u[1] + end +end +# only valid for specific history function +function f_dde_neves_thompson(::Type{Val{:analytic}}, t, u₀) + if u₀ == [1] + if 0 ≤ t ≤ 2*log(2) + return [2*exp(-t) - 1] + elseif t ≤ 2*log(6) + return [1 - 6*exp(-t)] + elseif t ≤ 2*log(66) + return [66*exp(-t) - 1] + end + end + error("This analytical solution is only valid on [0, 2ln(66)] and for history function ϕ(t) = 1 for t ≤ 0") +end + +""" + prob_dde_neves_thompson + +DDE with vanishing time dependent delay at ``t = 0`` (K. W. Neves and S. Thompson, Solution +of systems of functional differential equations with state dependent delays, 1992), given by + +```math +u'(t) = \\begin{cases} +-1 - u(t) & \\text{if } u(t/2) \\geq 0 \\ +1 - u(t) & \\text{if } u(t/2) < 0 +\\end{case} +``` + +for ``t \\in [0, 2*\\log 66]`` with history function + +```math +u(t) = 1 +``` + +for ``t \\leq 0``. +""" +prob_dde_neves_thompson = DDEProblem(f_dde_neves_thompson, t -> [1.], [1.], (0., 2*log(66)), + [], [(t, u) -> t/2]) + +function f_dde_paul1(t, u, h, du) + du[1] = - 2*h(t - 1 - abs(u[1]))[1]*(1 - u[1]^2) +end + +""" + prob_dde_paul1 + +DDE with state dependent delay (C. A. H. Paul, A test set of functional differential +equations, 1994), given by + +```math +u'(t) = - 2u(t - 1 - |u(t)|)(1 - u(t)^2) +``` + +for ``t \\in [0, 30]`` with history function + +```math +u(t) = 1/2 +``` + +for ``t \\leq 0``. +""" +prob_dde_paul1 = DDEProblem(f_dde_paul1, t -> [0.5], [0.5], (0., 30.), [], + [(t, u) -> 1 + abs(u[1])]) + +function f_dde_paul2(t, u, h, du) + h1 = h(t - u[2])[1] + du[1] = -2*h1 + du[2] = (abs(h1) - abs(u[1]))/(1 + abs(h1)) +end + +""" + prob_dde_paul2 + +DDE with state dependent delay (C. A. H. Paul, A test set of functional differential equations, 1994). +""" +prob_dde_paul2 = DDEProblem(f_dde_paul2, t -> [1; 0.5], [1; 0.5], (0., 30.), [], + [(t, u) -> u[2]]) + +function build_prob_dde_mahaffy(tspan, h, σ₀, T₁, γ, Q, k, a, K, r) + function f(t, u, h, du) + du[1] = σ₀*h(t-T₁)[2] - γ*u[1] - Q + du[2] = a/(1 + K*u[1]^r) - k*u[2] + du[3] = 1 - Q*exp(γ*u[3])/(σ₀*h(t-T₁-u[3])[2]) + end + + DDEProblem(f, h, h(0), tspan, [T₁], [(t, u) -> T₁ + u[3]]) +end + +function h_mahaffy1(t) + if t < -6 + return [3.325; 9.5; 120] + else + return [3.325; 10; 120] + end +end +prob_dde_mahaffy1 = build_prob_dde_mahaffy((0., 300.), h_mahaffy1, 0.0031, 6, 0.001, 0.0275, + 2.8, 6570, 0.0382, 6.96) + +prob_dde_mahaffy2 = build_prob_dde_mahaffy((0., 100.), t -> [3.5; 10; 50], 0.00372, 3, 0.1, + 0.00178, 6.65, 15600, 0.0382, 6.96) + +""" +Two models of hematopoiesis with state dependent delay (J. M. Mahaffy, J. Bélair and +M. C. Mackey, Hematopoietic model with moving boundary condition and state dependent delay, 1996). +""" +prob_dde_mahaffy1, prob_dde_mahaffy2 + +function f_dde_neves2(t, u, h, du) + t2 = exp(1 - u[2]) + du[1] = u[2] + du[2] = -h(t2)[2]*u[2]^2*t2 +end +# only valid for specific history function +function f_dde_neves2(::Type{Val{:analytic}}, t, u₀) + 0 < t ≤ 5 && u₀ == [log(0.1); 10] && return [log(t); 1/t] + error("This analytical solution is only valid on (0, 5) and for history function ϕ(t) = [ln(t); 1/t] for 0 < t ≤ 0.1") +end + +""" + prob_dde_neves2 + +DDE with vanishing state dependent delay at ``t = 1`` (K. W. Neves, Automatic integration +of functional differential equations: An approach, 1975). +""" +prob_dde_neves2 = DDEProblem(f_dde_neves2, t -> [log(t); 1/t], [log(0.1); 10], (0.1, 5.), + [], [(t, u) -> t - exp(1 - u[2])]) + +function build_f_dde_gatica(r₁, r₂, α, δ) + function f_dde_gatica(t, u, h, du) + u₁u₂ = u[1]*u[2] + r₁u₁u₂ = r₁*u₁u₂ + r₂u₃ = r₂*u[3] + uₕ = h(t - u[4]) + uₕ₁₂ = uₕ[1]*uₕ[2] + + du[1] = -r₁u₁u₂ + r₂u₃ + du[2] = -r₁u₁u₂ + α*r₁*uₕ₁₂ + du[3] = r₁u₁u₂ - r₂u₃ + du[4] = 1 + (3*δ - u₁u₂ - u[3])/(uₕ₁₂ + uₕ[3])*exp(δ*u[4]) + end +end + +""" + prob_dde_gatica + +Model of antigen antibody dynamics with fading memory, with vanishing state dependent delay +at ``t = 0`` (J. Gatica and P. Waltman, A threshold model of antigen antibody dynamics with fading memory, 1982). +""" +prob_dde_gatica = DDEProblem(build_f_dde_gatica(0.02, 0.005, 3, 0.01), t -> [5; 0.1; 0; 0], + [5; 0.1; 0; 0], (0., 40.), [], [(t, u) -> u[4]]) + +#= +Quorum Sensing model +=# +function build_prob_dde_qs(u₀, tspan, τ, D, γₛ, Kₘ, nₛ, a, αₐ, βₐ, C₁, n₁, γₐ, αC, R, γ₃, Kₑ, αₗ, C₂, n₂, γₗ) + function hill_growth(x, K, n) + if x > K + return 1/(1 + (K / x)^n) + end + + if K == 0 + return one(x) + end + + y = (x / K)^n + return y / (1 + y) + end + + S₀ = u₀[1] + + function f_dde_qs(t, u, h, du) + if u[1] < 0 + # eq. 1 and 2 not defined for u[1] < 0 + du[1] = 0 + du[2] = 0 + else + tmp = u[2] * hill_growth(u[1], Kₘ, nₛ) + du[1] = D * (S₀ - u[1]) - γₛ * tmp + du[2] = a * tmp - D * u[2] + end + + if u[4] < 0 + # eq. 3 not defined for u[4] < 0 + du[3] = 0 + du[4] = max(0, αC * (R - u[4]) * u[3] - γ₃ * u[4]) + else + du[3] = (αₐ + βₐ * hill_growth(u[4], C₁, n₁)) * u[2] - + (γₐ + D) * u[3] - αC * (R - u[4]) * u[3] + γ₃ * u[4] - + Kₑ * u[3] * u[5] + du[4] = αC * (R - u[4]) * u[3] - γ₃ * u[4] + end + + Cₜ = h(t - τ)[4] + if Cₜ < 0 + # eq. 5 not defined for Cₜ < 0 + du[5] = 0 + else + du[5] = αₗ * hill_growth(Cₜ, C₂, n₂) * u[2] - (γₗ + D) * u[5] + end + end + + DDEProblem(f_dde_qs, t -> u₀, u₀, tspan, [τ]) +end + +""" + prob_dde_qs + +Model of Quorum Sensing (QS) of Pseudomonas putida IsoF in continuous cultures with constant +delay (Buddrus-Schiemann et al., Analysis of N-Acylhomoserine Lactone Dynamics in Continuous +Cultures of Pseudomonas Putida IsoF By Use of ELISA and UHPLC/qTOF-MS-derived Measurements +and Mathematical Models, Analytical and Bioanalytical Chemistry, 2014). """ -prob_dde_2delays_long_scalar_notinplace = - DDEProblem(f_2delays_long_notinplace, t->0.0, 1.0, - (0.0, 100.0), [1//3, 1//5]) +prob_dde_qs = build_prob_dde_qs([1; 8.4e8; 2.5e-9; 7.6e-8; 5e-15], # initial values + (0., 45.), # time span + 2, # delay + 0.1, 1.3e-12, 0.38, 1.3, 0.66, 2.3e-19, 2.3e-18, 70e-9, + 2.3, 0.05, 4e4, 5e-7, 0.080, 1.5e-4, 1.1e-8, 70e-9, 2.5, + 0.005)