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

4 views (last 30 days)
Theresa
Theresa on 18 Feb 2013
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);

Answers (1)

Babak
Babak on 19 Feb 2013
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)

Community Treasure Hunt

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

Start Hunting!