Runge-kutta 4th order Method Help
显示 更早的评论
% if use ode23 work
% but usin rk4 dont, please help!
% I use matlab to test and implement routines in fortran 90
% Tks!
h = 1; % step day
ti = 1; % first day
tf = 365; % one year
t = ti:h:tf;
y = zeros(1, length(t));
R1(1) = 1.0; % initial resource 1
R2(1) = 0.3; % initial resource 2
C1(1) = 1.0; % initial consumer 1
C2(1) = 1.0; % initial consumer 2
P(1) = 1.0; % initial predator
% prey-predator model
% P
% / \
% C1 C2
% | |
% R1 R2
% Paremeters
p = 0.50;
xc = 0.15;
xp = 0.08;
yc = 2.30;
yp = 1.70;
r0 = 0.25;
c0 = 0.50;
%% Equations
% R1
F_R1 = (R1.*(1-R1))-(xc*yc*((C1.*R1)./(R1+r0)));
% R2
F_R2 = (R2.*(1-R2))-(xc*yc*((C2.*R2)./(R2+r0)));
% C1
F_C1 = (xc*yc*((C1.*R1)./(R1+r0)))-((p.*C1./(p.*C1+(1-p).*C2)).*xp.*yp.*((P.*C1)/(c0+C1)))-xc.*C1;
% C2
F_C2 = (xc*yc*((C2.*R2)./(R2+r0)))-((((1-p).*C2)./(p.*C1+(1-p).*C2)).*xp.*yp.*((P.*C2)./(c0+C2)))-xc.*C2;
% P
F_P = ((p*C1./(p*C1+(1-p)*C2)).*xp.*yp.*((P.*C1)/(c0+C1)))+((((1-p)*C2)./(p*C1+(1-p)*C2)).*xp.*yp.*((P.*C2)./(c0+C2)))-xp.*P;
for i=1:(length(t)-1)
t(i+1)=t(i)+h;
kR1_1 = F_R1(t(i),R1(i));
kR1_2 = F_R1(t(i)+0.5*h,R1(i)+0.5*h*kR1_1);
kR1_3 = F_R1((t(i)+0.5*h),(R1(i)+0.5*h*kR1_2));
kR1_4 = F_R1((t(i)+h),(R1(i)+kR1_3*h));
R1(i+1) = R1(i) + (1/6)*(kR1_1+2*kR1_2+2*kR1_3+kR1_4)*h; % main equation
kR2_1 = F_R2(t(i),R2(i));
kR2_2 = F_R2(t(i)+0.5*h,R2(i)+0.5*h*kR2_1);
kR2_3 = F_R2((t(i)+0.5*h),(R2(i)+0.5*h*kR2_2));
kR2_4 = F_R2((t(i)+h),(R2(i)+kR2_3*h));
R2(i+1) = R2(i) + (1/6)*(kR2_1+2*kR2_2+2*kR1_3+kR2_4)*h;
kC1_1 = F_C1(t(i),C1(i));
kC1_2 = F_C1(t(i)+0.5*h,C1(i)+0.5*h*kC1_1);
kC1_3 = F_C1((t(i)+0.5*h),(C1(i)+0.5*h*kC1_2));
kC1_4 = F_C1((t(i)+h),(C1(i)+kC1_3*h));
C1(i+1) = C1(i) + (1/6)*(kC1_1+2*kC1_2+2*kC1_3+kC1_4)*h;
kC2_1 = F_C2(t(i),C2(i));
kC2_2 = F_C2(t(i)+0.5*h,C2(i)+0.5*h*kC2_1);
kC2_3 = F_C2((t(i)+0.5*h),(C2(i)+0.5*h*kC2_2));
kC2_4 = F_C2((t(i)+h),(C2(i)+kC2_3*h));
C2(i+1) = C2(i) + (1/6)*(kC2_1+2*kC2_2+2*kC2_3+kC2_4)*h;
kP_1 = F_P(t(i),P(i));
kP_2 = F_P(t(i)+0.5*h,P(i)+0.5*h*kP_1);
kP_3 = F_P((t(i)+0.5*h),(P(i)+0.5*h*kP_2));
kP_4 = F_P((t(i)+h),(P(i)+kP_3*h));
P(i+1) = P(i) + (1/6)*(kP_1+2*kP_2+2*kP_3+kP_4)*h;
end
采纳的回答
James Tursa
2021-2-3
编辑：James Tursa
2021-2-3
All of your F's need to be functions that operate on the current full states. E.g., take this one
F_C1 = (xc*yc*((C1.*R1)./(R1+r0)))-((p.*C1./(p.*C1+(1-p).*C2)).*xp.*yp.*((P.*C1)/(c0+C1)))-xc.*C1;
It is obvious that the C1 derivative depends on the current states of C1, C2, R1, and P. So you need a function that takes these as inputs to compute the derivative. E.g., something like this
F_C1 = @(t,C1,C2,R1,P)(xc*yc*((C1.*R1)./(R1+r0)))-((p.*C1./(p.*C1+(1-p).*C2)).*xp.*yp.*((P.*C1)/(c0+C1)))-xc.*C1;
So modify your code to turn all of those F's into function handles as shown above.
Then you need to realize that all of your variables are interdependent, so using predictor methods will need to predict all of the states simultanously, not one at a time as you are coding. E.g., take these lines in your current code for C1:
kC1_1 = F_C1(t(i),C1(i));
kC1_2 = F_C1(t(i)+0.5*h,C1(i)+0.5*h*kC1_1);
kC1_3 = F_C1((t(i)+0.5*h),(C1(i)+0.5*h*kC1_2));
kC1_4 = F_C1((t(i)+h),(C1(i)+kC1_3*h));
The first line would change to this using our function handle from above:
kC1_1 = F_C1( t(i), C1(i), C2(i), R1(i), P(i) );
But what about the next line? It would have to change to something like this:
kC1_2 = F_C1( t(i)+0.5*h, C1(i)+0.5*h*kC1_1, C2(i)+0.5*h*kC2_1, R1(i)+0.5*h*kR1_1, P(i)+0.5*h*kP_1 );
where the C1, C2, R1, and P predictions are using their respective kC1_1 and kC2_1 and kR1_1 and kP_1 values. But some of these values haven't been computed yet for the current states, such as kP_1. Similar problems will exist in your other derivative calculations. Soooo ... you need to rearrange your code to compute all of the k_1 values first, then compute all of the k_2 values next, etc. E.g.,
kR1_1 = etc.
kR2_1 = etc.
kC1_1 = etc.
kC2_1 = etc.
kP_1 = etc.
kR1_2 = etc.
kR2_2 = etc.
kC1_2 = etc.
kC2_2 = etc.
kP_2 = etc.
:
etc.
Give a shot at making these changes and then come back to us with any further problems you might have.
Finally, at some point you will probably get tired of making lots of copy & paste of the RK4 code for each variable. One thing to look into is making your state a 5-element vector instead of using five different variables. That way you only need to write one set of RK4 code that will apply to this vector. This also makes it so that you can use the built-in functions such as ode45( ) to compare with.
6 个评论
Paulo Araujo Filho
2021-2-6
Paulo Araujo Filho
2021-2-6
编辑：Paulo Araujo Filho
2021-2-6
James Tursa
2021-2-7
编辑：James Tursa
2021-2-7
The problem is that in your cwd changed code, the F's are no longer function handles. You have inadvertently made them arrays that cannot take fractional inputs, hence the errors. You need to change this and make them function handles again. E.g., something like:
cwd.F_R1 = @(t,R1,C1)(R1.*(1-R1))-(xc*yc*((C1.*R1)./(R1+r0)));
And I think you can probably delete the cwd_rkfunc() function entirely since you will be initializing the function handles in the cwd_rkinitial() function.
Finally, it looks like you are missing an index on this line:
cwd.kR1_2 = cwd.F_R1( etc. )
which should be this instead
cwd.kR1_2(i) = cwd.F_R1( etc. )
And you should strongly consider pre-allocating all of your arrays inside your cwd_initial( ) function so that your code runs faster. You can use length(t) for this.
Paulo Araujo Filho
2021-2-9
James Tursa
2021-2-9
You keep changing the architecture from post to post, but I will try to keep up ...
This line
cwd = cwd_rkfunc(p,xc,xp,yc,yp,r0,c0);
wipes out EVERYTHING that was previously in cwd and replaces it with a cwd that contains only the function handles. You don't want that to happen. You want to create these function handles only once BEFORE the loop runs, and then use them inside the loop. So your code should look like this:
cwd = cwd_rkfunc(p,xc,xp,yc,yp,r0,c0); % do this first
% then allocate your states
cwd.R1 = zeros(1,length(t));
cwd.R2 = zeros(1,length(t));
cwd.C1 = zeros(1,length(t));
cwd.C2 = zeros(1,length(t));
cwd.P = zeros(1,length(t));
cwd.kR1_1 = zeros(1,length(t));
cwd.kR2_1 = zeros(1,length(t));
cwd.kC1_1 = zeros(1,length(t));
cwd.kC2_1 = zeros(1,length(t));
cwd.kP_1 = zeros(1,length(t));
cwd.kR1_2 = zeros(1,length(t));
cwd.kR2_2 = zeros(1,length(t));
cwd.kC1_2 = zeros(1,length(t));
cwd.kC2_2 = zeros(1,length(t));
cwd.kP_2 = zeros(1,length(t));
cwd.kR1_3 = zeros(1,length(t));
cwd.kR2_3 = zeros(1,length(t));
cwd.kC1_3 = zeros(1,length(t));
cwd.kC2_3 = zeros(1,length(t));
cwd.kP_3 = zeros(1,length(t));
cwd.kR1_4 = zeros(1,length(t));
cwd.kR2_4 = zeros(1,length(t));
cwd.kC1_4 = zeros(1,length(t));
cwd.kC2_4 = zeros(1,length(t));
cwd.kP_4 = zeros(1,length(t));
% then set your initial values
cwd.R1(1) = 1.0; % initial recourse 1
cwd.R2(1) = 1.0; % initial resourse 2
cwd.C1(1) = 1.0; % initial consumer 1
cwd.C2(1) = 1.0; % initial consumer 2
cwd.P(1) = 1.0; % initial predator
%% RK4
for i=1:(length(t)-1)
t(i+1)=t(i)+h;
% cwd = cwd_rkfunc(p,xc,xp,yc,yp,r0,c0); % move this line to above
cwd = cwd_rksteps(cwd,i,t,h);
cwd = cwd_rkint(cwd,i); % changed it to i
%cwd.R1(i)= R1; % delete this line
%cwd.R2(i)= R2; % delete this line
%cwd.C1(i)= C1; % delete this line
%cwd.C2(i)= C2; % delete this line
%cwd.P(i)= P; % delete this line
end
Then you will also need to change your cwd_rkint routine to use cwd properly:
function cwd = cwd_rkint(cwd,i)
cwd.R1(i+1) = cwd.R1(i) + (1/6)*(cwd.kR1_1(i)+2*cwd.kR1_2(i)+2*cwd.kR1_3(i)+cwd.kR1_4(i))*h;
cwd.R2(i+1) = cwd.R2(i) + (1/6)*(cwd.kR2_1(i)+2*cwd.kR2_2(i)+2*cwd.kR2_3(i)+cwd.kR2_4(i))*h;
cwd.C1(i+1) = cwd.C1(i) + (1/6)*(cwd.kC1_1(i)+2*cwd.kC1_2(i)+2*cwd.kC1_3(i)+cwd.kC1_4(i))*h;
cwd.C2(i+1) = cwd.C2(i) + (1/6)*(cwd.kC2_1(i)+2*cwd.kC2_2(i)+2*cwd.kC2_3(i)+cwd.kC2_4(i))*h;
cwd.P(i+1) = cwd.P(i) + (1/6)*(cwd.kP_1(i)+2*cwd.kP_2(i)+2*cwd.kP_3(i)+cwd.kP_4(i))*h;
end
Paulo Araujo Filho
2021-2-9
编辑：Paulo Araujo Filho
2021-2-9
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Introduction to Installation and Licensing 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
