- 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.

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

3 views (last 30 days)

Show older comments

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

##### 0 Comments

### Answers (1)

Sahas
on 5 Sep 2024

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”:

Hope this is beneficial!

### See Also

### Community Treasure Hunt

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

Start Hunting!