Hi, i am trying plot my thesis graph using dde23 but i keep on getting this error, derivative and history function have different lengths.. HELP PLEASE!
3 views (last 30 days)
Show older comments
function Hassan
lags = 0.2;
tspan = [0,60];
sol = dde23(@ddehs, lags, @hshist, tspan);
time = sol.x;
S = sol.y(1,:);
I = sol.y(2,:);
V = sol.y(3,:);
figure;
plot(time,S,'r')
xlabel('Time (Years)'); ylabel('Susceptible Human Population');
legend('S_h')
figure;
plot(time,I,'b')
xlabel('Time (Years'); ylabel('Infected Human Population');
legend('I_h')
figure;
plot(time,V,'b')
xlabel('Time (Years'); ylabel('Vaccinated Human Population');
legend('V_h')
end
function dydt = ddehs(~,y,z)
Kh = 10000; muH = 0.0066;
aH = 0.02; Bh = 0.3;
muD = 0.08; tau = 0.2;
wH = 1;
ylag1 = z(:,1);
dydt = [ muH*Kh - (muH + aH)*y(1) - Bh*y(1)*ylag1*exp(-muH*tau)
Bh*y(1)*ylag1*exp(-muD*tau) - (muH + wH)*y(2)
aH*y(1) - muH*y(3) ];
end
function s = hshist(~)
s = ones(3,1);
end
0 Comments
Accepted Answer
Torsten
on 26 Jan 2023
Moved: Torsten
on 26 Jan 2023
z(:,1) is a 3x1 vector with z(1,1) = y1(t-0.02), z(2,1) = y2(t-0.02) and z(3,1) = y3(t-0.02),
You have to decide which of the three delay terms you want to insert in dydx(1) and dydx(2).
9 Comments
Torsten
on 1 Feb 2023
Edited: Torsten
on 1 Feb 2023
lags = [0.2,2.0,6.0,12.0];
time = 0:0.1:60;
for i = 1:numel(lags)
tspan = [0,60];
sol = dde23(@ddehs, lags(i), @hshist, tspan);
yint = deval(sol,time);
S(i,:) = yint(1,:);
I(i,:) = yint(2,:);
V(i,:) = yint(3,:);
end
figure;
plot(time,S,'r')
xlabel('Time (Years)'); ylabel('Susceptible Human Population');
legend('S_h')
figure;
plot(time,I,'b')
xlabel('Time (Years)'); ylabel('Infected Human Population');
legend('I_h')
figure;
plot(time,V,'b')
xlabel('Time (Years)'); ylabel('Vaccinated Human Population');
legend('V_h')
function dydt = ddehs(~,y,z)
Kh = 10000; muH = 0.0066;
aH = 0.02; Bh = 0.3;
muD = 0.08; tau = 0.2;
wH = 1;
ylag1 = z(2,1);
dydt = [ muH*Kh - (muH + aH)*y(1) - Bh*y(1)*ylag1*exp(-muH*tau)
Bh*y(1)*ylag1*exp(-muD*tau) - (muH + wH)*y(2)
aH*y(1) - muH*y(3) ];
end
function s = hshist(~)
s = ones(3,1);
end
More Answers (1)
See Also
Categories
Find more on Introduction to Installation and Licensing in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!