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

5 次查看（过去 30 天）
SPerera
SPerera 2024-9-5
评论： SPerera 2024-9-6
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)');

请先登录，再回答此问题。

回答（1 个）

Sahas
Sahas 2024-9-5
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!
  1 个评论
SPerera
SPerera 2024-9-6
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?

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Programming 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by