Problem using 4th order Runge-Kutta method

I am solving 3 equations using Runge-Kutta method. But the solution doesn't match with exact solution. Please suggest where am I making mistake.

clear;

close all;

clc;

Fx = @(t,x,vx,y,vy,z,vz) [vx; 3*t; vy; t; vz; 7*t] ;

dt=0.1 ;

t=linspace(0, 10, 100);

h=dt;

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

vx=zeros(length(t),1);

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

vy=zeros(length(t),1);

z=zeros(length(t),1);

vz=zeros(length(t),1);

x(1) = 0 ;

vx(1) = 0 ;

y(1) = 0 ;

vy(1) = 0 ;

z(1) = 0 ;

vz(1) = 0 ;

for i=1: length(t)-1

K1 = Fx(t(i),x(i),vx(i),y(i),vy(i),z(i),vz(i));

K2 = Fx(t(i)+0.5*h, x(i)+0.5*h*K1(1), vx(i)+0.5*h*K1(2), y(i)+0.5*h*K1(3), vy(i)+0.5*h*K1(4), z(i)+0.5*h*K1(5), vz(i)+0.5*h*K1(6));

K3 = Fx(t(i)+0.5*h, x(i)+0.5*h*K2(1), vx(i)+0.5*h*K2(2), y(i)+0.5*h*K2(3), vy(i)+0.5*h*K2(4), z(i)+0.5*h*K2(5), vz(i)+0.5*h*K2(6));

K4 = Fx(t(i)+h, x(i)+h*K3(1), vx(i)+h*K3(2), y(i)+h*K3(3), vy(i)+h*K3(4), z(i)+h*K3(5), vz(i)+h*K3(6));

x(i+1) = x(i) + (1/6)*(K1(1)+2*K2(1)+2*K3(1)+2*K4(1))*h;

vx(i+1) = vx(i) + (1/6)*(K1(2)+2*K2(2)+2*K3(2)+2*K4(2))*h;

y(i+1) = y(i) + (1/6)*(K1(3)+2*K2(3)+2*K3(3)+2*K4(3))*h;

vy(i+1) = vy(i) + (1/6)*(K1(4)+2*K2(4)+2*K3(4)+2*K4(4))*h;

z(i+1) = z(i) + (1/6)*(K1(5)+2*K2(5)+2*K3(5)+2*K4(5))*h;

vz(i+1) = vz(i) + (1/6)*(K1(6)+2*K2(6)+2*K3(6)+2*K4(6))*h;

end

plot(t,x,t,y,t,z,'LineWidth',2)

hold on

plot(t,t.^3/2, t,(1/6)*t.^3,t, (7/6)*t.^3,'LineWidth',2)

xlabel('$t$','FontSize',20,'FontWeight','bold', 'Interpreter','latex');

ylabel('$x, y, z$', 'FontSize',20,'FontWeight','bold', 'Interpreter','latex');

set(gca,'FontSize',15);

legend('x','y','z','x_exact','y_exact','z_exact')

[James Tursa]

First thing I suggest you do is rewrite all of your code in a vectorized fashion. E.g., write your function handle as

Fx = @(t,y) [function of t and y(1), y(2), etc.] ;

Where y is a 6-element vector. I would suggest you reorder your terms as follows:

y(1) = x

y(2) = y

y(3) = z

y(4) = vx

y(5) = vy

y(6) = vz

Then your function handle becomes:

Fx = @(t,y) [y(4:6); 3*t; t; 7*t]

This will greatly simplify your downstream code. It will also allow you to directly use this function handle in a call to ode45( ) so you can compare your RK4 results with the MATLAB ode function results.

Answer 1

[James Tursa]

You have a factor of 2*K4 that shouldn't be there. It should just be K4. E.g.,

Fx = @(t,x,vx,y,vy,z,vz) [vx; 3*t; vy; t; vz; 7*t] ; 
dt=0.01 ;
t=linspace(0, 10, ceil(10/dt));
h=dt;
x=zeros(length(t),1);
vx=zeros(length(t),1);
y=zeros(length(t),1);
vy=zeros(length(t),1);
z=zeros(length(t),1);
vz=zeros(length(t),1);
x(1) = 0 ;
vx(1) = 0 ;
y(1) = 0 ;
vy(1) = 0 ;
z(1) = 0 ;
vz(1) = 0 ;
for i=1: length(t)-1                              
    
    K1 = Fx(t(i),x(i),vx(i),y(i),vy(i),z(i),vz(i));    
    K2 = Fx(t(i)+0.5*h, x(i)+0.5*h*K1(1), vx(i)+0.5*h*K1(2), y(i)+0.5*h*K1(3), vy(i)+0.5*h*K1(4), z(i)+0.5*h*K1(5), vz(i)+0.5*h*K1(6));   
    K3 = Fx(t(i)+0.5*h, x(i)+0.5*h*K2(1), vx(i)+0.5*h*K2(2), y(i)+0.5*h*K2(3), vy(i)+0.5*h*K2(4), z(i)+0.5*h*K2(5), vz(i)+0.5*h*K2(6));   
    K4 = Fx(t(i)+h, x(i)+h*K3(1), vx(i)+h*K3(2), y(i)+h*K3(3), vy(i)+h*K3(4), z(i)+h*K3(5), vz(i)+h*K3(6));
    
    x(i+1) = x(i) + (1/6)*(K1(1)+2*K2(1)+2*K3(1)+K4(1))*h;  
    vx(i+1) = vx(i) + (1/6)*(K1(2)+2*K2(2)+2*K3(2)+K4(2))*h;  
    
    y(i+1) = y(i) + (1/6)*(K1(3)+2*K2(3)+2*K3(3)+K4(3))*h;  
    vy(i+1) = vy(i) + (1/6)*(K1(4)+2*K2(4)+2*K3(4)+K4(4))*h; 
    
    z(i+1) = z(i) + (1/6)*(K1(5)+2*K2(5)+2*K3(5)+K4(5))*h;  
    vz(i+1) = vz(i) + (1/6)*(K1(6)+2*K2(6)+2*K3(6)+K4(6))*h;   
end
plot(t,x,t,y,t,z,'LineWidth',2) 
hold on
plot(t,t.^3/2, t,(1/6)*t.^3,t, (7/6)*t.^3,'LineWidth',2)
xlabel('$t$','FontSize',20,'FontWeight','bold', 'Interpreter','latex');
ylabel('$x, y, z$', 'FontSize',20,'FontWeight','bold', 'Interpreter','latex');
set(gca,'FontSize',15);
legend('x','y','z','x_exact','y_exact','z_exact')