Manual Runge-Kutta for system of two ODEs.

I am struggling to obtain the correct graph for the system of ODEs as follows:

x'=-y+6x, y'=-y+4x, between t=0,0.7

I can obtain the correct graph using Euler's method, as seen here:

But cannot do the same for a manual Runge Kutta method. And I don't want to use the integrated ode45 functions if I don't have to. What am I doing wrong? My code is below:

clear,clc

h = 0.1
t_beg = 0
t_end = 0.7
x_initial= 0.5
y_initial= -0.5


F_tx=@(x,y)(-x+6*y);
F_ty=@(x,y)(-y+4*x);

t=t_beg:h:t_end;

x=zeros(1,length(t));
x(1)=x_initial;

y=zeros(1,length(t));
y(1)=y_initial;

for i=1:(length(t)-1)
    
    kx1 = F_tx(t(i),x(i));
    kx2 = F_tx(t(i)+0.5*h,x(i)+0.5*h*kx1);
    kx3 = F_tx((t(i)+0.5*h),(x(i)+0.5*h*kx2));
    kx4 = F_tx((t(i)+h),(x(i)+kx3*h));
    
    x(i+1) = x(i) + (1/6)*(kx1+2*kx2+2*kx3+kx4)*h;
    
    ky1 = F_ty(t(i),y(i));
    ky2 = F_ty(t(i)+0.5*h,y(i)+0.5*h*ky1);
    ky3 = F_ty((t(i)+0.5*h),(y(i)+0.5*h*ky2));
    ky4 = F_ty((t(i)+h),(y(i)+ky3*h));
    
    y(i+1) = y(i) + (1/6)*(ky1+2*ky2+2*ky3+ky4)*h;
    
end

% plot(x,y)

figure(1)
plot(t,y)
hold on
plot(t,x)

Answer 1

[Alan Stevens]

You need to change the order within the loop to

for i=1:(length(t)-1)
    
    kx1 = F_tx(x(i),y(i));
    ky1 = F_ty(x(i),y(i));
    kx2 = F_tx(x(i)+0.5*h*kx1,y(i)+0.5*h*ky1);
    ky2 = F_ty(x(i)+0.5*h*kx1,y(i)+0.5*h*ky1);
    kx3 = F_tx((x(i)+0.5*h*kx2),y(i)+0.5*h*ky2);
    ky3 = F_ty((x(i)+0.5*h*kx2),(y(i)+0.5*h*ky2));
    kx4 = F_tx((x(i)+kx3*h),y(i)+ky3*h);
    ky4 = F_ty((x(i)+kx3*h),(y(i)+ky3*h));
    
    x(i+1) = x(i) + (1/6)*(kx1+2*kx2+2*kx3+kx4)*h; 
    y(i+1) = y(i) + (1/6)*(ky1+2*ky2+2*ky3+ky4)*h;
    
end

and note that the first argument is x not t.

[FPixelz]

How would I apply this to a set of three odes? Say, the rossler attactor where there is now x',y' and z'? Would kx1 need to look something like kx1 = F_tx(x(i),y(i),z(i)); even though there is not necessarily an x component in F_tx?

[Alan Stevens]

Presumably you mean even if there is not a z-component in F_tx. The answer is no, you only need include the necessary arguments.

[James Tursa]

For three variables x, y, z you still need to respect the order of the k evaluations. Do kx1, ky1, kz1 first. Then do kx2, ky2, kz2. Etc.

Or to get this same effect use the vector approach that Jan has posted.

[FPixelz]

Thanks for your help!

Answer 2

[Jan]

The diagram looks, yoike you are integrating:

@(t, y) [-y(1)+6*y(2); -y(2)+4*y(1)]

This code produces an equivalent output:

t0 = 0;
tF = 0.7;
x0 = 0.5;
y0 = -0.5;

[t, y] = ode45(@(t,y) [-y(1)+6*y(2); -y(2)+4*y(1)], ...
               [t0, tF], [x0, y0]);
figure;
plot(t,y);

But your function to be integrated is something else:

F_tx = @(x,y) (-x + 6 * y);
F_ty = @(x,y) (-y + 4 * x);

Here the function depends on the 1st input, which is t in my code. This is a confusion of "x/y" versus "t/y", whereby your "y" consists of the components x and y.

[EDITED] A working solution:

F = @(t, y) [-y(1) + 6 * y(2); ...
             -y(2) + 4 * y(1)];

y       = zeros(2, length(t));
y(:, 1) = [x_initial; y_initial];
for i=1:(length(t)-1)    
    kx1 = F(t(i),           y(:, i));
    kx2 = F(t(i) + 0.5 * h, y(:, i) + 0.5 * h * kx1);
    kx3 = F(t(i) + 0.5 * h, y(:, i) + 0.5 * h * kx2);
    kx4 = F(t(i) + h,       y(:, i) + kx3 * h);
    
    y(:, i+1) = y(:, i) + (kx1 + 2 * kx2 + 2 * kx3 + kx4) * h / 6;
end
figure()
plot(t, y)

By the way, compare the readability of the code, which contains spaces around the operators.

[FPixelz]

Thanks for the help, this has proved very useful to me.

Answer 3

[khalida]

I am struggling to obtain the correct algorithm for the system of ODEs as follows:

y'=z, z'=-((1+5x)/(2x(x+1))z-y^3+5x+11x^2+0.296x^9+0.666x^10+0.5x^11+0.125x^12), between x=0,0.7

y(0)=0, y'(0)=z(0)=0

i need the matlab code for this system of odes

[Jan]

Please do not append a new question in the section for answers of another question. Open a new thread instead and delete this message here. Post, what you have tried so far.