ode23, ode45 acting weird - why?
13 views (last 30 days)
Show older comments
William Rose
on 14 Apr 2021
Commented: David Goodmanson
on 18 Apr 2021
I am integrating a pair of simple non-stiff first order differential equations. The output, using ode45(), is as expected, when tspan=[0 41]: see plot below.
But when tspan=[0 40] (versus [0 41] above), there is a glitch around t=7 sec. It seems ode45 took a few unusually long steps around t=7. See plot below.
So I tried ode23(). The routine also works fine when tspan=[0 41]: see plot below, which basically matches plot 1 above, although the step sizes are a bit different.
However, ode23() completely fails when tspan=[0 40]. There is no error message, but the variables never leave their initial conditions, and ode23() takes only 10 steps to span the interval. See plot below. Note that the y-axis units on the bottom plot are x.
I am attaching the code, with ode45() and tspan=[0 41]. I never change function dydt01.m. The only things I changed were in the main program, cvmodel01.m: tspan=[0 41] or [0 40], and ode45() or ode23().
Why is this happening? The error with ode45() is subtle, so one could be unaware of it, in a more realistic simulation.
ode45 and ode23 use adaptive stepsize. Therefore they make initial guesses for stepsize. The initial guesses must be different when tEnd=40 than when tEnd=41. Does this difference lead to a glitch at t=7, in the case of ode45()?
In the case of ode23(), with tspan=[0 40], something particularly unlucky happens. What makes this simulation "go" is the beating of the heart, which is represented as a current source with output A*sin^2(pi*t).* ode23() decides to use step sizes of exactly 4 seconds, which seems surprisingly large, so it happens to land on the zeros of the current waveform every time, and it does not notice that there were four "missed" heartbeats in those 4 seconds. It is a form of aliasing. It is odd and unfortunate that ode23() takes such inappropriately large step sizes.
*This is not a realistic model of the heart. It is model 1 for teaching Matlab to biomedical engineering students.
0 Comments
Accepted Answer
David Goodmanson
on 14 Apr 2021
Edited: David Goodmanson
on 14 Apr 2021
Hi William,
Very interesting behavior. Different look but still impressive when the span ends at 40.1. I think the reason may be that your initial conditions are a little to perfect. With
y0(1)=Pmf*Ca; %arterial volume (mL)
y0(2)=Pmf*Cv; %venous volume (mL)
and with
Pa=y(1)/Ca;
Pv=y(2)/Cv;
for the initial conditions, the Ca's and Cv's cancel out, meaning that Pa-Pv = 0 and qR = 0. Since the sine term is also zero at t=0, dy/dt = 0 initially which may be why ode23 flatlines. Things work out fine for ode45 and
y0(1)=Pmf*Ca*.99 %arterial volume (mL)
y0(2)=Pmf*Cv; %venous volume (mL)
and a maybe even better approach is to keep th ic's the same as before and just use
tspan = [.01 40]
which kicks in the sine term at the start.
As to how those original ic's can affect the result that far into the time record, I don't know (yet?).
0 Comments
More Answers (1)
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!