# How to write a code for a bounded periodic function?

3 visualizaciones (últimos 30 días)
Anne Brinkman el 3 de Jun. de 2019
Respondida: Torsten el 4 de Jun. de 2019
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 comentariosMostrar -2 comentarios más antiguosOcultar -2 comentarios más antiguos

Iniciar sesión para comentar.

Torsten el 4 de Jun. de 2019
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 comentariosMostrar -2 comentarios más antiguosOcultar -2 comentarios más antiguos

Iniciar sesión para comentar.

### Categorías

Más información sobre Ordinary Differential Equations en Help Center y File Exchange.

R2017b

### Community Treasure Hunt

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

Start Hunting!

Translated by