42 views (last 30 days)

Show older comments

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.

David Goodmanson
on 14 Apr 2021 at 21:59

Edited: David Goodmanson
on 14 Apr 2021 at 22:16

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

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

Start Hunting!