Is there a ddesd() function?

[frostRed]

I need to solve a time vary delay differential equations like this and I find dde() can only solve the const delay differential equation. The following code can not run.

using Plots
pyplot()

using DifferentialEquations

tau1(t) = 0.2 * sin(t) ^ 2
# tau1 = 1
function rossler(t, u, h, du)
    du[1] = -u[2] - u[3] + 0.3 * h(t - tau1(t))[1]
    du[2] = u[1] + 0.2* u[2]
    du[3] = 0.2 - 5.7 * u[3] + u[1] * u[2]
end

u0 = [1.2;1.2;1.2]
lags1 = [tau1]
h1(t) = zeros(3) 
tspan = (0.0,100.0)
prob = DDEProblem(rossler, h1, u0, tspan, lags1)
alg = MethodOfSteps(DP8())
sol = solve(prob,alg)

plt = plot(sol,vars=(1,2,3), xlabel = "x", ylabel = "y", zlabel = "z")

[ChrisRackauckas]

Hey, I think our tutorial needs some updates, but you can solve it. If you see

http://docs.juliadiffeq.org/latest/types/dde_types.html

you'll see there's a spot for telling the solver what the state-dependent lags (as a tuple of functions lag(t,u)) are. Here's an example from the tests, where it sets a "state-dependent delay" as a function which returns a constant (of course, it doesn't need to be a constant, but this is just for testing...)

https://github.com/JuliaDiffEq/DelayDiffEq.jl/blob/master/test/dependent_delays.jl#L14

So that's how you set it. I'll get a tutorial written up: I thought we already did that but I guess not.

Other information:

If you don't set the state-dependent lags, then using RK4() is the same as MATLAB's ddesd. This is actually quickly mentioned at the bottom of the tutorial:

http://docs.juliadiffeq.org/latest/tutorials/dde_example.html#Undeclared-Delays-1

