ode23, ode45 acting weird - why?

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.

Answer 1

[David Goodmanson]

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

Answer 2

[William Rose]

@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.

[Paul]

You absolutely should not have to try different start and end times, which are fundamental to the problem you're trying to solve. Definitely be on the lookout for weird behavior, especially with these variable step solvers.

I don't have a lot of experience with ode45, et. al., but I do with Simulink and I never let the software choose the max step size. I try to set the max step size based on what I think the dynamics are of the system, look at the results, see if they seem to make sense, iterate as needed. In this problem the forcing function has a period of T = 1, so perhaps a max step size of 0.1 would be appropriate

[t,y]=ode45(@dydt01,[0 40],y0,odeset('MaxStep',0.1));

Solution looks good, and note that the plot of CO looks nice and smooth. Another useful parameter is the initial step size, which overcomes the issue with ode23 that is probably related to @David Goodmanson's observiation about the initial conditions and initial input:

[t,y]=ode23(@dydt01,[0 40],y0,odeset('InitialStep',0.1));

For this problem, it's probably best to set both. Check doc odeset for more info and discussion on these and other error control parameters.

[William Rose]

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!

[David Goodmanson]

Yes they are.