I have these warnings in my code. Can anybody help me? Thanks
1 view (last 30 days)
Show older comments
I attached the code and a .txt. Before using the code you have to save the .txt as .DAT since I put it like that on Matlab
Attencion! You need to change comas into points in the .DAT file
6 Comments
Walter Roberson
on 20 Oct 2022
Moving the reading out of the loop does not appear to present any problems for me.
rpms = [1600, 1800, 2000, 2200, 2400, 2600, 2800, 3000];
dz = 0.2; % Diferencia de cotas entre la entrada y la salida = 20 cm
rho = 1e3; % Densidad = 1000 kg/m3
g = 9.8; % Aceleracio´n de la gravedad
D = 0.2; % Longitud caracter´ıstica del impulsor
% Creamos vectores vac´ıos para almacenar los datos:
vecq = []; vech = []; vecwu = []; vecwt = []; veceta = [];
% CAUDAL ALTURA POT UTIL POT EJE RENDIMIENTO
run = readmatrix('https://www.mathworks.com/matlabcentral/answers/uploaded_files/1162238/robainutile.txt',...
'DecimalSeparator',','); %##### MODIFIED TO RUN HERE
% Convertimos de RPM a RAD/S
run(:,1) = run(:,1)*2*pi/60;
% Convertimos de BAR a PASCAL
run(:,2:3) = run(:,2:3)*1e5;
% La columna 4, representa el par en SI => no se toca
% Convertimos de L/MIN a M3/S
run(:,5) = run(:,5)*1e-3/60;
% Definimos un vector 25 caudal con 500 puntos:
q = linspace(0,max(run(:,5)),500);
rundata = run;
for i=length(rpms):-1:1
% Definimos el nombre de los ficheros que hay que leer
run = rundata;
% Ordenamos los datos en orden Q creciente
[aux, ind] = sort(run(:,5));
run = run(ind,:);
% Calculamos la altura manome´trica:
gHm = (run(:,3)-run(:,2))/rho + g*dz;
% Ajuste de gHm con una funcion polinomica de segundo orden
[p,S,mu] = polyfit(run(:,5),gHm,2);
% Interpolamos gHm para valores de caudal q:
[y,delta] = polyval(p,q',S,mu);
% Almacenamos los resultados en un vector vecq y vech
vecq = [vecq ; q']; vech = [vech; y];
% Calculamos la potencia útil:
Wu = rho*run(:,5).*gHm;
% Ajustamos y la interpolamos
[pwu,Swu,muwu] = polyfit(run(:,5),Wu,3);
Wufit=polyval(pwu,q',Swu,muwu);
%Almacenamos en vector
vecwu = [vecwu ; Wufit];
% Calculamos la potencia en el eje
Wt = run(:,1).*run(:,4);
% y la interpolamos
[pwt,Swt,muwt] = polyfit(run(:,5),Wt,1);
Wtfit=polyval(pwt,q',Swt,muwt);
% Almacenamos los resultados en un vector vecwt
vecwt = [vecwt ; Wtfit];
% Calculamos rendimiento
eta = Wu./Wt*100;%rho*q.*y' ./Wtfit'*100;
% Interpolamos rendimiento
[peta,Seta,mueta] = polyfit(run(:,5),eta,2);
etafit=polyval(peta,q',Seta,mueta);
% Almacenamos los resultados en un vector veceta
veceta = [veceta ; etafit];
% Hacemos las curvas dimensionales
figure(1); hold on;
xlabel('Q (m^3/s)');ylabel('gHm (m^2/s^2)');
plot(run(:,5),gHm,'ob'); hold on;
plot(q,y,'r-');
grid on;box on;
figure(2); hold on;
xlabel('Q (m^3/s)');ylabel('Wu (W)');
plot(run(:,5),Wu,'ob'); hold on;
plot(q,Wufit,'r-');
grid on;box on;
figure(3); hold on;
xlabel('Q (m^3/s)');ylabel('Wt (W)');
plot(run(:,5),Wt,'ob'); hold on;
plot(q,Wtfit,'r-');
grid on; box on;
figure(4); hold on;
xlabel('Q (m^3/s)');ylabel('\eta (%)');
plot(run(:,5),eta,'ob'); hold on;
plot(q,etafit,'r-');
xlabel('Q (m^3/s)');ylabel('\eta (%)');
grid on; box on;
end
m = 101;
n = 100;
ejeQ = linspace(0,max(vecq),m);
ejeWu = linspace(0,max(vecwu),n);
ejeWt = linspace(0,max(vecwt),n);
ejeHm = linspace(min(vech),max(vech),n);
[QQ,HH] = meshgrid(ejeQ,ejeHm);
[QQ,WW] = meshgrid(ejeQ,ejeWt);
ETAqh = griddata(vecq,vech,veceta,QQ,HH,'linear');
ETAqw = griddata(vecq,vecwt,veceta,QQ,WW,'linear');
indETAqh=isnan(ETAqh);
indETAqw=isnan(ETAqw);
figure;
plot3(vecq,vech/g,veceta,'k.')
hold on
ETAqh2=gridfit(vecq,vech,veceta,ejeQ,ejeHm);
ETAqh2(find(indETAqh==1)) = NaN;
surf(QQ,HH/g,ETAqh2);colormap(jet);colorbar;
camlight right; lighting phong; shading interp;
xlabel('Q (m^3/s)'); ylabel('H_m (m)'); zlabel('\eta (%)');
title('Rendimientos usando GRIDFIT')
box on
figure;
plot3(vecq,vecwt,veceta,'k.')
hold on
ETAqw2=gridfit(vecq,vecwt,veceta,ejeQ,ejeWt);
ETAqw2(find(indETAqw==1)) = NaN;
surf(QQ,WW,ETAqw2);colormap(jet);colorbar;
camlight right; lighting phong; shading interp;
xlabel('Q (m^3/s)'); ylabel('W_t (W)'); zlabel('\eta (%)');
title('Rendimientos usando GRIDFIT')
box on
figure;
[C,h] = contourf(QQ,HH/g,ETAqh2,[0 5 10 15 20 25 30 32.5 35]);colormap(jet);colorbar;
text_handle = clabel(C,h,'FontSize',14,'Color','k','Rotation',0,'LabelSpacing',500);
set(text_handle,'BackgroundColor',[1 1 1])
xlabel('Q (m^3/s)'); ylabel('H_m (m)');
title('Isorendimientos usando GRIDFIT')
figure;
[C,h] = contourf(QQ,WW,ETAqw2,[0 5 10 15 20 25 30 32.5 35]);colormap(jet);colorbar;
text_handle = clabel(C,h,'FontSize',14,'Color','k','Rotation',0,'LabelSpacing',500);
set(text_handle,'BackgroundColor',[1 1 1])
xlabel('Q (m^3/s)'); ylabel('W_t (W)');
title('Isorendimientos usando GRIDFIT')
Answers (1)
Walter Roberson
on 20 Oct 2022
for i=length(rpms):-1:1
The calculations in your loop do not depend upon i, so you are generating exactly the same data every time through the loop. That is why you are getting duplicates and not getting the output you want.
5 Comments
Walter Roberson
on 24 Oct 2022
I am going to assume that gridfit() is Surface Fitting using gridfit from the File Exchange, by @John D'Errico
Walter Roberson
on 24 Oct 2022
I do not know why you are getting the warning about rank deficient inside of gridfit()
@John D'Errico -- the problem is occuring on
else
% no mask
zgrid = reshape(A\rhs,ny,nx);
after 'gradients'. A is sparse 20700 by 10100 but rank(A) is 10099 -- one column has a problem.
See Also
Categories
Find more on Computational Geometry in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!