|
1 | 1 | using ModelingToolkit, JumpProcesses, LinearAlgebra, NonlinearSolve, Optimization, |
2 | 2 | OptimizationOptimJL, OrdinaryDiffEq, RecursiveArrayTools, SciMLBase, |
3 | 3 | SteadyStateDiffEq, StochasticDiffEq, DelayDiffEq, SymbolicIndexingInterface, |
4 | | - DiffEqCallbacks, Test, Plots |
| 4 | + DiffEqCallbacks, StochasticDelayDiffEq, Test, Plots |
5 | 5 | using ModelingToolkit: t_nounits as t, D_nounits as D |
6 | 6 |
|
7 | 7 | # Sets rnd number. |
|
952 | 952 | sol = solve(prob, MethodOfSteps(Tsit5())) |
953 | 953 | @test sol[sym] ≈ sol(sol.t .- sol.ps[delay]; idxs = original) |
954 | 954 | end |
| 955 | + |
| 956 | +@testset "SDDEs" begin |
| 957 | + function oscillator(; name, k = 1.0, τ = 0.01) |
| 958 | + @parameters k=k τ=τ |
| 959 | + @brownian a |
| 960 | + @variables x(..)=0.1 + t y(t)=0.1 + t jcn(t)=0.0 + t delx(t) |
| 961 | + eqs = [D(x(t)) ~ y + a, |
| 962 | + D(y) ~ -k * x(t - τ) + jcn, |
| 963 | + delx ~ x(t - τ)] |
| 964 | + return System(eqs, t; name = name) |
| 965 | + end |
| 966 | + systems = @named begin |
| 967 | + osc1 = oscillator(k = 1.0, τ = 0.01) |
| 968 | + osc2 = oscillator(k = 2.0, τ = 0.04) |
| 969 | + end |
| 970 | + eqs = [osc1.jcn ~ osc2.delx, |
| 971 | + osc2.jcn ~ osc1.delx] |
| 972 | + @named coupledOsc = System(eqs, t) |
| 973 | + @named coupledOsc = compose(coupledOsc, systems) |
| 974 | + sys = structural_simplify(coupledOsc) |
| 975 | + prob = SDDEProblem(sys, [], (0.0, 10.0); constant_lags = [sys.osc1.τ, sys.osc2.τ]) |
| 976 | + sym = sys.osc1.delx |
| 977 | + delay = sys.osc1.τ |
| 978 | + original = sys.osc1.x |
| 979 | + @test prob[sym] ≈ prob[original] .+ (prob.tspan[1] - prob.ps[delay]) |
| 980 | + integ = init(prob, ImplicitEM()) |
| 981 | + step!(integ, 10.0, true) |
| 982 | + @test integ[sym] ≈ SciMLBase.get_sol(integ)(integ.t - integ.ps[delay]; idxs = original) |
| 983 | + sol = solve(prob, ImplicitEM()) |
| 984 | + @test sol[sym] ≈ sol(sol.t .- sol.ps[delay]; idxs = original) |
| 985 | +end |
0 commit comments