Why did I get incorrect values when solving this equation using the fourth order Runge-Kutta Method?

This is my equation,

And this is my code,

B = 30;  % width of mouth, Ba + Bb = B
Y = 17600000;  % surface area of the lagoon
a = 0.1;  % tidal amplitude
g = 9.81;  % gravitational accelaration
Cd = 0.03;  % drag coefficient
L = 500;  % length of the channel, La = Lb = L
H = 3;  % mean depths of the mouth, Ha = Hb = H
Qf = 30;  % net fresh water supply
T = 12.42*3600;  % tidal period (the time it takes for a full tidal cycle)

etao_data = [0.020391;0.031391;0.043391]; % for 10min
etao_time = 0:300:10*60;
dt= 1;
t_s = 0:dt:10*60;
etao_star= interp1(etao_time,etao_data,t_s,'spline');
P = ((g * B^2 * T^2 * H^3) / (Y^2 * L * Cd * a))^5;
S = T * Qf / (a * Y);

n = length(t_s);
etal_star_17 = zeros(1, n);  % Initialize eta_l_star
etal_star_17(1) = 0;  % Initial condition

for i = 1:n-1
    eta_o = etao_star(i);
    eta_l = etal_star_17(i);
    eta_m = (eta_o +eta_l)/2;
    f = @(t, eta) P* ((1 + (a * eta_m / H)).^(1.5)) ./ (sqrt(1 + ((H/ (2 * Cd * L)) * (1+ (a * eta_m / H))))).* ((eta_o - eta_l) ./ sqrt(abs(eta_o - eta_l)))+ S;

    k1 = dt* f(t_s(i), eta_l);
    k2 = dt* f(t_s(i) + dt/2, eta_l * k1 / 2);
    k3 = dt* f(t_s(i) + dt/2, eta_l * k2 / 2);
    k4 = dt* f(t_s(i)+ dt , eta_l + k3);

    etal_star_17(i+1) = eta_l + (k1 + 2*k2 + 2*k3 + k4)/6;
end

[VBBV]

Try with those changes mentioned

[Torsten]

Your Runge-Kutta method is wrong:

    k2 = dt*f(t_s(i) + dt/2, eta_l + k1 / 2);
    k3 = dt*f(t_s(i) + dt/2, eta_l + k2 / 2); 

instead of

    k2 = dt* f(t_s(i) + dt/2, eta_l * k1 / 2);
    k3 = dt* f(t_s(i) + dt/2, eta_l * k2 / 2);

[VBBV]

@Sajani as @Torsten mentioned that is also one more change,

[Sajani]

@VBBV and @Torsten i found that,

k2 = dt* f(t_s(i) + dt/2, eta_l + dt * k1 / 2);
k3 = dt* f(t_s(i) + dt/2, eta_l + dt * k2 / 2);
k4 = dt* f(t_s(i)+ dt , eta_l + dt * k3);

So what is correct one?

[Torsten]

Either use

    k1 = f(t_s(i), eta_l);
    k2 = f(t_s(i) + dt/2, eta_l + dt * k1 / 2);
    k3 = f(t_s(i) + dt/2, eta_l + dt * k2 / 2);
    k4 = f(t_s(i)+ dt , eta_l + dt * k3);
    etal_star_17(i+1) = eta_l + dt * (k1 + 2*k2 + 2*k3 + k4)/6;  

or

    k1 = dt* f(t_s(i), eta_l);
    k2 = dt* f(t_s(i) + dt/2, eta_l + k1 / 2);
    k3 = dt* f(t_s(i) + dt/2, eta_l + k2 / 2);
    k4 = dt* f(t_s(i)+ dt , eta_l + k3);
    etal_star_17(i+1) = eta_l + (k1 + 2*k2 + 2*k3 + k4)/6;        

The version you found is wrong.

Further, your function f = @(t, eta) does not depend on the input variable "eta". Thus you will have to use

    f = @(t, eta) P* ((1 + (a * eta_m / H)).^(1.5)) ./ (sqrt(1 + ((H/ (2 * Cd * L)) * (1+ (a * eta_m / H))))).* ((eta_o - eta) ./ sqrt(abs(eta_o - eta)))+ S;    

instead of

    f = @(t, eta) P* ((1 + (a * eta_m / H)).^(1.5)) ./ (sqrt(1 + ((H/ (2 * Cd * L)) * (1+ (a * eta_m / H))))).* ((eta_o - eta_l) ./ sqrt(abs(eta_o - eta_l)))+ S;    

