Problem with time steps in Runge-Kutte method

I am solving 6 ODEs with Runge-Kutte method. But I get different results with different time steps. Please suggest how to deal with this. My code is following:

clear;
close all;
clc;

nu=0.89*10^-6 ;
g=9.81;
ae=10 ;     % aspect ratio
ds=1*10^-3 ;  % equivalent dia
ecc=((ae.^2)-1)./((ae.^2)+1) ;   % eccentricity
A=0.2;   % wave amplitude
k=0.33/A; % k*A <= 0.33 for Linear wave theory to hold correct
w1=sqrt(g*k); % angular frequency
T=2*pi/w1 ;    % time period
B=1.1 ;
theta_ini= 0  ;  % azimuthal angle matrix
phi_ini= pi/4 ;    % polar angle matrix
zh= 3 ;  
vp= 0.01*w1*A ;    % Particle intrinsic velocity (m/s)

% Time span for evolution of dynamics
t0 = 0; %seconds
tf =20; %seconds
tspan=[t0 tf] ;
tau_s=5 ;
h_mat=[0.1 0.05 0.01 0.005];

for jh=1:4
    h=h_mat(jh);
    % resistance tensor in body axes
    K11d=((4/3).*(ae.^(-1/3)).*(1-ae.^2))./(ae-((2*(ae.^2)-1).*log(ae+sqrt(ae.^2-1))./sqrt(ae.^2-1))) ;
    K22d=((8/3)*(ae.^(-1/3)).*((ae.^2)-1))./(ae+((2*(ae.^2)-3).*log(ae+sqrt(ae.^2-1))./sqrt(ae.^2-1))) ;
    K33d=K22d ;
    
    K33_ini=K11d*cos(phi_ini)*cos(phi_ini)+K22d*sin(phi_ini)*sin(phi_ini) ; % K33 @t=0
    
    x0=0;
    y0=0;
    z0=0;
    % Initial field velocity @t=0, x=0, z=0
    u0x=w1*A;  
    u0y=0;
    u0z=0;
    
    % Initial particle orientation
    p0x = sin(phi_ini).*cos(theta_ini);
    p0y = sin(phi_ini).*sin(theta_ini);
    p0z = cos(phi_ini);
    
    % Initial particle velocity
    v0x=u0x+vp*p0x ;
    v0y=u0y+vp*p0y ;
    v0z=-((B-1)*g*ds^2)/(18*nu*K33_ini)+u0z+vp*p0z ;
    
    nu=0.89*10^-6 ;
    c1=18*nu/(B*ds^2) ;
    g=9.81;
    
    % mass tensor in body axes
    alpha0=(1./((ae.^2)-1).^(3/2)).*(2*sqrt((ae.^2)-1)+ae.*log((ae-sqrt(ae.^2-1))./(ae+sqrt(ae.^2-1)))) ;
    beta0=(0.5./((ae.^2)-1).^(3/2)).*(2.*(ae.^2).*sqrt((ae.^2)-1)+ae.*log((ae-sqrt(ae.^2-1))./(ae+sqrt(ae.^2-1)))) ;
    gamma0=beta0 ;
    C11d=alpha0./(2-alpha0) ;
    C22d=beta0./(2-beta0) ;
    C33d=gamma0./(2-gamma0) ;
    
    t = tspan(1):h:tspan(2);
    
    x=zeros(length(t),1);
    vx=zeros(length(t),1);
    px=zeros(length(t),1);
    y=zeros(length(t),1);
    vy=zeros(length(t),1);
    py=zeros(length(t),1);
    z=zeros(length(t),1);
    vz=zeros(length(t),1);
    pz=zeros(length(t),1);
    dx=zeros(length(t),1);
    dvx=zeros(length(t),1);
    
    x(1) = x0 ;
    vx(1) = v0x ;
    y(1) = y0 ;
    vy(1) = v0y ;
    z(1) = z0 ;
    vz(1) = v0z ;
    px(1) = p0x ;
    py(1) = p0y ;
    pz(1) = p0z ;
    
    Fx = @(t,x,vx,y,vy,z,vz,px,py,pz) [vx; (1./(1+(1/B).*(C11d.*sin(acos(pz)).*sin(acos(pz))+C22d.*cos(acos(pz)).*cos(acos(pz))))).*((1/B).*((w1^2)*A*exp(k*z)*sin(k*x-w1*t)+...
        w1*A*exp(k*z)*cos(k*x-w1*t).*(-k*w1*A*exp(k*z)*sin(k*x-w1*t))+w1*A*exp(k*z)*sin(k*x-w1*t).*(k*w1*A*exp(k*z)*cos(k*x-w1*t)))+...
        (1/B)*((C11d.*sin(acos(pz)).*sin(acos(pz))+C22d.*cos(acos(pz)).*cos(acos(pz)))*(w1^2)*A*exp(k*z)*sin(k*x-w1*t))-...
        c1*((K11d.*sin(acos(pz)).*sin(acos(pz))+K22d.*cos(acos(pz)).*cos(acos(pz))).*(vx-w1*A*exp(k*z)*cos(k*x-w1*t)))); vy; (1./(1+(1/B).*C22d)).*(-c1.*(K22d.*vy));vz; (1./(1+(1/B).*(C11d.*cos(acos(pz)).*cos(acos(pz))+C22d.*sin(acos(pz)).*sin(acos(pz))))).*((1/B).*(-(w1^2)*A*exp(k*z)*cos(k*x-w1*t)+...
        w1*A*exp(k*z)*cos(k*x-w1*t).*(k*w1*A*exp(k*z)*cos(k*x-w1*t))+w1*A*exp(k*z)*sin(k*x-w1*t).*(k*w1*A*exp(k*z)*sin(k*x-w1*t)))+...
        (1/B)*((C11d.*cos(acos(pz)).*cos(acos(pz))+C22d.*sin(acos(pz)).*sin(acos(pz)))*(-(w1^2)*A*exp(k*z)*cos(k*x-w1*t)))-...
        c1.*((K11d.*cos(acos(pz)).*cos(acos(pz))+K22d.*sin(acos(pz)).*sin(acos(pz))).*(vz-w1*A*exp(k*z)*sin(k*x-w1*t)))-(1-1/B).*g); -px.*pz/tau_s+ecc.*((-k*w1*A*exp(k*z)*sin(k*x-w1*t)).*px+(k*w1*A*exp(k*z)*cos(k*x-w1*t)).*pz)-ecc.*px.*((px.*(-k*w1*A*exp(k*z)*sin(k*x-w1*t)).*px+...
        pz.*(k*w1*A*exp(k*z)*cos(k*x-w1*t)).*px)+(px.*(k*w1*A*exp(k*z)*cos(k*x-w1*t)).*pz+pz.*(k*w1*A*exp(k*z)*sin(k*x-w1*t)).*pz));-py.*pz/tau_s-ecc.*py.*((px.*(-k*w1*A*exp(k*z)*sin(k*x-w1*t)).*px+pz.*(k*w1*A*exp(k*z)*cos(k*x-w1*t)).*px)+...
        (px.*(k*w1*A*exp(k*z)*cos(k*x-w1*t)).*pz+pz.*(k*w1*A*exp(k*z)*sin(k*x-w1*t)).*pz));(1-pz.*pz)/tau_s+ecc.*((k*w1*A*exp(k*z)*cos(k*x-w1*t)).*px+(k*w1*A*exp(k*z)*sin(k*x-w1*t)).*pz)-ecc.*pz.*((px.*(-k*w1*A*exp(k*z)*sin(k*x-w1*t)).*px+...
        pz.*(k*w1*A*exp(k*z)*cos(k*x-w1*t)).*px)+(px.*(k*w1*A*exp(k*z)*cos(k*x-w1*t)).*pz+pz.*(k*w1*A*exp(k*z)*sin(k*x-w1*t)).*pz))];
    
    for i=1: length(t)-1
        
        K1 = Fx(t(i),x(i),vx(i),y(i),vy(i),z(i),vz(i),px(i),py(i),pz(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), px(i)+0.5*h*K1(7), py(i)+0.5*h*K1(8), pz(i)+0.5*h*K1(9));  
        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), px(i)+0.5*h*K2(7), py(i)+0.5*h*K2(8), pz(i)+0.5*h*K2(9));  
        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), px(i)+h*K3(7), py(i)+h*K3(8), pz(i)+h*K3(9));
        
        px(i+1) = px(i) + (1/6)*(K1(7)+2*K2(7)+2*K3(7)+K4(7))*h;  
        py(i+1) = py(i) + (1/6)*(K1(8)+2*K2(8)+2*K3(8)+K4(8))*h;  
        pz(i+1) = pz(i) + (1/6)*(K1(9)+2*K2(9)+2*K3(9)+K4(9))*h;  
        vx(i+1) = vx(i) + (1/6)*(K1(2)+2*K2(2)+2*K3(2)+K4(2))*h;  
        vy(i+1) = vy(i) + (1/6)*(K1(4)+2*K2(4)+2*K3(4)+K4(4))*h;
        vz(i+1) = vz(i) + (1/6)*(K1(6)+2*K2(6)+2*K3(6)+K4(6))*h;
        
        vx(i+1) = vx(i+1) + vp*px(i+1) ;    
        vy(i+1) = vy(i+1) + vp*py(i+1) ;
        vz(i+1) = vz(i+1) + vp*pz(i+1) ;
        x(i+1) = x(i)+vx(i+1)*h;  
        y(i+1) = y(i)+vy(i+1)*h;  
        z(i+1) = z(i)+vz(i+1)*h;   
    end
    
    % % plotting of results
    plot(k*x, k*z, '-r','LineWidth',2)
    xlabel('$kx$','FontSize',20,'FontWeight','bold', 'Interpreter','latex');
    ylabel('$kz$','FontSize',20,'FontWeight','bold', 'Interpreter','latex');
    set(gca,'FontSize',15);
    hold on
