Why doesn't my ODE solution to a simple exponential decay model go to zero?

Hi,

I'm trying to model a simple exponential decay (this is part of a much larger model, but the issue traces back to this so I've simplified for clarity). The model should decay to zero, but it never gets to zero. Instead, the solution bounces around zero and goes negative sometimes. Is there anything I can do to eliminate this behavior, and get a simple exponential decay model to decay to zero?

I've tried changing the tolerance and that doesn't work. I've also tried using the NonNegative function in my larger model, and I get very, very different results. In this small example, the solution goes to zero, then goes up to a very small number again (why doesn't it stay at zero)?

Because this is part of a larger model, I can't just use the exact ODE solution.

Thanks for any help you can provide.

Sample code:

function [t y] = simtest()

options = odeset('AbsTol', 1e-12, 'RelTol', 1e-12);%
[t y] = ode45(@test, [0 100], 1,options);
figure()
subplot(1,2,1)
plot(t, y(:,1))
subplot(1,2,2)
plot(t, y(:,1))
set(gca, 'YScale', 'log');

function dydt = test(t,y)
dydt = -y(1);

Answer 1

[Babak]

The problem comes from this line

 set(gca, 'YScale', 'log');

Even if you do your own log10 for the y axis data by defining

 y_log10 = log10(y);

then

 plot(t,y_log10);

you still see similar behavior which is because of the maximum error resolution of MATLAB. that is the nuber of digits it saves the data. 10^(-15) is really a very high resolution which is higher resolution than what you have assigned (10^-12) so you will see gitter for log10(y)<10^(-12)