Why did I get the invalid Q values (inf, -inf, NaN) from following code?

This is my equation,

I got following by applying implicit finite differencce method,

Here and

And this id my code

t = 0:300:(8*60+15)*60; % time 
x = [0,3750,7500,8000]; % distance 
[xq,tq] = meshgrid(0:100:8000,0:60:(8*60+15)*60);

eta = load('Eta.txt');
eta_intp = interp2(x,t,eta,xq,tq,'spline');

% Fit a polynomial to the interpolated data

% Flatten the xq and tq grids to create predictor matrix
x_t_flat = [xq(:),tq(:)]; 

% Flatten the eta_interpolated matrix to create response vector
eta_flat = eta_intp(:); 

% Fit a polynomial function
degree = 9;
model = polyfitn(x_t_flat,eta_flat,degree);

% Define a function to evaluate the polynomial at each (xq, tq) point
eta_poly = @(xq_val, tq_val) polyvaln(model,[xq_val,tq_val]);

% Use arrayfun to evaluate the polynomial at all grid points
eta_poly_values = arrayfun(eta_poly,xq,tq);

figure (3)
surf(xq,tq,eta_poly_values,'Edgecolor','none');
xlabel('Distance (m)');
ylabel('Time (s)');
zlabel('Water Surface Elevation (m)');
title('Interpolated Eta Values');

dx = 100;
dt = 60;

L = 8000;
T1 = (8*60+15)*60;
M = round(L/dx +1); % space nodes 
N = round(T1/dt +1); % time nodes

g = 9.81; % gravitational accelaration 
B = 2000; % cross sectional width
H = 4; % undisturbed depth
Cd = 0.0022; % drag coefficient
beta = dt/dx;

A = zeros(N,M);

for n = 1:N
    for i = 1:M
        A(n,i) = (B*(H+eta_intp(n,i)));
    end
end

Q = zeros(N,M);
Q(:,1) = 0; % left boundary condition Q(0,t)=0
% no right boundary condition
Q(1,2:M) = 0; % initial condition without left boundary condition Q(x,0)=0

for n = 1:N-1
    for i = 2:M
        Q(n+1,i) = (-g*A(n,i)*(eta_poly_values(n,i)-eta_poly_values(n,i-1)) -(Q(n,i)*(Q(n,i)-Q(n,i-1))/(A(n,i))))*beta -((B*dt*Cd/A(n,i))*abs(Q(n,i)/A(n,i)) -1)*Q(n,i);
    end
end

figure (5)
surf(xq,tq,Q,'Edgecolor','none')
title('Q values using 1-D Momentum Equation')
xlabel('Distance (m)');
ylabel('Time (s)');
zlabel('Q (m^3/s)');

Answer 1

[Sahas]

Hi @SPerera,

As per my understanding, the equation mentioned gives the “flow discharge”, represented by “Q”, for the next time step and its effects from various components such as gravity, inertia, and friction.

The following reasons might be the cause for getting invalid values for “Q”:

• Ensure that “A(n, i)”, which is in the denominator of the inertial term of the flow discharge equation, is not zero. Otherwise, it will result in “NaN”.
• Ensure that the initial condition for the flow discharge equation, Q(n, i-1), is defined for all “i” to avoid boundary condition problems which can also lead to a “NaN”.
• Check the inputs for “Q” value calculations are not invalid. You can do this by using “isnan” or “”isinf” functions on “eta_poly_values” and “A” matrices.
• Sometimes a high degree polynomial (degree 9 polynomial in the code provided) can lead to overfitting during interpolation operations, resulting in “inf” value. Try fitting with a lower degree polynomial by adjusting the “degree” variable in the code.

Hope this is beneficial!

[SPerera]

Hi @Sahas,

Yes, this is 1-D momentum equation. I checked the validity of “eta_poly_values” and “A” matrices. I got zero matrices as the results. Also I adjusted the degree to 2, but got the inavalid values. Actually I got higher values for Q. Values increase rapidly and become "inf". But it should be less than 50 . Are there any other issues in my code?