MATLAB Answers

0

How to write a code for a bounded periodic function?

Anne Brinkman さんによって質問されました 2019 年 6 月 3 日
最新アクティビティ 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.IMG_20190603_151631.jpg

  0 件のコメント

サインイン to comment.

1 件の回答

Torsten
回答者: Torsten
2019 年 6 月 4 日
 採用された回答

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

  0 件のコメント

サインイン to comment.



Translated by