or something even more complicated since eta_m also depends on eta_l:

    f = @(t, eta) P* ((1 + (a * (eta_o + eta)/2 / H)).^(1.5)) ./ (sqrt(1 + ((H/ (2 * Cd * L)) * (1+ (a * (eta_o + eta)/2 / H))))).* ((eta_o - eta) ./ sqrt(abs(eta_o - eta)))+ S;

[Sajani]

@Torsten Thank you, I got it.

Answer 1

[VBBV]

@Sajani there are some missing values for the equation. Also the exponent in equation is 0.5

B = 30;  % width of mouth, Ba + Bb = B
Y = 17600000;  % surface area of the lagoon
a = 0.1;  % tidal amplitude
g = 9.81;  % gravitational accelaration
Cd = 0.03;  % drag coefficient
L = 500;  % length of the channel, La = Lb = L
H = 3;  % mean depths of the mouth, Ha = Hb = H
Qf = 30;  % net fresh water supply
T = 12.42*3600;  % tidal period (the time it takes for a full tidal cycle)

etao_data = [0.020391;0.031391;0.043391]; % for 10min
etao_time = 0:300:10*60;
dt= 1;
t_s = 0:dt:10*60;
etao_star= interp1(etao_time,etao_data,t_s,'spline');
P = ((g * B^2 * T^2 * H^3) / (Y^2 * L * Cd * a))^0.5;
                                               % ---vv

A = 3; % a value 
S = T * Qf / (a * A);

n = length(t_s);
etal_star_17 = zeros(1, n);  % Initialize eta_l_star
etal_star_17(1) = 0;  % Initial condition

for i = 1:n-1
    eta_o = etao_star(i);
    eta_l = etal_star_17(i);
    eta_m = (eta_o +eta_l)/2;
    f = @(t, eta) P* ((1 + (a * eta_m / H)).^(1.5)) ./ (sqrt(1 + ((H/ (2 * Cd * L)) * (1+ (a * eta_m / H))))).* ((eta_o - eta_l) ./ sqrt(abs(eta_o - eta_l)))+ S;
                                           
    k1 = dt* f(t_s(i), eta_l);
    k2 = dt* f(t_s(i) + dt/2, eta_l * k1 / 2);
    k3 = dt* f(t_s(i) + dt/2, eta_l * k2 / 2);
    k4 = dt* f(t_s(i)+ dt , eta_l + k3);

    etal_star_17(i+1) = eta_l + (k1 + 2*k2 + 2*k3 + k4)/6;
end
%Y = 17600000;  % surface area of the lagoona = 0.1;  % tidal amplitudeg = 9.81;  % gravitational accelarationCd = 0.03;  % drag coefficientL = 500;  % length of the channel, La = Lb = LH = 3;  % mean depths of the mouth, Ha = Hb = HQf = 30;  % net fresh water supplyT = 12.42*3600;  % tidal period (the time it takes for a full tidal cycle)etao_data = [0.020391;0.031391;0.043391]; % for 10minetao_time = 0:300:10*60;dt= 1;t_s = 0:dt:10*60;etao_star= interp1(etao_time,etao_data,t_s,'spline');P = ((g * B^2 * T^2 * H^3) / (Y^2 * L * Cd * a))^5;S = T * Qf / (a * Y);n = length(t_s);etal_star_17 = zeros(1, n);  % Initialize eta_l_staretal_star_17(1) = 0;  % Initial conditionfor i = 1:n-1  eta_o = etao_star(i);  eta_l = etal_star_17(i);  eta_m = (eta_o +eta_l)/2;  f = @(t, eta) P* ((1 + (a * eta_m / H)).^(1.5)) ./ (sqrt(1 + ((H/ (2 * Cd * L)) * (1+ (a * eta_m / H))))).* ((eta_o - eta_l) ./ sqrt(abs(eta_o - eta_l)))+ S;  k1 = dt* f(t_s(i), eta_l);  k2 = dt* f(t_s(i) + dt/2, eta_l * k1 / 2);  k3 = dt* f(t_s(i) + dt/2, eta_l * k2 / 2);  k4 = dt* f(t_s(i)+ dt , eta_l + k3);  etal_star_17(i+1) = eta_l + (k1 + 2*k2 + 2*k3 + k4)/6;end

[Sajani]

@VBBV Thank you. I got it. Then what about the k2, k3 & k4?