ode23, ode45 acting weird - why?

13 views (last 30 days)
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.

Accepted Answer

David Goodmanson
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?).

More Answers (1)

William Rose
William Rose on 15 Apr 2021
@David Goodmanson, Thank you for your thoughts and for running the code. I am glad to know that you got similar results. As for "too perfect" initial conditions, yes, the system is in equilibrium initially. This is a common and natural starting point, like a mass on a spring that is not stretched or compressed initally, but with a cyclic external force. Yes, tweaking the start time or the end time or the initial conditions a bit can make it work right. It is interesting and strange that the most natural start time, zero, and the most natural start point, equilibrium, can lead to a total failure by ode23 to find the correct solution, and a more subtle funny result from ode45. I will take it as a lesson to try different solvers, to try slightly different end or start times, and to watch out for weird behavior.
Thanks again.
  3 Comments
William Rose
William Rose on 18 Apr 2021
Thank you @Paul for looking at this and for your insight. I didn't know about MaxStep or InitialStep options. Those are both good ideas!

Sign in to comment.

Tags

Products


Release

R2018b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!