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);

Warning: Failure at t=3.074452e-03. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (1.092265e-17) at time t.

% 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

[Torsten]

Why do you use

k = utheta / Ue;

and not

k = f(5)

? That's inconsistent.

Further if

H = (1.3 + 1.3 * (0.7 - k) + 3 * (0.7 - k)^2)^(2/3)

, then

dHdx = -2/3*(1.3 + 1.3 * (0.7 - k) + 3 * (0.7 - k)^2)^(-1/3) * (1.3 + 3*2*(0.7 - k)) * dk/dx

Does this equal your dHdx ?

[Sam Chak]

On top of that, there is another inconsistency:

dUedx = -1*Ue*f(6)*(1 + N2)*(Ue/f(1))*f(2);

but in the paper

[Sam Chak]

Also, since H is a function of k in Eq. (15), then

.

The 6th state equation is probably unnecessary because

can be directly derived

and, will also be a function of k.

Because both depend on the state k, only is necessary, which is the 5th state equation in Eq. (13):

I understand why you coded .

[Sam Chak]

Also, below Eq. (10), you can see how U is defined. This U is needed in Eq. (16), but you coded it differently.

Note that depends on the solution of , and depends on both solutions of and , and in Eq. (11) requires the solution of .

(Edited to correct mistakes pointed out by @Torsten).

Since the authors didn't make as the 6th ODE, most likely they used the approximation of .

Try fixing the code below:

function dfdx = undularODE(x, f)
    % definitions
    h       = f(1);
    hx      = f(2);
    S       = f(3);
    theta   = f(4);
    k       = f(5);

    % Independent parameters
    g       = 9.81;
    % N     = 0.142857;     % Probably not needed?
    q       = 0.12;         % m^2 / s
    K       = 2;
    s0      = 1/163;        % bottom slope
    utheta  = 1.1;          % m/s velocity at momentum thickness 
    Nu      = 1.12*10^-6;   % kinematic viscosity
    Re      = 178458;                                           % Reynolds number
    b       = 1;
    F       = 1.47;
    
    % Dependent parameters
    ep1     = hx^2/(1 + hx^2);                              % epsilon 1
    Cf      = 0.075/(log(Re) - 2)^2;                        % skin friction coeff
    ff      = 4*Cf*(1 + N)^2;                               % friction factor
    sf      = (ff/8)*(1 + 2*h/b)*F^2;                       % friction slope 
    % k       = utheta / Ue;
    H       = (1.3 + 1.3*(0.7 - k) + 3*(0.7 - k)^2)^(2/3);  % shape factor
    N       = (H - 1)/2;                                    % calculated value of N 
    beta    = (1 + N)^2/(1 + 2*N);
    Ue      = (1 + N)*(q/h);                                % max velocity at boundary layer edge
    Rtheta  = (Ue * theta) / Nu;                            % momentum thickness Reynolds number
    
    dHdx    = ...;      % depends on the solution  of k(x)
    dUedx   = ...;      % depends on the solutions of k(x) and h(x)
    gamma   = (theta / Ue) * dUedx * Rtheta^(1/4);          % Buri Shape Factor
    % dNdx(end+1)= 0.5*dHdx;
    
    num2    = (g*h)/(beta*q^2)*(S - h^2/2) - 1 + ep1*(N + 1)/(N + 3);
    den2    = h/(1 + hx^2)*(2/(K + N + 2) - (1 - 2*N)/(K + 2*N + 2));

    % Define the system of ODEs
    dhdx    = hx;
    dhxdx   = num2/den2;
    dSdx    = h*(s0 - sf);
    dtheta  = - theta*(dHdx/(2 + 2*N) - hx/h)*(H + 2) + Cf/2;
    dkdx    = 1.46*(1 - k^2)/((1.5 + k)*theta*Rtheta^(1/4))*(gamma + 0.118*(0.67 - k));
    
    % eq6     = 0.5*dHdx;

    % Return the derivatives
    dfdx    = [dhdx
               dhxdx
               dSdx
               dtheta
               dkdx];
end

[Torsten]

N = (H-1)/2

instead of

N       = 0.142857;

And I don't think that derivatives depend on initial values:

dHdx    = ...;      % depends on the initial value  of k 
dUedx   = ...;      % depends on the initial values of k and h

[Sam Chak]

@Torsten is right. I incorrectly phrased the comments. When typing those comments, my mind was thinking whether there is an algebraic loop dependency between k and because .

[Sam Chak]

Please review the computation of Cf in the code against the formula in Eq. (14):

[Joseph Improta]

@Torsten @Sam Chak thank you for your suggestions. Im gonna go though them and work on myt code. I will get back to you with what I find. I will say this is how my mind interpreted it initally: to calculate dH/dx. I am doing the chain rule with dH/dk and dk/dx to obtain dH/dx. To calculate dk/dx though, I need gamma which needs dUe/dx which then needs dN/dx. The only way I have found to calculate dN/dx is to do the chain rule with dN/dH and dH/dx.

[Torsten]

Insert the expression for gamma in equation (13). Then use

dN/dx = 0.5*dH/dx = 0.5*(-2/3*(1.3 + 1.3 * (0.7 - k) + 3 * (0.7 - k)^2)^(-1/3) * (1.3 + 3*2*(0.7 - k)) * dk/dx)

and insert it in the expression for dUe/dx.

You get an equation of the form

