ODE: how to plot result and piecewise function?

2 次查看（过去 30 天）
Marco Sammito
Marco Sammito 2016-12-3
评论： Star Strider 2016-12-3
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

请先登录，再回答此问题。

回答（1 个）

Star Strider
Star Strider 2016-12-3
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.
  2 个评论
Marco Sammito
Marco Sammito 2016-12-3
How can I stop the integration, as you said, and use the final values as new initial values?
Star Strider
Star Strider 2016-12-3
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.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Programming 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by