Skip to content

Commit 02831fc

Browse files
ChrisRackauckasChrisRackauckas
authored
Merge pull request #1532 from ChrisRackauckas-Claude/stiff-dde-benchmarks
Add stiff DDE benchmark suite from RADAR5 test problems
2 parents a7f23db + e7b6eab commit 02831fc

File tree

4 files changed

+685
-0
lines changed

4 files changed

+685
-0
lines changed

‎benchmarks/StiffDDE/Enzyme_wpd.jmd‎

Lines changed: 172 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,172 @@
1+
---
2+
title: RADAR5 Enzyme Kinetics Work-Precision Diagrams
3+
author: Chris Rackauckas
4+
---
5+
6+
# RADAR5 Enzyme Kinetics
7+
8+
This is a stiff delay differential equation model of enzyme kinetics with an inhibitor
9+
molecule, taken from the RADAR5 test suite by Guglielmi and Hairer. The system has
10+
dimension 4 with one constant delay.
11+
12+
The model describes the dynamics of a four-component system where the inhibitor enters
13+
through a delayed feedback on the fourth component:
14+
15+
```math
16+
U(t) = \frac{1}{1 + K_1 \, u_4(t - \tau)^3}
17+
```
18+
```math
19+
u_1'(t) = 10.5 - u_1(t) \, U(t)
20+
```
21+
```math
22+
u_2'(t) = u_1(t) \, U(t) - u_2(t)
23+
```
24+
```math
25+
u_3'(t) = u_2(t) - u_3(t)
26+
```
27+
```math
28+
u_4'(t) = u_3(t) - 0.5 \, u_4(t)
29+
```
30+
31+
for $t \in [0, 160]$ with constant history $\phi(t) = (60, 10, 10, 20)$ for $t \leq 0$,
32+
where $\tau = 4$ and $K_1 = 0.0005$.
33+
34+
## References
35+
36+
Guglielmi, N. and Hairer, E. (2005). Users' Guide for the code RADAR5 - Version 2.1.
37+
38+
```julia
39+
using DelayDiffEq, DiffEqDevTools, Plots
40+
gr()
41+
```
42+
43+
## Problem Definition
44+
45+
```julia
46+
function f_enzyme!(du, u, h, p, t)
47+
τ, K1 = p
48+
y4_delayed = h(p, t - τ; idxs = 4)
49+
U = 1.0 / (1.0 + K1 * y4_delayed^3)
50+
du[1] = 10.5 - u[1] * U
51+
du[2] = u[1] * U - u[2]
52+
du[3] = u[2] - u[3]
53+
du[4] = u[3] - 0.5 * u[4]
54+
nothing
55+
end
56+
57+
function h_enzyme(p, t; idxs::Union{Nothing, Int} = nothing)
58+
if idxs === nothing
59+
[60.0, 10.0, 10.0, 20.0]
60+
elseif idxs == 1
61+
60.0
62+
elseif idxs == 2
63+
10.0
64+
elseif idxs == 3
65+
10.0
66+
elseif idxs == 4
67+
20.0
68+
else
69+
error("index must be between 1 and 4")
70+
end
71+
end
72+
73+
τ = 4.0
74+
K1 = 0.0005
75+
prob = DDEProblem(f_enzyme!, h_enzyme, (0.0, 160.0), (τ, K1); constant_lags = [τ])
76+
```
77+
78+
## Reference Solution
79+
80+
```julia
81+
sol = solve(prob, MethodOfSteps(Rodas5P()); reltol = 1e-14, abstol = 1e-14)
82+
test_sol = TestSolution(sol)
83+
plot(sol; title = "Enzyme Kinetics Solution")
84+
```
85+
86+
## High Tolerances
87+
88+
### Rosenbrock methods
89+
90+
```julia
91+
abstols = 1.0 ./ 10.0 .^ (4:7)
92+
reltols = 1.0 ./ 10.0 .^ (1:4)
93+
94+
setups = [Dict(:alg => MethodOfSteps(Rosenbrock23())),
95+
Dict(:alg => MethodOfSteps(Rodas4())),
96+
Dict(:alg => MethodOfSteps(Rodas5())),
97+
Dict(:alg => MethodOfSteps(Rodas5P()))]
98+
names = ["Rosenbrock23", "Rodas4", "Rodas5", "Rodas5P"]
99+
wp = WorkPrecisionSet(prob, abstols, reltols, setups;
100+
names = names, appxsol = test_sol, maxiters = Int(1e5), error_estimate = :final)
101+
plot(wp; title = "Enzyme: Rosenbrock Methods (final error)")
102+
```
103+
104+
```julia
105+
wp = WorkPrecisionSet(prob, abstols, reltols, setups;
106+
names = names, appxsol = test_sol, maxiters = Int(1e5), error_estimate = :L2)
107+
plot(wp; title = "Enzyme: Rosenbrock Methods (L2 error)")
108+
```
109+
110+
### SDIRK methods
111+
112+
```julia
113+
setups = [Dict(:alg => MethodOfSteps(TRBDF2())),
114+
Dict(:alg => MethodOfSteps(SDIRK2())),
115+
Dict(:alg => MethodOfSteps(KenCarp4())),
116+
Dict(:alg => MethodOfSteps(Kvaerno4())),
117+
Dict(:alg => MethodOfSteps(Kvaerno5()))]
118+
names = ["TRBDF2", "SDIRK2", "KenCarp4", "Kvaerno4", "Kvaerno5"]
119+
wp = WorkPrecisionSet(prob, abstols, reltols, setups;
120+
names = names, appxsol = test_sol, maxiters = Int(1e5), error_estimate = :final)
121+
plot(wp; title = "Enzyme: SDIRK Methods (final error)")
122+
```
123+
124+
```julia
125+
wp = WorkPrecisionSet(prob, abstols, reltols, setups;
126+
names = names, appxsol = test_sol, maxiters = Int(1e5), error_estimate = :L2)
127+
plot(wp; title = "Enzyme: SDIRK Methods (L2 error)")
128+
```
129+
130+
### Stiff vs Non-Stiff Comparison
131+
132+
```julia
133+
setups = [Dict(:alg => MethodOfSteps(Rodas5P())),
134+
Dict(:alg => MethodOfSteps(TRBDF2())),
135+
Dict(:alg => MethodOfSteps(KenCarp4())),
136+
Dict(:alg => MethodOfSteps(Tsit5())),
137+
Dict(:alg => MethodOfSteps(DP5())),
138+
Dict(:alg => MethodOfSteps(OwrenZen5()))]
139+
names = ["Rodas5P", "TRBDF2", "KenCarp4", "Tsit5", "DP5", "OwrenZen5"]
140+
wp = WorkPrecisionSet(prob, abstols, reltols, setups;
141+
names = names, appxsol = test_sol, maxiters = Int(1e6), error_estimate = :final)
142+
plot(wp; title = "Enzyme: Stiff vs Non-Stiff (final error)")
143+
```
144+
145+
## Low Tolerances
146+
147+
```julia
148+
abstols = 1.0 ./ 10.0 .^ (8:11)
149+
reltols = 1.0 ./ 10.0 .^ (5:8)
150+
151+
setups = [Dict(:alg => MethodOfSteps(Rosenbrock23())),
152+
Dict(:alg => MethodOfSteps(Rodas4())),
153+
Dict(:alg => MethodOfSteps(Rodas5())),
154+
Dict(:alg => MethodOfSteps(Rodas5P())),
155+
Dict(:alg => MethodOfSteps(TRBDF2())),
156+
Dict(:alg => MethodOfSteps(KenCarp4()))]
157+
names = ["Rosenbrock23", "Rodas4", "Rodas5", "Rodas5P", "TRBDF2", "KenCarp4"]
158+
wp = WorkPrecisionSet(prob, abstols, reltols, setups;
159+
names = names, appxsol = test_sol, maxiters = Int(1e5), error_estimate = :final)
160+
plot(wp; title = "Enzyme: Low Tolerances (final error)")
161+
```
162+
163+
```julia
164+
wp = WorkPrecisionSet(prob, abstols, reltols, setups;
165+
names = names, appxsol = test_sol, maxiters = Int(1e5), error_estimate = :L2)
166+
plot(wp; title = "Enzyme: Low Tolerances (L2 error)")
167+
```
168+
169+
```julia, echo = false
170+
using SciMLBenchmarks
171+
SciMLBenchmarks.bench_footer(WEAVE_ARGS[:folder], WEAVE_ARGS[:file])
172+
```