end

Answer 1

[James Tursa]

"... I get different results with different time steps ..."

In general, that is the expected result. I suppose the real question is are any of your results valid? To answer that, I suggest you rewrite your derivative equations in a vector format so that it can be fed to one of the MATLAB solvers such as ode45( ) for comparison. E.g., create this function handle to use with ode45( ):

F = @(t,y) Fx(t,y(1),y(2),y(3),y(4),y(5),y(6),y(7),y(8),y(9))

[AkB]

Thanks. I will use a vector form. Please suggest me one thing. My function looks combersome with many eqns like this:

Fx = @(t,x,vx,y,vy) [vx; ((w1^2)*A*exp(k*z)*sin(k*x-w1*t)+w1*A*exp(k*z)*cos(k*x-w1*t).*(-k*w1*A*exp(k*z)*sin(k*x-w1*t))+w1*A*exp(k*z)*sin(k*x-w1*t).*(k*w1*A*exp(k*z)*cos(k*x-w1*t)))+(1/B)*((C11d.*sin(acos(pz)).*sin(acos(pz))+C22d.*cos(acos(pz)).*cos(acos(pz)))*(w1^2)*A*exp(k*z)*sin(k*x-w1*t)); vy; (1./(1+(1/B).*C22d)).*(-c1.*(K22d.*vy))];

