Solve a differential equation system with 4th order Runge-Kutta method with three equations

10 Ansichten (letzte 30 Tage)
Hi, I'm trying to solve a Lotka-Volterra system with three equations via the 4th order Runge-Kutta method but I'm not being able to solve it.
This is the function I was using to try to solve it.
Any feedback would be helpfull
Thanks!
function [] = runge_kutta_sis(X,Y,Z,x0,y0,z0,a,b,h) %where X Y and Z are my ecuations, x0 y0 and z0 are the initial values. a is initial time and b is final time
t = a:h:b;
n = length(t);
x=[x0]; y=[y0]; z=[z0];
for i=1:n-1
j1 = h*A(x(i),y(i),z(i));
k1 = h*C(x(i),y(i),z(i));
l1 = h*L(x(i),y(i),z(i));
j2 = h*A(x(i)+j1/2,y(i)+k1/2,z(i)+l1/2,t(i)+h/2);
k2 = h*C(x(i)+j1/2,y(i)+k1/2,z(i)+l1/2,t(i)+h/2);
l2 = h*L(x(i)+j1/2,y(i)+k1/2,z(i)+l1/2,t(i)+h/2);
j3 = h*A(x(i)+j2/2,y(i)+k2/2,z(i)+l2/2,t(i)+h/2);
k3 = h*C(x(i)+j2/2,y(i)+k2/2,z(i)+l2/2,t(i)+h/2);
l3 = h*L(x(i)+j2/2,y(i)+k2/2,z(i)+l2/2,t(i)+h/2);
j4 = h*A(x(i)+j3,y(i)+k3,z(i)+l3,t(i)+h);
k4 = h*C(x(i)+j3,y(i)+k3,z(i)+l3,t(i)+h);
l4 = h*L(x(i)+j3,y(i)+k3,z(i)+l3,t(i)+h);
x(i+1)= x(i)+(h/6)*(j1+2*j2+2*j3*j4);
y(i+1) = y(i)+(h/6)*(k1+2*k2+2*k3*k4);
z(i+1) = z(i)+(h/6)*(l1+2*l2+2*l3*l4);
end
hold on
plot(t, x)
plot(t, y)
plot(t, z)

Akzeptierte Antwort

James Tursa
James Tursa am 24 Nov. 2020
Not sure if you actually use t in your derivative functions, but you are missing that input from your j1, k1, and l1 code:
j1 = h*A(x(i),y(i),z(i),t(i)); % added t(i)
k1 = h*C(x(i),y(i),z(i),t(i)); % added t(i)
l1 = h*L(x(i),y(i),z(i),t(i)); % added t(i)
Note that all of your j, k, l variables already have the h factor applied to them. So this code should not again multiply by h:
x(i+1) = x(i)+(1/6)*(j1+2*j2+2*j3*j4); % changed h to 1
y(i+1) = y(i)+(1/6)*(k1+2*k2+2*k3*k4); % changed h to 1
z(i+1) = z(i)+(1/6)*(l1+2*l2+2*l3*l4); % changed h to 1

Weitere Antworten (0)

Kategorien

Mehr zu Programming finden Sie in Help Center und File Exchange

Produkte

Community Treasure Hunt

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

Start Hunting!

Translated by