‎benchmarks/StiffDDE/Hepatitis_wpd.jmd‎

Lines changed: 208 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,208 @@
1+
---
2+
title: RADAR5 Hepatitis B Work-Precision Diagrams
3+
author: Chris Rackauckas
4+
---
5+
6+
# RADAR5 Hepatitis B Virus Infection
7+
8+
This is a stiff delay differential equation model of acute hepatitis B virus infection,
9+
taken from the RADAR5 test suite by Guglielmi and Hairer. The system has dimension 10
10+
with 5 constant delays ($\tau_1 = 0.6$, $\tau_2 = 0.6$, $\tau_3 = 2.0$, $\tau_4 = 2.0$,
11+
$\tau_5 = 3.0$ days).
12+
13+
The model describes the dynamics of the immune response to hepatitis B virus infection,
14+
including free virus ($u_1$), infected hepatocytes ($u_2$), antibody-bound virus ($u_3$),
15+
CTL effector cells ($u_4$), and various immune cell populations ($u_5$--$u_{10}$). The
16+
delays represent the time lags in the immune response cascade.
17+
18+
The problem is extremely stiff due to rate constants spanning over 50 orders of magnitude
19+
(from $O(10^{-22})$ to $O(10^{33})$), making it an excellent test case for stiff DDE
20+
solvers.
21+
22+
## References
23+
24+
Guglielmi, N. and Hairer, E. (2005). Users' Guide for the code RADAR5 - Version 2.1.
25+
26+
```julia
27+
using DelayDiffEq, DiffEqDevTools, Plots
28+
gr()
29+
```
30+
31+
## Problem Definition
32+
33+
```julia
34+
function f_hepatitis!(du, u, h, p, t)
35+
α = p
36+
xi(z) = 1.0 - z / α[7]
37+
38+
# Delayed values for 5 delays
39+
y4_τ1 = h(p, t - α[41]; idxs = 4)
40+
y4_τ2 = h(p, t - α[42]; idxs = 4)
41+
y4_τ3 = h(p, t - α[43]; idxs = 4)
42+
y4_τ4 = h(p, t - α[44]; idxs = 4)
43+
y4_τ5 = h(p, t - α[45]; idxs = 4)
44+
45+
y5_τ1 = h(p, t - α[41]; idxs = 5)
46+
y5_τ3 = h(p, t - α[43]; idxs = 5)
47+
48+
y6_τ2 = h(p, t - α[42]; idxs = 6)
49+
y6_τ4 = h(p, t - α[44]; idxs = 6)
50+
y6_τ5 = h(p, t - α[45]; idxs = 6)
51+
52+
y7_τ3 = h(p, t - α[43]; idxs = 7)
53+
54+
y8_τ4 = h(p, t - α[44]; idxs = 8)
55+
y8_τ5 = h(p, t - α[45]; idxs = 8)
56+
57+
du[1] = α[1] * u[2] + α[2] * α[3] * u[2] * u[7] - α[4] * u[1] * u[10] -
58+
α[5] * u[1] - α[6] * u[1] * (α[7] - u[2] - u[3])
59+
du[2] = α[8] * u[1] * (α[7] - u[2] - u[3]) - α[3] * u[2] * u[7] - α[9] * u[2]
60+
du[3] = α[3] * u[2] * u[7] + α[9] * u[2] - α[10] * u[3]
61+
du[4] = α[11] * α[12] * u[1] - α[13] * u[4]
62+
du[5] = α[14] * (xi(u[3]) * α[15] * y4_τ1 * y5_τ1 - u[4] * u[5]) -
63+
α[16] * u[4] * u[5] * u[7] + α[17] * (α[18] - u[5])
64+
du[6] = α[19] * (xi(u[3]) * α[20] * y4_τ2 * y6_τ2 - u[4] * u[6]) -
65+
α[21] * u[4] * u[6] * u[8] + α[22] * (α[23] - u[6])
66+
du[7] = α[24] * (xi(u[3]) * α[25] * y4_τ3 * y5_τ3 * y7_τ3 -
67+
u[4] * u[5] * u[7]) - α[26] * u[2] * u[7] + α[27] * (α[28] - u[7])
68+
du[8] = α[29] * (xi(u[3]) * α[30] * y4_τ4 * y6_τ4 * y8_τ4 -
69+
u[4] * u[6] * u[8]) + α[31] * (α[32] - u[8])
70+
du[9] = α[33] * xi(u[3]) * α[34] * y4_τ5 * y6_τ5 * y8_τ5 +
71+
α[35] * (α[36] - u[9])
72+
du[10] = α[37] * u[9] - α[38] * u[10] * u[1] - α[39] * u[10]
73+
74+
nothing
75+
end
76+
77+
α = zeros(45)
78+
# Kinetic parameters
79+
α[1] = 83.0; α[2] = 5.0; α[3] = 6.6e14; α[4] = 3.0e11; α[5] = 0.4
80+
α[6] = 2.5e7; α[7] = 0.5e-12; α[8] = 2.3e9; α[9] = 0.052; α[10] = 0.15
81+
α[11] = 9.4e9; α[12] = 1.0e-15; α[13] = 1.2; α[14] = 2.7e16; α[15] = 2.0
82+
α[16] = 5.3e27; α[17] = 1.0; α[18] = 1.0e-18; α[19] = 2.7e16; α[20] = 2.0
83+
α[21] = 8.0e28; α[22] = 1.0; α[23] = 1.0e-19; α[24] = 5.3e33; α[25] = 16.0
84+
α[26] = 1.6e14; α[27] = 0.4; α[28] = 1.0e-18; α[29] = 8.0e32; α[30] = 16.0
85+
α[31] = 0.1; α[32] = 1.0e-18; α[33] = 1.7e30; α[34] = 3.0; α[35] = 0.4
86+
α[36] = 4.3e-22; α[37] = 0.85e7; α[38] = 8.6e11; α[39] = 0.043
87+
# Delays
88+
α[41] = 0.6; α[42] = 0.6; α[43] = 2.0; α[44] = 2.0; α[45] = 3.0
89+
90+
function h_hepatitis(p, t; idxs::Union{Nothing, Int} = nothing)
91+
α = p
92+
u0 = [2.9e-16, 0.0, 0.0, 0.0, α[18], α[23], α[28], α[32], α[36],
93+
α[37] * α[36] / α[39]]
94+
if idxs === nothing
95+
u0
96+
else
97+
u0[idxs]
98+
end
99+
end
100+
101+
prob = DDEProblem(f_hepatitis!, h_hepatitis, (0.0, 130.0), α;
102+
constant_lags = [0.6, 2.0, 3.0])
103+
```
104+
105+
## Reference Solution
106+
107+
```julia
108+
sol = solve(prob, MethodOfSteps(Rodas5P());
109+
reltol = 1e-12, abstol = 1e-20, dt = 1e-6)
110+
test_sol = TestSolution(sol)
111+
plot(sol; title = "Hepatitis B Solution")
112+
```
113+
114+
We also plot the first three components (virus dynamics) separately:
115+
116+
```julia
117+
plot(sol; idxs = [1, 2, 3], title = "Hepatitis B: Virus Dynamics")
118+
```
119+
120+
## High Tolerances
121+
122+
### Rosenbrock methods
123+
124+
```julia
125+
abstols = 1.0 ./ 10.0 .^ (4:7)
126+
reltols = 1.0 ./ 10.0 .^ (1:4)
127+
128+
setups = [Dict(:alg => MethodOfSteps(Rosenbrock23())),
129+
Dict(:alg => MethodOfSteps(Rodas4())),
130+
Dict(:alg => MethodOfSteps(Rodas5())),
131+
Dict(:alg => MethodOfSteps(Rodas5P()))]
132+
names = ["Rosenbrock23", "Rodas4", "Rodas5", "Rodas5P"]
133+
wp = WorkPrecisionSet(prob, abstols, reltols, setups;
134+
names = names, appxsol = test_sol, maxiters = Int(1e5), error_estimate = :final,
135+
dt = 1e-6)
136+
plot(wp; title = "Hepatitis B: Rosenbrock Methods (final error)")
137+
```
138+
139+
```julia
140+
wp = WorkPrecisionSet(prob, abstols, reltols, setups;
141+
names = names, appxsol = test_sol, maxiters = Int(1e5), error_estimate = :L2,
142+
dt = 1e-6)
143+
plot(wp; title = "Hepatitis B: Rosenbrock Methods (L2 error)")
144+
```
145+
146+
### SDIRK methods
147+
148+
```julia
149+
setups = [Dict(:alg => MethodOfSteps(TRBDF2())),
150+
Dict(:alg => MethodOfSteps(KenCarp4())),
151+
Dict(:alg => MethodOfSteps(Kvaerno4()))]
152+
names = ["TRBDF2", "KenCarp4", "Kvaerno4"]
153+
wp = WorkPrecisionSet(prob, abstols, reltols, setups;
154+
names = names, appxsol = test_sol, maxiters = Int(1e5), error_estimate = :final,
155+
dt = 1e-6)
156+
plot(wp; title = "Hepatitis B: SDIRK Methods (final error)")
157+
```
158+
159+
```julia
160+
wp = WorkPrecisionSet(prob, abstols, reltols, setups;
161+
names = names, appxsol = test_sol, maxiters = Int(1e5), error_estimate = :L2,
162+
dt = 1e-6)
163+
plot(wp; title = "Hepatitis B: SDIRK Methods (L2 error)")
164+
```
165+
166+
### All stiff methods comparison
167+
168+
```julia
169+
setups = [Dict(:alg => MethodOfSteps(Rosenbrock23())),
170+
Dict(:alg => MethodOfSteps(Rodas5P())),
171+
Dict(:alg => MethodOfSteps(TRBDF2())),
172+
Dict(:alg => MethodOfSteps(KenCarp4()))]
173+
names = ["Rosenbrock23", "Rodas5P", "TRBDF2", "KenCarp4"]
174+
wp = WorkPrecisionSet(prob, abstols, reltols, setups;
175+
names = names, appxsol = test_sol, maxiters = Int(1e5), error_estimate = :final,
176+
dt = 1e-6)
177+
plot(wp; title = "Hepatitis B: All Stiff Methods (final error)")
178+
```
179+
180+
## Low Tolerances
181+
182+
```julia
183+
abstols = 1.0 ./ 10.0 .^ (8:11)
184+
reltols = 1.0 ./ 10.0 .^ (5:8)
185+
186+
setups = [Dict(:alg => MethodOfSteps(Rosenbrock23())),
187+
Dict(:alg => MethodOfSteps(Rodas4())),
188+
Dict(:alg => MethodOfSteps(Rodas5P())),
189+
Dict(:alg => MethodOfSteps(TRBDF2())),
190+
Dict(:alg => MethodOfSteps(KenCarp4()))]
191+
names = ["Rosenbrock23", "Rodas4", "Rodas5P", "TRBDF2", "KenCarp4"]
192+
wp = WorkPrecisionSet(prob, abstols, reltols, setups;
193+
names = names, appxsol = test_sol, maxiters = Int(1e5), error_estimate = :final,
194+
dt = 1e-6)
195+
plot(wp; title = "Hepatitis B: Low Tolerances (final error)")
196+
```
197+
198+
```julia
199+
wp = WorkPrecisionSet(prob, abstols, reltols, setups;
200+
names = names, appxsol = test_sol, maxiters = Int(1e5), error_estimate = :L2,
201+
dt = 1e-6)
202+
plot(wp; title = "Hepatitis B: Low Tolerances (L2 error)")
203+
```
204+
205+
```julia, echo = false
206+
using SciMLBenchmarks
207+
SciMLBenchmarks.bench_footer(WEAVE_ARGS[:folder], WEAVE_ARGS[:file])
208+
```

0 commit comments

Comments
 (0)