Skip to content

Commit 84401ac

Browse files
update the problems to not use overloads
1 parent edb605f commit 84401ac

File tree

2 files changed

+82
-96
lines changed

2 files changed

+82
-96
lines changed

src/dde_premade_problems.jl

Lines changed: 24 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -25,7 +25,7 @@ function f_1delay(du, u, h, p, t)
2525
du[1] = - h(p, t - oneunit(t))[1] / oneunit(t)
2626
end
2727

28-
function f_1delay(::Type{Val{:analytic}}, u₀, p, t)
28+
function f_1delay_analytic(u₀, p, t)
2929
z = t/oneunit(t)
3030

3131
if z < 0
@@ -63,8 +63,10 @@ function f_1delay(::Type{Val{:analytic}}, u₀, p, t)
6363
end
6464
end
6565

66+
ff_1delay = DDEFunction(f_1delay,analytic=f_1delay_analytic)
67+
6668
build_prob_dde_1delay(u₀, ::T=u₀) where {T} =
67-
DDEProblem(f_1delay, [u₀], (p, t)->[zero(u₀)], (zero(T), T(10)),
69+
DDEProblem(ff_1delay, [u₀], (p, t)->[zero(u₀)], (zero(T), T(10)),
6870
constant_lags = [oneunit(T)])
6971

7072
"""
@@ -99,12 +101,12 @@ function f_1delay_notinplace(u, h, p, t)
99101
- h(p, t - oneunit(t)) ./ oneunit(t)
100102
end
101103

102-
f_1delay_notinplace(::Type{Val{:analytic}}, u0, p, t) = f_1delay(Val{:analytic}, u0, p, t)
104+
ff_1delay_notinplace = DDEFunction(f_1delay_notinplace,analytic=f_1delay_analytic)
103105

104106
#### Vectorized history function
105107

106108
build_prob_dde_1delay_notinplace(u₀, ::T=u₀) where {T} =
107-
DDEProblem(f_1delay_notinplace, [u₀], (p, t)->[zero(u₀)], (zero(T), T(10)),
109+
DDEProblem(ff_1delay_notinplace, [u₀], (p, t)->[zero(u₀)], (zero(T), T(10)),
108110
constant_lags = [oneunit(T)])
109111

110112
"""
@@ -118,7 +120,7 @@ prob_dde_1delay_notinplace = build_prob_dde_1delay_notinplace(1.0)
118120
#### Scalar history function
119121

120122
build_prob_dde_1delay_scalar_notinplace(u₀, ::T=u₀) where {T} =
121-
DDEProblem(f_1delay_notinplace, u₀, (p, t) -> zero(u₀), (zero(T), T(10)),
123+
DDEProblem(ff_1delay_notinplace, u₀, (p, t) -> zero(u₀), (zero(T), T(10)),
122124
constant_lags = [oneunit(T)])
123125

124126
"""
@@ -137,7 +139,7 @@ function f_2delays(du, u, h, p, t::T) where T
137139
du[1] = (- h(p, t - T(1//3))[1] - h(p, t - T(1//5))[1]) / oneunit(t)
138140
end
139141

140-
function f_2delays(::Type{Val{:analytic}}, u₀, p, t)
142+
function f_2delays_analytic(u₀, p, t)
141143
z = t/oneunit(t)
142144

143145
if z < 0
@@ -173,8 +175,10 @@ function f_2delays(::Type{Val{:analytic}}, u₀, p, t)
173175
end
174176
end
175177

178+
ff_2delays = DDEFunction(f_2delays,analytic=f_2delays_analytic)
179+
176180
build_prob_dde_2delays(u₀, ::T=u₀) where {T} =
177-
DDEProblem(f_2delays, [u₀], (p, t) -> [zero(u₀)], (zero(T), oneunit(T)),
181+
DDEProblem(ff_2delays, [u₀], (p, t) -> [zero(u₀)], (zero(T), oneunit(T)),
178182
constant_lags = [T(1//3), T(1//5)])
179183

180184
"""
@@ -209,13 +213,12 @@ function f_2delays_notinplace(u, h, p, t::T) where T
209213
(- h(p, t - T(1//3)) .- h(p, t - T(1//5))) ./ oneunit(t)
210214
end
211215

212-
f_2delays_notinplace(::Type{Val{:analytic}}, u0, p, t) =
213-
f_2delays(Val{:analytic}, u0, p, t)
216+
ff_2delays_notinplace = DDEFunction(f_2delays_notinplace,analytic=f_2delays_analytic)
214217

215218
#### Vectorized history function
216219

217220
build_prob_dde_2delays_notinplace(u₀, ::T=u₀) where {T} =
218-
DDEProblem(f_2delays_notinplace, [u₀], (p, t) -> [zero(u₀)], (zero(T), oneunit(T)),
221+
DDEProblem(ff_2delays_notinplace, [u₀], (p, t) -> [zero(u₀)], (zero(T), oneunit(T)),
219222
constant_lags = [T(1//3), T(1//5)])
220223

221224
"""
@@ -229,7 +232,7 @@ prob_dde_2delays_notinplace = build_prob_dde_2delays_notinplace(1.0)
229232
#### Scalar history function
230233

231234
build_prob_dde_2delays_scalar_notinplace(u₀, ::T=u₀) where {T} =
232-
DDEProblem(f_2delays_notinplace, u₀, (p, t) -> zero(u₀), (zero(T), oneunit(T)),
235+
DDEProblem(ff_2delays_notinplace, u₀, (p, t) -> zero(u₀), (zero(T), oneunit(T)),
233236
constant_lags = [T(1//3), T(1//5)])
234237

235238
"""
@@ -416,18 +419,20 @@ function f_dde_neves1(du, u, h, p, t)
416419
du[1] = 1 - h(p, exp(1-1/t))[1]
417420
end
418421
# only valid for specific history function
419-
function f_dde_neves1(::Type{Val{:analytic}}, u₀, p, t )
422+
function f_dde_neves1_analytic(u₀, p, t )
420423
0 < t 10 && u₀ == [log(0.1)] && return [log(t)]
421424
error("This analytical solution is only valid on (0, 10] and for history function ϕ(t) = ln(t) for 0 < t ≤ 0.1")
422425
end
423426

427+
ff_dde_neves1 = DDEFunction(f_dde_neves1,analytic=f_dde_neves1_analytic)
428+
424429
"""
425430
prob_dde_neves1
426431
427432
DDE with vanishing time dependent delay at ``t = 1`` (K. W. Neves, Automatic integratorion
428433
of functional differential equations: An approach, 1975).
429434
"""
430-
prob_dde_neves_1 = DDEProblem(f_dde_neves1, [log(0.1)], (p, t) -> [log(t)], (0.1, 10.),
435+
prob_dde_neves_1 = DDEProblem(ff_dde_neves1, [log(0.1)], (p, t) -> [log(t)], (0.1, 10.),
431436
dependent_lags = [(u,p,t) -> t - exp(1 - 1/t)])
432437

433438
function f_dde_neves_thompson(du, u, h, p, t)
@@ -438,7 +443,7 @@ function f_dde_neves_thompson(du, u, h, p, t)
438443
end
439444
end
440445
# only valid for specific history function
441-
function f_dde_neves_thompson(::Type{Val{:analytic}}, u₀, p, t )
446+
function f_dde_neves_thompson_analytic(u₀, p, t )
442447
if u₀ == [1]
443448
if 0 t 2*log(2)
444449
return [2*exp(-t) - 1]
@@ -451,6 +456,9 @@ function f_dde_neves_thompson(::Type{Val{:analytic}}, u₀, p, t )
451456
error("This analytical solution is only valid on [0, 2ln(66)] and for history function ϕ(t) = 1 for t ≤ 0")
452457
end
453458

459+
ff_dde_neves_thompson = DDEFunction(f_dde_neves_thompson,
460+
analytic=f_dde_neves_thompson_analytic)
461+
454462
"""
455463
prob_dde_neves_thompson
456464
@@ -472,7 +480,8 @@ u(t) = 1
472480
473481
for ``t \\leq 0``.
474482
"""
475-
prob_dde_neves_thompson = DDEProblem(f_dde_neves_thompson, [1.], (p, t) -> [1.], (0., 2*log(66)),
483+
prob_dde_neves_thompson = DDEProblem(ff_dde_neves_thompson, [1.],
484+
(p, t) -> [1.], (0., 2*log(66)),
476485
constant_lags = [(u,p,t) -> t/2])
477486

478487
function f_dde_paul1(du, u, h, p, t)

src/sde_premade_problems.jl

Lines changed: 58 additions & 81 deletions
Original file line numberDiff line numberDiff line change
@@ -10,9 +10,9 @@ export prob_sde_linear_stratonovich, prob_sde_2Dlinear_stratonovich
1010

1111
### SDE Examples
1212

13-
f = (u,p,t) -> 1.01u
14-
σ = (u,p,t) -> 0.87u
15-
(ff::typeof(f))(::Type{Val{:analytic}},u0,p,t,W) = u0.*exp.(0.63155t+0.87W)
13+
f_linear(u,p,t) = 1.01u
14+
σ_linear(u,p,t) = 0.87u
15+
linear_analytic(u0,p,t,W) = @.(u0*exp(0.63155t+0.87W))
1616

1717
@doc doc"""
1818
```math
@@ -25,24 +25,15 @@ u(u0,p,t,W_t)=u0\\exp((α-\\frac{β^2}{2})t+βW_t)
2525
```
2626
2727
"""
28-
prob_sde_linear = SDEProblem(f,σ,1/2,(0.0,1.0))
29-
30-
f = (u,p,t) -> 1.01u
31-
σ = (u,p,t) -> 0.87u
32-
(ff::typeof(f))(::Type{Val{:analytic}},u0,p,t,W) = u0.*exp.(1.01t+0.87W)
33-
prob_sde_linear_stratonovich = SDEProblem(f,σ,1/2,(0.0,1.0))
34-
35-
f = (du,u,p,t) -> begin
36-
for i = 1:length(u)
37-
du[i] = 1.01*u[i]
38-
end
39-
end
40-
σ = (du,u,p,t) -> begin
41-
for i in 1:length(u)
42-
du[i] = .87*u[i]
43-
end
44-
end
45-
(ff::typeof(f))(::Type{Val{:analytic}},u0,p,t,W) = u0.*exp.(0.63155*t+0.87*W)
28+
prob_sde_linear = SDEProblem(SDEFunction(f_linear,σ_linear,
29+
analytic=linear_analytic),σ_linear,1/2,(0.0,1.0))
30+
31+
linear_analytic_strat(u0,p,t,W) = @.(u0*exp(1.01t+0.87W))
32+
prob_sde_linear_stratonovich = SDEProblem(SDEFunction(f_linear,σ_linear,
33+
analytic=linear_analytic_strat),
34+
σ_linear,1/2,(0.0,1.0))
35+
f_linear_iip(du,u,p,t) = @.(du = 1.01*u)
36+
σ_linear_iip(du,u,p,t) = @.(du = 0.87*u)
4637
@doc doc"""
4738
8 linear SDEs (as a 4x2 matrix):
4839
@@ -55,24 +46,17 @@ where β=1.01, α=0.87, and initial condtion u0=1/2 with solution
5546
u(u0,p,t,W_t)=u0\\exp((α-\\frac{β^2}{2})t+βW_t)
5647
```
5748
"""
58-
prob_sde_2Dlinear = SDEProblem(f,σ,ones(4,2)/2,(0.0,1.0))
59-
60-
f = (du,u,p,t) -> begin
61-
for i = 1:length(u)
62-
du[i] = 1.01*u[i]
63-
end
64-
end
65-
σ = (du,u,p,t) -> begin
66-
for i in 1:length(u)
67-
du[i] = .87*u[i]
68-
end
69-
end
70-
(ff::typeof(f))(::Type{Val{:analytic}},u0,p,t,W) = u0.*exp.(1.01*t+0.87*W)
71-
prob_sde_2Dlinear_stratonovich = SDEProblem(f,σ,ones(4,2)/2,(0.0,1.0))
72-
73-
f = (u,p,t) -> -.25*u*(1-u^2)
74-
σ = (u,p,t) -> .5*(1-u^2)
75-
(ff::typeof(f))(::Type{Val{:analytic}},u0,p,t,W) = ((1+u0).*exp.(W)+u0-1)./((1+u0).*exp.(W)+1-u0)
49+
prob_sde_2Dlinear = SDEProblem(SDEFunction(f_linear_iip,σ_linear_iip,
50+
analytic=linear_analytic),
51+
σ_linear_iip,ones(4,2)/2,(0.0,1.0))
52+
prob_sde_2Dlinear_stratonovich = SDEProblem(SDEFunction(f_linear_iip,σ_linear_iip,
53+
analytic=linear_analytic_strat),
54+
σ_linear_iip,ones(4,2)/2,(0.0,1.0))
55+
56+
f_cubic(u,p,t) = -.25*u*(1-u^2)
57+
σ_cubic(u,p,t) = .5*(1-u^2)
58+
cubic_analytic(u0,p,t,W) = @. ((1+u0)*exp(W)+u0-1)/((1+u0)*exp(W)+1-u0)
59+
ff_cubic = SDEFunction(f_cubic,σ_cubic,analytic = cubic_analytic)
7660
@doc doc"""
7761
```math
7862
du_t = \\frac{1}{4}u(1-u^2)dt + \\frac{1}{2}(1-u^2)dW_t
@@ -84,11 +68,12 @@ and initial condtion u0=1/2, with solution
8468
u(u0,p,t,W_t)=\\frac{(1+u0)\\exp(W_t)+u0-1}{(1+u0)\\exp(W_t)+1-u0}
8569
```
8670
"""
87-
prob_sde_cubic = SDEProblem(f,σ,1/2,(0.0,1.0))
71+
prob_sde_cubic = SDEProblem(ff_cubic,σ_cubic,1/2,(0.0,1.0))
8872

89-
f = (u,p,t) -> -0.01*sin.(u).*cos.(u).^3
90-
σ = (u,p,t) -> 0.1*cos.(u).^2
91-
(ff::typeof(f))(::Type{Val{:analytic}},u0,p,t,W) = atan.(0.1*W + tan.(u0))
73+
f_wave(u,p,t) = @. -0.01*sin(u)*cos(u)^3
74+
σ_wave(u,p,t) = @. 0.1*cos(u)^2
75+
wave_analytic(u0,p,t,W) = @. atan(0.1*W + tan(u0))
76+
ff_wave = SDEFunction(f_wave,σ_wave,analytic=wave_analytic)
9277
@doc doc"""
9378
```math
9479
du_t = -\\frac{1}{100}\sin(u)\cos^3(u)dt + \\frac{1}{10}\cos^{2}(u_t) dW_t
@@ -100,13 +85,13 @@ and initial condition `u0=1.0` with solution
10085
u(u0,p,t,W_t)=\\arctan(\\frac{W_t}{10} + \\tan(u0))
10186
```
10287
"""
103-
prob_sde_wave = SDEProblem(f,σ,1.,(0.0,1.0))
88+
prob_sde_wave = SDEProblem(ff_wave,σ_wave,1.0,(0.0,1.0))
10489

105-
f = (u,p,t) -> p[2]./sqrt.(1+t) - u./(2*(1+t))
106-
σ = (u,p,t) -> p[1]*p[2]./sqrt.(1+t)
90+
f_additive(u,p,t) = @. p[2]/sqrt(1+t) - u/(2*(1+t))
91+
σ_additive(u,p,t) = @. p[1]*p[2]/sqrt(1+t)
10792
p = (0.1,0.05)
108-
(ff::typeof(f))(::Type{Val{:analytic}},u0,p,t,W) = u0./sqrt.(1+t) + p[2]*(t+p[1]*W)./sqrt.(1+t)
109-
93+
additive_analytic(u0,p,t,W) = @. u0/sqrt(1+t) + p[2]*(t+p[1]*W)/sqrt(1+t)
94+
ff_additive = SDEFunction(f_additive,σ_additive,analytic=additive_analytic)
11095
@doc doc"""
11196
Additive noise problem
11297
@@ -120,40 +105,24 @@ and initial condition u0=1.0 with α=0.1 and β=0.05, with solution
120105
u(u0,p,t,W_t)=\\frac{u0}{\\sqrt{1+t}} + \\frac{β(t+αW_t)}{\\sqrt{1+t}}
121106
```
122107
"""
123-
prob_sde_additive = SDEProblem(f,σ,1.,(0.0,1.0),p)
124-
125-
const sde_wave_αvec = [0.1;0.1;0.1;0.1]
126-
const sde_wave_βvec = [0.5;0.25;0.125;0.1115]
127-
f = (du,u,p,t) -> begin
128-
for i in 1:length(u)
129-
du[i] = sde_wave_βvec[i]/sqrt(1+t) - u[i]/(2*(1+t))
130-
end
131-
end
132-
133-
σ = (du,u,p,t) -> begin
134-
for i in 1:length(u)
135-
du[i] = sde_wave_αvec[i]*sde_wave_βvec[i]/sqrt(1+t)
136-
end
137-
end
138-
(ff::typeof(f))(::Type{Val{:analytic}},u0,p,t,W) = u0./sqrt(1+t) + sde_wave_βvec.*(t+sde_wave_αvec.*W)./sqrt(1+t)
108+
prob_sde_additive = SDEProblem(ff_additive,σ_additive,1.0,(0.0,1.0),p)
139109

110+
f_additive_iip(du,u,p,t) = @.(du = p[2]/sqrt(1+t) - u/(2*(1+t)))
111+
σ_additive_iip(du,u,p,t) = @.(du = p[1]*p[2]/sqrt(1+t))
112+
ff_additive_iip = SDEFunction(f_additive_iip,σ_additive_iip,analytic=additive_analytic)
113+
p = ([0.1;0.1;0.1;0.1],[0.5;0.25;0.125;0.1115])
140114
@doc doc"""
141115
A multiple dimension extension of `additiveSDEExample`
142116
143117
"""
144-
prob_sde_additivesystem = SDEProblem(f,σ,[1.;1.;1.;1.],(0.0,1.0))
118+
prob_sde_additivesystem = SDEProblem(ff_additive_iip,σ_additive_iip,[1.;1.;1.;1.],(0.0,1.0),p)
145119

146-
f = @ode_def_nohes LorenzSDE begin
120+
f_lorenz = @ode_def_bare LorenzSDE begin
147121
dx = σ*(y-x)
148122
dy = x*-z) - y
149123
dz = x*y - β*z
150124
end σ ρ β
151-
152-
σ = (du,u,p,t) -> begin
153-
for i in 1:3
154-
du[i] = 3.0 #Additive
155-
end
156-
end
125+
σ_lorenz(du,u,p,t) = @.(du = 3.0)
157126
@doc doc"""
158127
Lorenz Attractor with additive noise
159128
@@ -167,17 +136,18 @@ dz &= (x*y - β*z)dt + αdW_t \\\\
167136
168137
with ``σ=10``, ``ρ=28``, ``β=8/3``, ``α=3.0`` and inital condition ``u0=[1;1;1]``.
169138
"""
170-
prob_sde_lorenz = SDEProblem(f,σ,ones(3),(0.0,10.0),(10.0,28.0,2.66))
139+
prob_sde_lorenz = SDEProblem(f_lorenz,σ_lorenz,ones(3),(0.0,10.0),(10.0,28.0,2.66))
171140

172141

173-
f = (u,p,t) -> (1/3)*u^(1/3) + 6*u^(2/3)
174-
σ = (u,p,t) -> u^(2/3)
175-
(ff::typeof(f))(::Type{Val{:analytic}},u0,p,t,W) = (2t + 1 + W/3)^3
142+
f_nltest(u,p,t) = (1/3)*u^(1/3) + 6*u^(2/3)
143+
σ_nltest(u,p,t) = u^(2/3)
144+
analytic_nltest(u0,p,t,W) = (2t + 1 + W/3)^3
145+
ff_nltest = SDEFunction(f_nltest,σ_nltest,analytic=analytic_nltest)
176146
@doc doc"""
177147
Runge–Kutta methods for numerical solution of stochastic differential equations
178148
Tocino and Ardanuy
179149
"""
180-
prob_sde_nltest = SDEProblem(f,σ,1.0,(0.0,10.0))
150+
prob_sde_nltest = SDEProblem(ff_nltest,σ_nltest,1.0,(0.0,10.0))
181151

182152
function oval2ModelExample(;largeFluctuations=false,useBigs=false,noiseLevel=1)
183153
#Parameters
@@ -331,21 +301,28 @@ stiff_quad_f_ito(u,p,t) = -(p[1]+(p[2]^2)*u)*(1-u^2)
331301
stiff_quad_f_strat(u,p,t) = -p[1]*(1-u^2)
332302
stiff_quad_g(u,p,t) = p[2]*(1-u^2)
333303

334-
function stiff_quad_f_ito(::Type{Val{:analytic}},u0,p,t,W)
304+
function stiff_quad_f_ito_analytic(u0,p,t,W)
335305
α = p[1]
336306
β = p[2]
337307
exp_tmp = exp(-2*α*t+2*β*W)
338308
tmp = 1 + u0
339309
(tmp*exp_tmp + u0 - 1)/(tmp*exp_tmp - u0 + 1)
340310
end
341-
function stiff_quad_f_strat(::Type{Val{:analytic}},u0,p,t,W)
311+
312+
ff_stiff_quad_ito = SDEFunction(stiff_quad_f_ito,stiff_quad_g,
313+
analytic=stiff_quad_f_ito_analytic)
314+
315+
function stiff_quad_f_strat_analytic(u0,p,t,W)
342316
α = p[1]
343317
β = p[2]
344318
exp_tmp = exp(-2*α*t+2*β*W)
345319
tmp = 1 + u0
346320
(tmp*exp_tmp + u0 - 1)/(tmp*exp_tmp - u0 + 1)
347321
end
348322

323+
ff_stiff_quad_strat = SDEFunction(stiff_quad_f_strat,stiff_quad_g,
324+
analytic=stiff_quad_f_strat_analytic)
325+
349326
@doc doc"""
350327
The composite Euler method for stiff stochastic
351328
differential equations
@@ -363,7 +340,7 @@ Stiffness of Euler is determined by α+β²<1
363340
Higher α or β is stiff, with α being deterministic stiffness and
364341
β being noise stiffness (and grows by square).
365342
"""
366-
prob_sde_stiffquadito = SDEProblem(stiff_quad_f_ito,stiff_quad_g,0.5,(0.0,3.0),(1.0,1.0))
343+
prob_sde_stiffquadito = SDEProblem(ff_stiff_quad_ito,stiff_quad_g,0.5,(0.0,3.0),(1.0,1.0))
367344

368345
@doc doc"""
369346
The composite Euler method for stiff stochastic
@@ -382,7 +359,7 @@ Stiffness of Euler is determined by α+β²<1
382359
Higher α or β is stiff, with α being deterministic stiffness and
383360
β being noise stiffness (and grows by square).
384361
"""
385-
prob_sde_stiffquadstrat = SDEProblem(stiff_quad_f_strat,stiff_quad_g,0.5,(0.0,3.0),(1.0,1.0))
362+
prob_sde_stiffquadstrat = SDEProblem(ff_stiff_quad_strat,stiff_quad_g,0.5,(0.0,3.0),(1.0,1.0))
386363

387364
@doc doc"""
388365
Stochastic Heat Equation with scalar multiplicative noise

0 commit comments

Comments
 (0)