Instead, I would like to re-write it as below:

dvxdt=((w1^2)*A*exp(k*z)*sin(k*x-w1*t)+w1*A*exp(k*z)*cos(k*x-w1*t).*(-k*w1*A*exp(k*z)*sin(k*x-w1*t))+w1*A*exp(k*z)*sin(k*x-w1*t).*(k*w1*A*exp(k*z)*cos(k*x-w1*t)))+(1/B)*((C11d.*sin(acos(pz)).*sin(acos(pz))+C22d.*cos(acos(pz)).*cos(acos(pz)))*(w1^2)*A*exp(k*z)*sin(k*x-w1*t));

dvydt=(1./(1+(1/B).*C22d)).*(-c1.*(K22d.*vy));

Fx = @(t,x,vx,y,vy) [vx; dvxdt; vy; dvydt];

However, I get errors saying my variables are undefined. Please suggest.

[Torsten]

Do it this way:

https://de.mathworks.com/matlabcentral/answers/1941229-help-on-solving-nonlinear-coupled-equations-with-runge-kutta-approach?s_tid=srchtitle

[Ban]

My eqns are pretty messy. So in terms of marix it looks very confusing. I am trying to write somrthing like this:

fx = @(t,x,vx,y,vy,z,vz) C11d.*sin(acos(pz)) ;

