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

17 次查看（过去 30 天）
Miguel Bublitz San Roman
回答： James Tursa 2020-11-24
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)

请先登录，再回答此问题。

采纳的回答

James Tursa
James Tursa 2020-11-24
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

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Programming 的更多信息

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by