Attempt at making a Runge Kutta 4 script
1 次查看(过去 30 天)
显示 更早的评论
Trying to create a simulation a of population using RK4 method, but I'm struggling to get it to run properly. The idea is to get y1 and y2 with 11 numbers, which will represent population amounts.
This is the main script:
close all
clear
alpha = 10^-3;
gamma = 6e-3;
c = 3;
a = 2;
Dt = 0.5;
lambda = sqrt(-6);
y1 = 600;
y2 = 1000;
t = 1:Dt:6;
for i = 1:10
k11 = Dt*dy1(i,y1,y2);
k12 = Dt*dy2(i,y1,y2);
k21 = Dt*(dy1(i,y1,y2) + 0.5*k11);
k22 = Dt*(dy2(i,y1,y2) + 0.5*k12);
k31 = Dt*(dy1(i,y1,y2) + 0.5*k21);
k32 = Dt*(dy2(i,y1,y2) + 0.5*k22);
k41 = Dt*(dy1(i,y1,y2) + k31);
k42 = Dt*(dy2(i,y1,y2) + k32);
y1(i+1) = y1(i) + (1/6)*(k11 + k21 + k31 + k41);
y2(i+1) = y2(i) + (1/6)*(k12 + k22 + k32 + k42);
end
dy1, dy2 and W are these funtions:
function g = dy1(t,y1,y2)
alpha = 10^-3;
gamma = 6e-3;
c = 3;
a = 2;
lambda = sqrt(-6);
%for t = 1:11
g = (a - alpha.*y2).*y1 - W(t).*y1;
%end
end
function h = dy2(t,y1,y2)
alpha = 10^-3;
gamma = 6e-3;
c = 3;
a = 2;
%y2 = a/alpha;
%y1 = c/gamma;
Dt = 0.5;
lambda = sqrt(-6);
y1 = 600;
y2 = 1000;
h = (-c + gamma*y1)*y2;
end
function f = W(t)
f = [];
if 0 <= t && t < (3/12)
f = 0;
end
if (3/12) <= t && t <= (8/12)
f = 1;
end
if (8/12) < t && t < 1
f = 0;
end
if t >= 1
t = t - 1;
end
end
No watter how i try to change the script or the functions, I keep getting differnt kinds of errors. I feel like the general idea is there, but I just don't understand what exactly I'm doing that matlab doesn't like. Any help or insight would be greatly appreciated.
2 个评论
Jan
2019-6-11
编辑:Jan
2019-6-11
The readers cannot guess, what you want to achieve. Why do you restrict the values of t at the end of the W(t) function? t is not used afterwards anymore, so you can omit the change to t-1.
You define t = 1:Dt:6, but do not use these values later on. I guess you want to replace the loop for i = 1:10 by
tVec = 1:Dt:6,
for k = 1:numel(tVec)
t = tVec(k);
k11 = Dt*dy1(t,y1,y2);
% ^ t instead of i
end
"I keep getting differnt kinds of errors" - then post the code and a copy of the complete error message. It is easier to fix an error than to guess, what the error is.
回答(0 个)
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!