dk/dx = U*0.5*(-2/3*(1.3 + 1.3 * (0.7 - k) + 3 * (0.7 - k)^2)^(-1/3) * (1.3 + 3*2*(0.7 - k)) * dk/dx) - (1+N)*U/h * dh/dx

Insert f(2) for dh/dx and solve for dk/dx.

[Sam Chak]

Hi @Joseph Improta

Following @Torsten's approach, by tracing the circular dependency from and find the explicit solution for , you should be able to break the algebraic loop. Logically, you should use the newly derived expression of for the 5th ODE, replacing Eq. (13).

where

[Sam Chak]

Hi @Joseph Improta

To reduce human error in hand calculation, you should attempt finding the explicit solution for using the 'syms', 'diff()' and the 'solve()' commands. Check if I entered the equations correctly.

In the code, U is , R is , h1 is h and h2 is .

syms k(x) U(x) N(x) H theta R q h1 h2 dk

N    = (H - 1)/2;
H    = (1.3 + 1.3*(0.7 - k) + 3*(0.7 - k)^2)^(2/3);
dNdH = 1/2;                                 % dN/dH
dHdk = diff(H, k)                           % dH/dk

dHdk(x) = 

dN   = dNdH*dHdk*dk;                        % dN/dx
dU   = (q/h1)*dN - (1 + N)*q/(h1^2)*h2;     % dUe/dx
Gamma= (theta/U)*dU*R;

eqn  = dk == 1.46*(1 - k^2)/((1.5 + k)*theta*R)*(Gamma + 0.118*(0.67 - k));
dk   = solve(eqn, dk)                       % dk/dx

dk = 

[Joseph Improta]

@Torsten @Sam ChakI edited the function with what you guys said. I am confused on some parts of it but here is what it looks like now:

function dfdx = undularODE(x, f)
    % definitions
    h       = f(1);
    hx      = f(2);
    S       = f(3);
    theta   = f(4);
    k       = f(5);

    % Independent parameters
    g       = 9.81;
    % N     = 0.142857;     % Probably not needed?
    q       = 0.12;         % m^2 / s
    K       = 2;
    s0      = 1/163;        % bottom slope
    utheta  = 1.1;          % m/s velocity at momentum thickness 
    Nu      = 1.12*10^-6;   % kinematic viscosity
    Re      = 178458;                                           % Reynolds number
    b       = 1;
    F       = 1.47;
    
    % Dependent parameters
    ep1     = hx^2/(1 + hx^2);                              % epsilon 1
    Cf      = 0.075/(log(Re) - 2)^2;                        % skin friction coeff
                         % friction slope 
    
    H       = (1.3 + 1.3*(0.7 - k) + 3*(0.7 - k)^2)^(2/3);  % shape factor
    N       = (H - 1)/2;  % calculated value of N 
    ff      = 4*Cf*(1 + N)^2;                               % friction factor
    sf      = (ff/8)*(1 + 2*h/b)*F^2;  

    beta    = (1 + N)^2/(1 + 2*N);
    Ue      = (1 + N)*(q/h); 
    k       = utheta / Ue; % max velocity at boundary layer edge
    Rtheta  = (Ue * theta) / Nu;                            % momentum thickness Reynolds number
    
    dHdx    = -2/3*(1.3 + 1.3 * (0.7 - k) + 3 * (0.7 - k)^2)^(-1/3) * (1.3 + 3*2*(0.7 - k)) * f(5);      % depends on the solution  of k(x)
    dNdx = 0.5*dHdx;
    dUedx   = Ue * dNdx;      % depends on the solutions of k(x) and h(x)
    gamma   = (theta / Ue) * dUedx * Rtheta^(1/4);          % Buri Shape Factor
    % dNdx(end+1)= 0.5*dHdx;
    
    num2    = (g*h)/(beta*q^2)*(S - h^2/2) - 1 + ep1*(N + 1)/(N + 3);
    den2    = h/(1 + hx^2)*(2/(K + N + 2) - (1 - 2*N)/(K + 2*N + 2));

    % Define the system of ODEs
    dhdx    = hx;
    dhxdx   = num2/den2;
    dSdx    = h*(s0 - sf);
    dtheta  = - theta*(dHdx/(2 + 2*N) - hx/h)*(H + 2) + Cf/2;
    dkdx    =(((73*k^2)/50 - 73/50)*((59*k)/500 + (Rtheta*f(2)*q*theta*(H/2 + 1/2))/((f(1))^2*Ue) - 3953/50000))/(Rtheta*theta*(k + 3/2)*((q*((73*k^2)/50 - 73/50)*(6*k - 11/2))/(3*f(1)*Ue*(k+ 3/2)*(3*(k - 7/10)^2 - (13*k)/10 + 221/100)^(1/3)) + 1));
    
    % eq6     = 0.5*dHdx;

    % Return the derivatives
    dfdx    = [dhdx
               dhxdx
               dSdx
               dtheta
               dkdx];
end

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];

% Solve the ODE
options = odeset('RelTol',1e-5,'AbsTol',1e-7);
[x, f] = ode23tb(@undularODE, [x0 xfinal], f0, options);

Warning: Failure at t=7.481200e-02. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (2.657856e-16) at time t.

% Plot the results
figure;
plot(x, f);
legend('h', 'hx', 'S', 'theta', 'k','N');

Warning: Ignoring extra legend entries.

xlabel('x');
ylabel('y');
title('Solution');

Also I don't know which solver I should be using. I started with ode45