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