ODE: how to plot result and piecewise function?

Hi. I have the following problem

where

I am able to plot t versus x. What I would like to is to plot t versus x but also t versus f (in the same graph with two y axes. To do so the command is plotyy, but I am not even able to plot t versus f on a different graph). I tried

plot(t,f(t,A,B,C))

after

    fval = 0;
  end

but that did not work.

function piecewise()
x0    = 1;
tspan = linspace(0,3,1000);
[T,X] = ode45(@DE, tspan, x0);
plot(T,X)
end
%
function dX = DE(t,x)
A = log(2);
B = log(3);
C =1;
dX = -13.925*A*B*C*x + f(t,A,B,C);
end
%
function fval = f(t,A,B,C)
if     (t <= log(2))
    fval = exp(-t);
elseif (t > log(2)) && (t <= log(3))
    fval = 4;
else
    fval = 0;
end
end

Answer 1

[Star Strider]

This is my approach:

tv = linspace(0, 3, 1000);
x0 = 1;

A = log(2);
B = log(3);
C = 1;

f = @(t,A,B) (exp(-C*t) .* (t <= A)) + (4*((A < t) & (t <= B)));

ODE = @(t,x,A,B,C,f)  -13.925*A*B*C*x + f(t,A,B);

[t,x] = ode15s(@(t,x)ODE(t,x,A,B,C,f), tv, x0);

ft = f(t,A,B);

figure(1)
plotyy(t,x,  t,ft)
grid
xlabel('\itt\rm')
ylabel('\itx\rm')
text([A B], [1 1]*0.7, {'\leftarrowA', '\leftarrowB'})
legend('x(t)', 'f(t)')

It uses ‘logical indexing’ to do the conditional operations. It is normally not recommended to integrate across discontinuities (the usual procedure is to stop the integration, use the previous final values as the new initial conditions, and integrate stepwise between discontinuities to the end), but it seems that we can get away with it here because the steps are ‘almost’ continuous.

[Marco Sammito]

How can I stop the integration, as you said, and use the final values as new initial values?

[Star Strider]

This works:

tv = linspace(0, 3, 1000);
x0 = 1;

A = log(2);
B = log(3);
C = 1;

ODE = @(t,x,A,B,C,f)  -13.925*A*B*C*x + f(t);

ti = tv(tv <= A);                                               % Time Interval
fi = @(t) exp(-C*t);                                            % Function For This Interval
[t{1},x{1}] = ode15s(@(t,x)ODE(t,x,A,B,C,fi), ti, x0);

x0 = x{1}(end);
ti = tv((A < tv) & (tv <= B));
fi = @(t) 4;
[t{2},x{2}] = ode15s(@(t,x)ODE(t,x,A,B,C,fi), ti, x0);

x0 = x{2}(end);
ti = tv(tv > B);
fi = @(t) 0;
[t{3},x{3}] = ode15s(@(t,x)ODE(t,x,A,B,C,fi), ti, x0);

tt = cell2mat(t');
xt = cell2mat(x');

ft = f(tv,A,B);

figure(1)
plotyy(tt,xt,  tv,ft)
grid
xlabel('\itt\rm')
ylabel('\itx\rm')
text([A B], [1 1]*0.7, {'\leftarrowA', '\leftarrowB'})
legend('x(t)', 'f(t)')

It’s easier to break these into three separate sections than to use a loop with all the associated logic tests to choose the various functions and time intervals. You can certainly do that if you want to, but it’s even less efficient.

This produces the same plot as my earlier — and much more efficient — code.