fy = @(t,x,vx,y,vy,z,vz) C22d.*(K22d.*vy)) ;

fz = @(t,x,vx,y,vy,z,vz) C11d.*cos(acos(pz)) ;

Fx = @(t,x,vx,y,vy,z,vz) [vx; fx; vy; fy; vz; fz; 0; 0; 0];

It gives this error:

The following error occurred converting from double to function_handle:

Too many output arguments.

[Torsten]

Write a function F in which you compute the time derivatives of your solution variables and call this function by the Runge-Kutta solver. Just look at the code I linked to. Your equations don't need to have a matrix structure for this:

function dZ = F(t,Z)

x = Z(1);

vx = Z(2);

y = Z(3);

vy = Z(4);

z = Z(5);

vz = Z(6);

px = Z(7);

py = Z(8);

pz = Z(9);

% Compute time derivatives dx, dvx, dy, dvy, dz, dvz, dpx, dpy dpz without using function handles

...

% Return the vector of time derivatives [dx dvx dy dvy dz dvz dpx dpy dpz] to runge-kutta code

dZ = [dx dvx dy dvy dz dvz dpx dpy dpz].'

end

But I still don't understand why you don't use ODE45. Do you think your Runge-Kutta solver will outperform MATLAB software ?

Answer 2

[Torsten]

        vx(i+1) = vx(i) + (1/6)*(K1(2)+2*K2(2)+2*K3(2)+K4(2))*h;  
        vy(i+1) = vy(i) + (1/6)*(K1(4)+2*K2(4)+2*K3(4)+K4(4))*h;
        vz(i+1) = vz(i) + (1/6)*(K1(6)+2*K2(6)+2*K3(6)+K4(6))*h;
        
        vx(i+1) = vx(i+1) + vp*px(i+1) ;    
        vy(i+1) = vy(i+1) + vp*py(i+1) ;
        vz(i+1) = vz(i+1) + vp*pz(i+1) ;
        x(i+1) = x(i)+vx(i+1)*h;  
        y(i+1) = y(i)+vy(i+1)*h;  
        z(i+1) = z(i)+vz(i+1)*h;      

This is not RK4 for vx, vy, vz, x, y and z.

If you use the usual update, you get the same results independent of h:

        vx(i+1) = vx(i) + (1/6)*(K1(2)+2*K2(2)+2*K3(2)+K4(2))*h;  
        vy(i+1) = vy(i) + (1/6)*(K1(4)+2*K2(4)+2*K3(4)+K4(4))*h;
        vz(i+1) = vz(i) + (1/6)*(K1(6)+2*K2(6)+2*K3(6)+K4(6))*h;
        x(i+1) = x(i) + (1/6)*(K1(1)+2*K2(1)+2*K3(1)+K4(1))*h;  
        y(i+1) = y(i) + (1/6)*(K1(3)+2*K2(3)+2*K3(3)+K4(3))*h;
        z(i+1) = z(i) + (1/6)*(K1(5)+2*K2(5)+2*K3(5)+K4(5))*h;