Using ode45 in Matlab

Question

I'm trying to simulate the time behavior for a physical process governed by a system of ODEs. When I switch the width of the input pulse from 20 to 19, there is no depletion of the y(1) state, which doesn't make sense physically. What am I doing wrong? Am I using ode45 incorrectly?

function test

width = 20;
center = 100;

tspan = 0:0.1:center+50*(width/2);

[t,y] = ode45(@ODEsystem,tspan,[1 0 0 0]);

plot(t,y(:,1),'k*',t,y(:,2),'k:',t,y(:,3),'k--',t,y(:,4),'k');
hold on;
axis([center-3*(width/2) center+50*(width/2) -0.1 1.1])
xlabel('Time')
ylabel('Relative values')
legend({'y1','y2','y3','y4'});

    function dy = ODEsystem(t,y)
        k1 = 0.1;
        k2 = 0.000333;
        k3 = 0.1;

        dy = zeros(size(y));

        % rectangular pulse
        I = rectpuls(t-center,width);

        % ODE system
        dy(1) = -k1*I*y(1);
        dy(2) = k1*I*y(1) - k2*y(2);
        dy(3) = k2*y(2) - k3*I*y(3);
        dy(4) = k3*I*y(3);
    end
end

score 2 · Accepted Answer · edited May 23 '17 at 11:43

You are changing the parameters of your your ODEs discontinuously in time. This results in a very stiff system and less accurate, or even completely wrong, results. In this case, because the your ODE is so simple when I = 0, an adaptive solver like ode45 will take very large steps. Thus, there's a high probability that it will step right over the region where you inject the impulse and never see it. See my answer here if you're confused as to why the code in your question misses the pulse even though you've specified tspan to have (output) steps of just 0.1.

In general it is a bad idea to have any discontinuities (if statements, abs, min/max, functions like rectpuls, etc.) in your integration function. Instead, you need to break up the integration and calculate your results piecewise in time. Here's a modified version of your code that implements this:

function test_fixed

width = 19;
center = 100;

t = 0:0.1:center+50*(width/2);
I = rectpuls(t-center,width); % Removed from ODE function, kept if wanted for plotting

% Before pulse
tspan = t(t<=center-width/2);
y0 = [1 0 0 0];
[~,out] = ode45(@(t,y)ODEsystem(t,y,0),tspan,y0); % t pre-calculated, no need to return
y = out;

% Pulse
tspan = t(t>=center-width/2&t<=center+width/2);
y0 = out(end,:); % Initial conditions same as last stage from previous integration
[~,out] = ode45(@(t,y)ODEsystem(t,y,1),tspan,y0);
y = [y;out(2:end,:)]; % Append new data removing identical initial condition

% After pulse
tspan = t(t>=center+width/2);
y0 = out(end,:);
[~,out] = ode45(@(t,y)ODEsystem(t,y,0),tspan,y0);
y = [y;out(2:end,:)];

plot(t,y(:,1),'k*',t,y(:,2),'k:',t,y(:,3),'k--',t,y(:,4),'k');
hold on;
axis([center-3*(width/2) center+50*(width/2) -0.1 1.1])
xlabel('Time')
ylabel('Relative values')
legend({'y1','y2','y3','y4'});

    function dy = ODEsystem(t,y,I)
        k1 = 0.1;
        k2 = 0.000333;
        k3 = 0.1;

        dy = zeros(size(y));

        % ODE system
        dy(1) = -k1*I*y(1);
        dy(2) = k1*I*y(1) - k2*y(2);
        dy(3) = k2*y(2) - k3*I*y(3);
        dy(4) = k3*I*y(3);
    end
end

See also my answer to a similar question.

Thanks! Since it's stiff, shouldn't a different ode solver be used, or does it become non-stiff when you break it up piecewise? Also, is there an alternative way to do it, since it gives me an error if I try simulating a vanishingly small width (a delta function)? — user4038089, May 01 '15 at 00:52
No, a stiff solver may be slightly more reliable, but it's not guaranteed to not suffer from the same issues when your system has such discontinuities. ODEs fundamentally require some degree of smoothness so it's also wrong mathematically. What's wrong with the method in this answer? I don't know what you mean by "vanishingly small width." You're dealing with discrete floating point values so the width must be some finite value. — horchler, May 01 '15 at 01:33

Using ode45 in Matlab

1 Answers1