Merge pull request #73 from SciML/sdde

stochastic delay differential equation analyticless convergence testing

Author: ChrisRackauckas

.travis.yml

@@ -4,7 +4,7 @@ os:
   - linux
 #  - osx
 julia:
-  - 1.3
+  - 1
 #  - nightly
 #matrix:
 #  allow_failures:

Project.toml

@@ -16,9 +16,9 @@ RootedTrees = "47965b36-3f3e-11e9-0dcf-4570dfd42a8c"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [compat]
-DelayDiffEq = "5"
-DiffEqBase = "6.11"
-DiffEqNoiseProcess = "3.3, 4.0, 5.0"
+DelayDiffEq = "5.20"
+DiffEqBase = "6.40"
+DiffEqNoiseProcess = "5.0"
 DiffEqProblemLibrary = "4.5"
 NLsolve = "4.2"
 RecipesBase = "0.7, 0.8, 1.0"
@@ -33,7 +33,8 @@ OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 ParameterizedFunctions = "65888b18-ceab-5e60-b2b9-181511a3b968"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"
+StochasticDelayDiffEq = "29a0d76e-afc8-11e9-03a4-eda52ae4b960"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["DelayDiffEq", "DiffEqProblemLibrary", "OrdinaryDiffEq", "ParameterizedFunctions", "Random", "StochasticDiffEq", "Test"]
+test = ["DelayDiffEq", "DiffEqProblemLibrary", "OrdinaryDiffEq", "ParameterizedFunctions", "Random", "StochasticDiffEq", "StochasticDelayDiffEq", "Test"]

src/DiffEqDevTools.jl

@@ -10,6 +10,7 @@ import Base: length
 
 import DiffEqBase: AbstractODEProblem, AbstractDDEProblem,
        AbstractODESolution, AbstractRODEProblem, AbstractSDEProblem,
+       AbstractSDDEProblem,
        AbstractDAEProblem, @def, ConvergenceSetup, DEAlgorithm,
        ODERKTableau, AbstractTimeseriesSolution, ExplicitRKTableau,
        ImplicitRKTableau

src/convergence.jl

@@ -57,7 +57,7 @@ function test_convergence(dts::AbstractArray,prob::Union{AbstractRODEProblem,Abs
 end
 
 function analyticless_test_convergence(dts::AbstractArray,
-                          prob::Union{AbstractRODEProblem,AbstractSDEProblem},
+                          prob::Union{AbstractRODEProblem,AbstractSDEProblem,AbstractSDDEProblem},
                           alg,test_dt;trajectories=100,
                           save_everystep=true,timeseries_steps=1,
                           timeseries_errors=save_everystep,adaptive=false,
@@ -76,9 +76,18 @@ function analyticless_test_convergence(dts::AbstractArray,
         brownian_values2 = cumsum([[zeros(size(prob.noise_rate_prototype,2))];[sqrt(test_dt)*randn(size(prob.noise_rate_prototype,2)) for i in 1:length(t)-1]])
       end
       np = NoiseGrid(t,brownian_values,brownian_values2)
-      _prob = SDEProblem(prob.f,prob.g,prob.u0,prob.tspan,
+
+      if prob isa AbstractSDDEProblem
+        _prob = SDDEProblem(prob.f,prob.g,prob.u0,prob.h,prob.tspan,prob.p,
+                            noise = np, noise_rate_prototype = prob.noise_rate_prototype,
+                            constant_lags = prob.constant_lags, dependent_lags = prob.dependent_lags,
+                            neutral = prob.neutral, order_discontinuity_t0 = prob.order_discontinuity_t0, prob.kwargs...);
+      else
+        _prob = SDEProblem(prob.f,prob.g,prob.u0,prob.tspan,prob.p,
                          noise=np,
                          noise_rate_prototype=prob.noise_rate_prototype);
+      end
+
       true_sol = solve(_prob,alg;adaptive=adaptive,dt=test_dt);
 
       for i in 1:length(dts)

test/analyticless_convergence_tests.jl

@@ -1,4 +1,4 @@
-using OrdinaryDiffEq, ParameterizedFunctions, Test, Random
+using OrdinaryDiffEq, StochasticDelayDiffEq, ParameterizedFunctions, Test, Random
 using ParameterizedFunctions.ModelingToolkit # macro hygiene
 f = @ode_def LotkaVolterra begin
   dx = 1.5x - x*y
@@ -39,3 +39,25 @@ test_dt = 1/2^8
 sim2 = analyticless_test_convergence(dts,prob,SRIW1(),test_dt,trajectories=100, use_noise_grid=false)
 @test abs(sim2.𝒪est[:final]-1.5) < 0.3
 @show sim2.𝒪est[:final]
+
+### SDDE
+
+function hayes_modelf(du,u,h,p,t)
+    τ,a,b,c,α,β,γ = p
+    du .= a.*u .+ b .* h(p,t-τ) .+ c
+end
+function hayes_modelg(du,u,h,p,t)
+    τ,a,b,c,α,β,γ = p
+    du .= α.*u .+ β.*h(p,t-τ) .+ γ
+end
+h(p,t) = (ones(1) .+ t);
+tspan = (0.,10.)
+
+pmul = [1.0,-4.,-2.,10.,-1.3,-1.2, 1.1]
+padd = [1.0,-4.,-2.,10.,-0.0,-0.0, 0.1]
+
+prob = SDDEProblem(hayes_modelf, hayes_modelg, [1.], h, tspan, pmul; constant_lags = (pmul[1],));
+dts = (1/2).^(7:-1:3)
+test_dt = 1/2^8
+sim2 = analyticless_test_convergence(dts,prob,RKMil(),test_dt,trajectories=100, use_noise_grid=false)
+@test abs(sim2.𝒪est[:final]-1.0) < 0.3