Edit: The feature that is critical is the notch, which all the small dt results and the stiff results show, but not the high dt results. This notch is a response to the driving function I, which looks like this:
ODEs:
function wk5(dP, P, params, t)
#=
The 5-Element WK with serial L
set Ls = 0 for 4-Element WK parallel
Formulation for DifferentialEquations.jl
P: solution vector (pressures p1 and p2)
params: parameter array
[ Rc, Rp, C, Lp, Ls ]
I need to find a way to tranfer the function name as well
for the time being we have to have the function in "I"
=#
# Split parameter array:
Rc, Rp, C, Lp, Ls = params
dP[1] = (
-Rc / Lp * P[1]
+ (Rc / Lp - 1 / Rp / C) * P[2]
+ Rc * (1 + Ls / Lp) * didt(I, t)
+ I(t) / C
)
dP[2] = -1 / Rp / C * P[2] + I(t) / C
return dP[1], dP[2]
end
Driving function (flow rate over time):
# Generic Input Waveform
# max volume flow in ml/s
max_i = 425
# min volume flow in m^3/s
min_i = 0.0
T = 0.9
# Syst. Time in s
systTime = 2 / 5 * T
# Dicrotic notch time
dicrTime = 0.02
function I_generic(t)
# implicit conditional using boolean multiplicator
# sine waveform
(
(max_i - min_i) * sin(pi / systTime * (t % T))
* (t % T < (systTime + dicrTime) )
+ min_i
)
end
Central Difference (I used to use Calculus, but put this quick hack in to keep it consistent with the Matlab code:
function didt(I, t)
dt = 1e-3;
didt = (I(t+dt) - I(t-dt)) / (2 * dt);
return didt
end
And finally, solving the thing:
# Initial condition and time span
P0 = [0.0, 0.0]
tspan = (0.0, 30.0)
# Set parameters for Windkessel Model
Rc = 0.033
Rp = 0.6
C = 1.25
# L for serial model!
Ls = 0.01
# L for parallel
Lp = 0.02
p5 = [Rc, Rp, C, Lp, Ls]
I = I_generic
problem = ODEProblem(wk5, P0, tspan, p5)
dtmax = 1e-4
@time solutionTsit = solve(problem)
@time solutionTsitLowdt = solve(problem, Tsit5(), dtmax=dtmax)
@time solutionBS3 = solve(problem, BS3())
@time solutionBS3Lowdt = solve(problem, BS3(), dtmax=dtmax);
@time solutionDP5 = solve(problem,DP5())
@time solutionDP5Lowdt = solve(problem,DP5(), dtmax=dtmax);
@time solutionStiff = solve(problem, alg_hints=[:stiff]);
The oscillations in this case are low frequency, like this (Tsit5 default):
But I also get something like this on the other model. But this may actually be a different problem, since this doesn’t have the massive notch, and shows a high frequency oscillation.
This seems to be what I expect of a high order non-stiff solver, but, again, Matlab ode45 doesn’t show that.
Edit: Chris’ hunch that it is the time step seems to be correct, the Matlab solution has around 5985 steps for the 30 seconds (5278 for the stiff solver ode15s), while Julia has 193, 901, 173 for Tsit5, BS3, and DP5, respectively, and 4078 for :stiff.
It is to be expected that the time step chosen for dI/dt (function didt) is critical and that the notch can get missed (and it does). However, with only 200 time steps over > 20 cycles, it is not surprising that small scale features in the driving function are smoothed out (in particular since it’s dI/dt driving the inertia terms in the pressure response). The real driving functions will have a notch like that, but it will be smoothed out and differentiable, but still smaller than the chosen time step for the integration - so the problem shows on some of the measures curves as well (not as pronounced, and - worryingly - not as reproducible).
Am I correct in assuming that DP5 and ode45 are using the same algorithm?