How to write a code for a bounded periodic function?
4 次查看(过去 30 天)
显示 更早的评论
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](https://www.mathworks.com/matlabcentral/answers/uploaded_files/222515/IMG_20190603_151631.jpeg)
![IMG_20190603_151631.jpg](https://www.mathworks.com/matlabcentral/answers/uploaded_files/222515/IMG_20190603_151631.jpeg)
0 个评论
采纳的回答
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 个评论
更多回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Numerical Integration and Differential Equations 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!