How to write a code for a bounded periodic function?

Question

Anne Brinkman 2019-6-3

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/465318-how-to-write-a-code-for-a-bounded-periodic-function

回答： Torsten 2019-6-4

close all
syms W(t);
W(t) = piecewise(3/12<=t<=8/12, 1,t>=1, W(t-1) ,0);
t(1) = 0;
y1(1) = 600;
y2(1) = 1000;
h = 0.5;
t_end = 5;
t = 0:h:t_end;
a = 2;
c = 3;
alfa = 10^-3;
gamma = 6*10^-3;
f = @(t,y1,y2) (a-alfa*y2)*y1-W(t)*y1;
g = @(t,y1,y2) (-c+gamma*y1)*y2;
for i = 1:(length(t)-1)
     
    k1 = f(t(i),y1(i),y2(i));
    l1 = g(t(i),y1(i),y2(i));
    k2 = f(t(i)+h/2,(y1(i)+0.5*k1*h),(y2(i)+(0.5*l1*h)));
    l2 = g(t(i)+h/2,(y1(i)+0.5*k1*h),(y2(i)+(0.5*l1*h)));
    
    k3 = f(t(i)+h/2,(y1(i)+0.5*k2*h),(y2(i)+(0.5*l2*h)));
    l3 = g(t(i)+h/2,(y1(i)+0.5*k2*h),(y2(i)+(0.5*l2*h)));
    
    k4 = f(t(i)+h,(y1(i)+k3*h),(y2(i)+l3*h));
    l4 = g(t(i)+h,(y1(i)+k3*h),(y2(i)+l3*h));
    
    y1(i+1) = y1(i) + (k1  +2*k2   +2*k3    +k4)*(h/6);
    y2(i+1) = y2(i) + (l1  +2*l2   +2*l3    +l4)*(h/6);
    
    ystar=[y1', y2'];
    
end

This code is for the implementation of a predator-prey model/ volterra-lotka model using the time integration method RK4. There is a t dependent (W(t)) variable in the function. How can these bounded values be implemented on matlab? Especially the W(t-1) part, this part will look at the previous boundaries.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

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

Answer 1

Torsten 2019-6-4

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/465318-how-to-write-a-code-for-a-bounded-periodic-function#answer_377777

t(1) = 0;
y1(1) = 600;
y2(1) = 1000;
h = 0.5;
t_end = 5;
t = 0:h:t_end;
a = 2;
c = 3;
alfa = 10^-3;
gamma = 6*10^-3;
f = @(t,y1,y2) (a-alfa*y2)*y1-W(t)*y1;
g = @(t,y1,y2) (-c+gamma*y1)*y2;
for i = 1:(length(t)-1)
    k1 = f(t(i),y1(i),y2(i));
    l1 = g(t(i),y1(i),y2(i));
    k2 = f(t(i)+h/2,(y1(i)+0.5*k1*h),(y2(i)+(0.5*l1*h)));
    l2 = g(t(i)+h/2,(y1(i)+0.5*k1*h),(y2(i)+(0.5*l1*h)));
    k3 = f(t(i)+h/2,(y1(i)+0.5*k2*h),(y2(i)+(0.5*l2*h)));
    l3 = g(t(i)+h/2,(y1(i)+0.5*k2*h),(y2(i)+(0.5*l2*h)));
    k4 = f(t(i)+h,(y1(i)+k3*h),(y2(i)+l3*h));
    l4 = g(t(i)+h,(y1(i)+k3*h),(y2(i)+l3*h));
    y1(i+1) = y1(i) + (k1  +2*k2   +2*k3    +k4)*(h/6);
    y2(i+1) = y2(i) + (l1  +2*l2   +2*l3    +l4)*(h/6);
    ystar=[y1', y2'];
end
function Wvalue = W(t)
  W0 = some value;
  Wvalue = 0.0;
  tt = t - floor(t);
  if tt >= 3/12 && tt <= 8/12
    Wvalue = W0;
  end
end