While this works, ddesd is notoriously bad (even Shampine's paper shows this) because it isn't able to exactly handle the discontinuities, which is mentioned here:

http://docs.juliadiffeq.org/latest/solvers/dde_solve.html#Lag-Handling-1

If the tolerance is low enough, that's okay, but I'd recommend setting the state-dependent delays if you want the answer to be more accurate more efficiently so it can do discontinuity tracking.

Long story short, let me get that tutorial up to speed.

[ChrisRackauckas]

The docs are updated. Here's how you'd do this problem:

using Plots
using DifferentialEquations
plotly()

tau1(t,u) = 0.2 * sin(t) ^ 2
# tau1 = 1
function rossler(t, u, h, du)
    du[1] = -u[2] - u[3] + 0.3 * h(t - tau1(t,u))[1]
    du[2] = u[1] + 0.2* u[2]
    du[3] = 0.2 - 5.7 * u[3] + u[1] * u[2]
end

u0 = [1.2;1.2;1.2]
d_lags = [tau1]
h1(t) = zeros(3)
tspan = (0.0,100.0)

# ddesd version
prob = DDEProblem(rossler, h1, u0, tspan)
alg = MethodOfSteps(RK4())
sol = solve(prob,alg,abstol=1e-8,reltol=1e-8)
plot(sol,tspan=(0.0,11.5))

# delay tracking version
prob = DDEProblem(rossler, h1, u0, tspan, nothing, d_lags)
alg = MethodOfSteps(Tsit5())
sol = solve(prob,alg)

# Random stiff solver with delay tracking
alg = MethodOfSteps(SDIRK2())
sol = solve(prob,alg)

However, note that all 3 algorithms are blowing up in the exact same spot in the exact same way. This indicates there may be something up with the problem specification. Out of curiosity, is this an example equation from somewhere? MethodOfSteps(RK4()) should be the same as Shampine's ddesd, so if you have a MATLAB script to compare it to that would be useful just to make sure it's the right equation. Otherwise, if there is something up, @devmotion have you seen a Rossler equation in the DDE tests?

[frostRed]

I have read about the content of the dependent lags. The example of tutorials and above three version codes can't run due to the same error: MethodError: no method matching start(::Void)

Here is dependent lags example in tutorials:

using DifferentialEquations
const p0 = 0.2; const q0 = 0.3; const v0 = 1; const d0 = 5
const p1 = 0.2; const q1 = 0.3; const v1 = 1; const d1 = 1
const d2 = 1; const beta0 = 1; const beta1 = 1; 

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

# lags = [tau]
dependent_lags = (tau,)

h(t) = ones(3)
tspan = (0.0,10.0)
u0 = [1.0,1.0,1.0]

prob = DDEProblem(bc_model,h,u0,tspan,nothing,dependent_lags)
alg = MethodOfSteps(Tsit5())
sol = solve(prob,alg)

Here is my problem's code:

---

using Plots
using DifferentialEquations
plotly()

tau1(t,u) = 1
# tau1(t,u) = 0.2 * sin(t) ^ 2
function rossler(t, u, h, du)
    du[1] = -u[2] - u[3] + 0.3 * h(t - tau1(t,u))[1]
    du[2] = u[1] + 0.2* u[2]
    du[3] = 0.2 - 5.7 * u[3] + u[1] * u[3]
end

u0 = [1.2;1.2;1.2]
d_lags = (tau1, )
h1(t) = zeros(3)
tspan = (0.0,100.0)

# delay tracking version
prob = DDEProblem(rossler, h1, u0, tspan, nothing, d_lags)
alg = MethodOfSteps(Tsit5())
sol = solve(prob,alg)

---

but the example of test work well

using Plots
using DifferentialEquations
plotly()

tau(t, u) = 0.2 * sin(t) ^ 2

function f_1delay(t, u, h, du)
    du[1] = - h(t - tau(t, u))[1] / oneunit(t)
end
prob = DDEProblem(f_1delay, t -> [0.0], [1.0], (0., 10.), [],
                   [tau])
alg = MethodOfSteps(BS3())
sol = solve(prob, alg)
plot(sol)

[ChrisRackauckas]

You need to Pkg.update()

[devmotion]

Before the latest release of DelayDiffEq it was not possible to use nothing as constant delays (hence, as in our test example, one had to use an empty vector []). However, there is still a check missing in the current version which causes the issue you see. So until this problem is fixed, the easiest solution is to just use [] instead of nothing, i.e. try

prob = DDEProblem(rossler, h1, u0, tspan, [], d_lags)

instead of

prob = DDEProblem(rossler, h1, u0, tspan, nothing, d_lags)

[ChrisRackauckas]

@devmotion I pushed a change last night so nothing should work, and ran those examples? But yes, it's a small syntactic difference and [] works in place of nothing.

[devmotion]

Strange, I updated to the newest version and got the same error because of https://github.com/JuliaDiffEq/DelayDiffEq.jl/blob/master/src/integrator_utils.jl#L359. But a simple test if !(typeof(constant_lags) <: Void) could fix it.

[ChrisRackauckas]

huh, wonder why it was working for me without that. I added a constant_lags != nothing (note that lowers to the same thing at compile-time since it's a singleton type).

[ChrisRackauckas]

Alright, the tag is going through now for the constant_lags is nothing fix.

https://github.com/JuliaLang/METADATA.jl/pull/11840

In the meantime, is that equation with that delay supposed to give a stable equation? I don't think it does, because I can solve it with high accuracy and it always hits a snag at the same point:

using Plots
using DelayDiffEq
plotly()

tau1(t,u) = 1
# tau1(t,u) = 0.2 * sin(t) ^ 2
function rossler(t, u, h, du)
    du[1] = -u[2] - u[3] + 0.3 * h(t - tau1(t,u))[1]
    du[2] = u[1] + 0.2* u[2]
    du[3] = 0.2 - 5.7 * u[3] + u[1] * u[2]
end

u0 = [1.2;1.2;1.2]
d_lags = (tau1, )
h1(t) = zeros(3)
tspan = (0.0,100.0)

# delay tracking version
prob = DDEProblem(rossler, h1, u0, tspan, nothing, d_lags)
alg = MethodOfSteps(Tsit5())
sol = solve(prob,alg,abstol=1e-8,reltol=1e-8)

plot(sol,plotdensity=10000,tspan=(0.0,30.0))

Where is the test equation coming from?

[frostRed]

I made a mistake in entering the equation. It should be du[3] = 0.2 - 5.7 * u[3] + u[1] * u[3].

By the way, about the equation you can see rossler attractor. This equation is time-varying delay rossler attractor.

Thank you very much.

[ChrisRackauckas]

Cool thanks for the report and the example. Glad to help. Let us know if you have any other problems!

[PeterXiaoGuo]

@ChrisRackauckas Hi Chris, I face a problem using DiffernetialEquaitons.jl's dde solver and Shampine's ddesd solver in matlab cannot give me enough accuracy.

Could you please give an example how to use time-dependent delay, e.g. delay = sin(t)?

I guess some other algorithm like Rodas4() or Rosenbrock23() may helps. At least, I tried ddesd and it cannot provide expected result.

And after I directly run your above code in Julia, it throws an error:

MethodError: no method matching DDEProblem(::typeof(rossler), ::typeof(h1), ::Array{Float64,1}, ::Tuple{Float64,Float64}, ::Nothing, ::Tuple{typeof(tau1)})

I can run DiffernetialEquaitons.jl's documentation dde sample code but why still face error here?

Thank you very much!

[devmotion]

The problem is that there were some breaking changes in DelayDiffEq and so Chris' code above is completely outdated: the order of the arguments in the DDE function and the history function is incorrect and delays have to be specified as keywords. The following code works on my computer:

using Plots
using DelayDiffEq
plotly()

tau1(u, p, t) = 1 # to get delay = sin(t) you have to replace it by tau1(u, p, t) = sin(t)
function rossler(du, u, h, p, t)
    du[1] = -u[2] - u[3] + 0.3 * h(p, t - tau1(u, p, t))[1]
    du[2] = u[1] + 0.2* u[2]
    du[3] = 0.2 - 5.7 * u[3] + u[1] * u[2]
end

u0 = [1.2;1.2;1.2]
d_lags = (tau1, )
h1(p, t) = zeros(3)
tspan = (0.0,100.0)

# delay tracking version
prob = DDEProblem(rossler, u0, h1, tspan, dependent_lags = d_lags)
alg = MethodOfSteps(Tsit5())
sol = solve(prob,alg,abstol=1e-8,reltol=1e-8)

plot(sol,plotdensity=10000,tspan=(0.0,30.0))

You can find more details about solving DDE problems in the documentation (see e.g. the DDE tutorial and in particular the section about state-dependent delays, the documentation about DDE problems, and the section about DDE solvers). I hope that helps!

[PeterXiaoGuo]

@devmotion Hi David, thanks for your reply!! It works now!

In case anyone facing the same problem: Replace tau in http://docs.juliadiffeq.org/latest/tutorials/dde_example.html#State-Dependent-Delay-Discontinuity-Tracking-1 example with tau(u, p, t) when define differential model, otherwise, an error will appear.