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

3 次查看（过去 30 天）
Sajani
Sajani 2025-4-6
评论： Sajani 2025-4-7
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
  6 个评论
显示 4更早的评论
Torsten
Torsten 2025-4-6
编辑：Torsten 2025-4-6
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;

请先登录，再进行评论。

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

采纳的回答

VBBV
VBBV 2025-4-6
编辑：VBBV 2025-4-6
@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

更多回答（0 个）

类别

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

标签

Community Treasure Hunt

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

Start Hunting!

Translated by