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

13 vues (au cours des 30 derniers jours)
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)

Réponse acceptée

James Tursa
James Tursa le 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

Plus de réponses (0)

Catégories

En savoir plus sur Programming dans Help Center et File Exchange

Produits

Community Treasure Hunt

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

Start Hunting!

Translated by