Unstable 6 equation ode
Hello all, I am working on solving a 6 equation ODE to model the shape of a free surface in a channel. The paper that I am referencing only uses 5 euqations but i have no idea how they are solving for one of the variables so I added the variable and equation to the ODE. I have tried almost solver avalible in matlab and am very stuck on what to do next. Ive been in contact with the papers authors but they havent been a huge help. If anyone has a suggestion, it would be a huge help.
Here is the code:
clear; close all;
% Initial conditions and parameters
x0 = 0; %m
xfinal = 1; %m
h0 = 0.11; %m
hx0 = 0;
hxx0 = 0;
s0 = 1/163; % bottom slope
F1 = 1.07; % Froude number
N = 0.142857;
K = 2; % curvature
q = 0.12; % m^2 / s
g = 9.81; % gravity
k = 0;
theta = 0.00009;
beta0 = (1 + N)^2 / (1 + 2 * N);
ep0 = (h0 * hxx0) / (1 + hx0^2);
ep1 = hx0^2 / (1 + hx0^2);
S0 = (h0^2 / 2) + beta0 * (q^2 / (g * h0)) * (1 * ep0 * ((2 / (K + N + 2)) – ((1 – 2 * N) / (K + 2 * N + 2))) – ep1 * ((N + 1) / (N + 3)));
f0 = [h0; hx0; S0; theta; k;0];
% Solve the ODE
options = odeset(‘RelTol’,1e-5,’AbsTol’,1e-7);
[x, f] = ode23tb(@undularODE, [x0 xfinal], f0, options);
% Plot the results
figure;
plot(x, f);
legend(‘h’, ‘hx’, ‘S’, ‘theta’, ‘k’,’N’);
xlabel(‘x’);
ylabel(‘y’);
title(‘Solution’);
Here is the UdularODE function:
function dfdx = undularODE(x, f)
% Parameters
g = 9.81;
N = 0.142857;
q = 0.12; % m^2 / s
K = 2;
s0 = 1/163; % bottom slope
utheta = 1.1; % m/s velocity at momentum thickness
theta = f(4);
Nu = 1.12*10^-6; % kinematic viscosity
% Computed parameters
beta = (1 + N)^2 / (1 + 2 * N);
ep1 = f(2)^2 / (1 + f(2)^2); % epsilon 1
Re = 178458; % Reynolds number
Cf = 0.075 / (log(Re) – 2)^2; % skin friction coeff
ff = 4 * Cf * (1 + N)^2; % friction factor
sf = ff / 8 * (1 + 2 * (f(1) / 1)) * 1.47^2; % friction slope
Ue = (1 + N) * (q / f(1)); % max velocity at boundary layer edge
k = utheta / Ue;
H = (1.3 + 1.3 * (0.7 – k) + 3 * (0.7 – k)^2)^(2/3); % shape factor
Rtheta = (Ue * theta) / Nu; % momentum thickness Reynolds number
N2 = (H-1)/2; % calculated value of N
dUedx = -1*Ue*f(6)*(1+N2)*(Ue/f(1))*f(2);%dUedx
gamma = (theta / Ue) * dUedx * Rtheta^(1/4); % Buri Shape Factor
dHdx = (3*k^2 – 5.5*k + 3.68)^(2/3) * ( (1.46*(1-k^2))/( (1.5+k)*theta*Rtheta^(1/4)* (gamma +0.118*(0.67-k)) ) );
%dNdx(end+1)= 0.5*dHdx;
% Define the system of ODEs
hx = f(2);
eq2 = (((g * f(1)) / (beta * q^2)) * (f(3) – (f(1)^2 / 2)) – 1 + ep1 * ((N + 1) / (N + 3))) / …
((f(1) / (1 + f(2)^2)) * (2 / (K + N + 2)) – ((1 – 2 * N) / (K + 2 * N + 2)));
eq3 = f(1) * (s0 – sf);
eq4 = -1 * f(4)*((1 / (2 + 2 * N)) * (dHdx) – (1 / f(1)) * hx) * (H + 2) + Cf / 2;
eq5 = 1.46 * ((1 – f(5)^2) / ((1.5 + f(5)) * theta * Rtheta^(1/4))) * (gamma + 0.118 * (0.67 – f(5)));
eq6 = 0.5*dHdx;
% Return the derivatives
dfdx = [hx;
eq2;
eq3;
eq4;
eq5;
eq6];
%keyboard
endHello all, I am working on solving a 6 equation ODE to model the shape of a free surface in a channel. The paper that I am referencing only uses 5 euqations but i have no idea how they are solving for one of the variables so I added the variable and equation to the ODE. I have tried almost solver avalible in matlab and am very stuck on what to do next. Ive been in contact with the papers authors but they havent been a huge help. If anyone has a suggestion, it would be a huge help.
Here is the code:
clear; close all;
% Initial conditions and parameters
x0 = 0; %m
xfinal = 1; %m
h0 = 0.11; %m
hx0 = 0;
hxx0 = 0;
s0 = 1/163; % bottom slope
F1 = 1.07; % Froude number
N = 0.142857;
K = 2; % curvature
q = 0.12; % m^2 / s
g = 9.81; % gravity
k = 0;
theta = 0.00009;
beta0 = (1 + N)^2 / (1 + 2 * N);
ep0 = (h0 * hxx0) / (1 + hx0^2);
ep1 = hx0^2 / (1 + hx0^2);
S0 = (h0^2 / 2) + beta0 * (q^2 / (g * h0)) * (1 * ep0 * ((2 / (K + N + 2)) – ((1 – 2 * N) / (K + 2 * N + 2))) – ep1 * ((N + 1) / (N + 3)));
f0 = [h0; hx0; S0; theta; k;0];
% Solve the ODE
options = odeset(‘RelTol’,1e-5,’AbsTol’,1e-7);
[x, f] = ode23tb(@undularODE, [x0 xfinal], f0, options);
% Plot the results
figure;
plot(x, f);
legend(‘h’, ‘hx’, ‘S’, ‘theta’, ‘k’,’N’);
xlabel(‘x’);
ylabel(‘y’);
title(‘Solution’);
Here is the UdularODE function:
function dfdx = undularODE(x, f)
% Parameters
g = 9.81;
N = 0.142857;
q = 0.12; % m^2 / s
K = 2;
s0 = 1/163; % bottom slope
utheta = 1.1; % m/s velocity at momentum thickness
theta = f(4);
Nu = 1.12*10^-6; % kinematic viscosity
% Computed parameters
beta = (1 + N)^2 / (1 + 2 * N);
ep1 = f(2)^2 / (1 + f(2)^2); % epsilon 1
Re = 178458; % Reynolds number
Cf = 0.075 / (log(Re) – 2)^2; % skin friction coeff
ff = 4 * Cf * (1 + N)^2; % friction factor
sf = ff / 8 * (1 + 2 * (f(1) / 1)) * 1.47^2; % friction slope
Ue = (1 + N) * (q / f(1)); % max velocity at boundary layer edge
k = utheta / Ue;
H = (1.3 + 1.3 * (0.7 – k) + 3 * (0.7 – k)^2)^(2/3); % shape factor
Rtheta = (Ue * theta) / Nu; % momentum thickness Reynolds number
N2 = (H-1)/2; % calculated value of N
dUedx = -1*Ue*f(6)*(1+N2)*(Ue/f(1))*f(2);%dUedx
gamma = (theta / Ue) * dUedx * Rtheta^(1/4); % Buri Shape Factor
dHdx = (3*k^2 – 5.5*k + 3.68)^(2/3) * ( (1.46*(1-k^2))/( (1.5+k)*theta*Rtheta^(1/4)* (gamma +0.118*(0.67-k)) ) );
%dNdx(end+1)= 0.5*dHdx;
% Define the system of ODEs
hx = f(2);
eq2 = (((g * f(1)) / (beta * q^2)) * (f(3) – (f(1)^2 / 2)) – 1 + ep1 * ((N + 1) / (N + 3))) / …
((f(1) / (1 + f(2)^2)) * (2 / (K + N + 2)) – ((1 – 2 * N) / (K + 2 * N + 2)));
eq3 = f(1) * (s0 – sf);
eq4 = -1 * f(4)*((1 / (2 + 2 * N)) * (dHdx) – (1 / f(1)) * hx) * (H + 2) + Cf / 2;
eq5 = 1.46 * ((1 – f(5)^2) / ((1.5 + f(5)) * theta * Rtheta^(1/4))) * (gamma + 0.118 * (0.67 – f(5)));
eq6 = 0.5*dHdx;
% Return the derivatives
dfdx = [hx;
eq2;
eq3;
eq4;
eq5;
eq6];
%keyboard
end Hello all, I am working on solving a 6 equation ODE to model the shape of a free surface in a channel. The paper that I am referencing only uses 5 euqations but i have no idea how they are solving for one of the variables so I added the variable and equation to the ODE. I have tried almost solver avalible in matlab and am very stuck on what to do next. Ive been in contact with the papers authors but they havent been a huge help. If anyone has a suggestion, it would be a huge help.
Here is the code:
clear; close all;
% Initial conditions and parameters
x0 = 0; %m
xfinal = 1; %m
h0 = 0.11; %m
hx0 = 0;
hxx0 = 0;
s0 = 1/163; % bottom slope
F1 = 1.07; % Froude number
N = 0.142857;
K = 2; % curvature
q = 0.12; % m^2 / s
g = 9.81; % gravity
k = 0;
theta = 0.00009;
beta0 = (1 + N)^2 / (1 + 2 * N);
ep0 = (h0 * hxx0) / (1 + hx0^2);
ep1 = hx0^2 / (1 + hx0^2);
S0 = (h0^2 / 2) + beta0 * (q^2 / (g * h0)) * (1 * ep0 * ((2 / (K + N + 2)) – ((1 – 2 * N) / (K + 2 * N + 2))) – ep1 * ((N + 1) / (N + 3)));
f0 = [h0; hx0; S0; theta; k;0];
% Solve the ODE
options = odeset(‘RelTol’,1e-5,’AbsTol’,1e-7);
[x, f] = ode23tb(@undularODE, [x0 xfinal], f0, options);
% Plot the results
figure;
plot(x, f);
legend(‘h’, ‘hx’, ‘S’, ‘theta’, ‘k’,’N’);
xlabel(‘x’);
ylabel(‘y’);
title(‘Solution’);
Here is the UdularODE function:
function dfdx = undularODE(x, f)
% Parameters
g = 9.81;
N = 0.142857;
q = 0.12; % m^2 / s
K = 2;
s0 = 1/163; % bottom slope
utheta = 1.1; % m/s velocity at momentum thickness
theta = f(4);
Nu = 1.12*10^-6; % kinematic viscosity
% Computed parameters
beta = (1 + N)^2 / (1 + 2 * N);
ep1 = f(2)^2 / (1 + f(2)^2); % epsilon 1
Re = 178458; % Reynolds number
Cf = 0.075 / (log(Re) – 2)^2; % skin friction coeff
ff = 4 * Cf * (1 + N)^2; % friction factor
sf = ff / 8 * (1 + 2 * (f(1) / 1)) * 1.47^2; % friction slope
Ue = (1 + N) * (q / f(1)); % max velocity at boundary layer edge
k = utheta / Ue;
H = (1.3 + 1.3 * (0.7 – k) + 3 * (0.7 – k)^2)^(2/3); % shape factor
Rtheta = (Ue * theta) / Nu; % momentum thickness Reynolds number
N2 = (H-1)/2; % calculated value of N
dUedx = -1*Ue*f(6)*(1+N2)*(Ue/f(1))*f(2);%dUedx
gamma = (theta / Ue) * dUedx * Rtheta^(1/4); % Buri Shape Factor
dHdx = (3*k^2 – 5.5*k + 3.68)^(2/3) * ( (1.46*(1-k^2))/( (1.5+k)*theta*Rtheta^(1/4)* (gamma +0.118*(0.67-k)) ) );
%dNdx(end+1)= 0.5*dHdx;
% Define the system of ODEs
hx = f(2);
eq2 = (((g * f(1)) / (beta * q^2)) * (f(3) – (f(1)^2 / 2)) – 1 + ep1 * ((N + 1) / (N + 3))) / …
((f(1) / (1 + f(2)^2)) * (2 / (K + N + 2)) – ((1 – 2 * N) / (K + 2 * N + 2)));
eq3 = f(1) * (s0 – sf);
eq4 = -1 * f(4)*((1 / (2 + 2 * N)) * (dHdx) – (1 / f(1)) * hx) * (H + 2) + Cf / 2;
eq5 = 1.46 * ((1 – f(5)^2) / ((1.5 + f(5)) * theta * Rtheta^(1/4))) * (gamma + 0.118 * (0.67 – f(5)));
eq6 = 0.5*dHdx;
% Return the derivatives
dfdx = [hx;
eq2;
eq3;
eq4;
eq5;
eq6];
%keyboard
end fluid, ode, ode45, unstable MATLAB Answers — New Questions
