Filter löschen
Filter löschen

MATLAB error using fzero function

8 Ansichten (letzte 30 Tage)
Alina Abdikadyr
Alina Abdikadyr am 23 Jan. 2023
Kommentiert: Star Strider am 27 Jan. 2023
My following code generates the plot of V and D values in figure 1. In the graph, the parabolas and straight lines intersect, and I need to find the roots from the plot. So I tried to use fzero function, but the error appeared: Operands to the logical AND (&&) and OR (||) operators must be convertible to logical scalar values. Use the ANY or ALL functions to reduce operands to logical scalar values.
Error in fzero (line 326) elseif ~isfinite(fx) || ~isreal(fx)
Error in HW1 (line 35) x=fzero(fun,1);
My code is:
clear all; close all
W = 10000; %[N]
S = 40; %[m^2]
AR = 7;
cd0 = 0.01;
k = 1 / pi / AR;
clalpha = 2*pi;
Tsl=800;
figure(1);hold on; xlabel('V');ylabel('D')
for h=0:1:8;
i=0;
for alpha = 1:0.25:12
i=i+1;
rho(i)=1.2*exp(-h/10.4);
cl(i) = clalpha * alpha * pi/180;
V(i) = sqrt(2*W/rho(i)/S/cl(i));
L(i) = 0.5 * rho(i) * V(i) * V(i) * S * cl(i);
cd(i) = cd0 + k * cl(i) * cl(i);
D(i) = 0.5 * rho(i) * V(i) * V(i) * S * cd(i);
clcd(i) = cl(i)/cd(i);
p(i) = D(i)*V(i);
ang(i) = alpha;
T(i)=Tsl*(rho(i)/1.2).^0.75;
end
figure(1); plot(V,D); hold on
plot(V,T);
end
fun = @(V) 0.5*V.*V.*rho.*S.*cd-T;
x=fzero(fun,1);
Operands to the logical AND (&&) and OR (||) operators must be convertible to logical scalar values. Use the ANY or ALL functions to reduce operands to logical scalar values.

Error in fzero (line 326)
elseif ~isfinite(fx) || ~isreal(fx)
Probably, I should not use the fzero function, but the task is to find the roots of V from a plot (figure 1). There are parabolas and straight lines respectively.

Akzeptierte Antwort

Star Strider
Star Strider am 23 Jan. 2023
I have no idea what intersections you want to find, so this is a guess.
It should at least get you started —
% clear all; close all
W = 10000; %[N]
S = 40; %[m^2]
AR = 7;
cd0 = 0.01;
k = 1 / pi / AR;
clalpha = 2*pi;
Tsl=800;
figure(1);hold on; xlabel('V');ylabel('D')
hv=0:1:8;
alphav = 1:0.25:12;
for k1 = 1:numel(hv)
h = hv(k1);
i=0;
for alpha = 1:0.25:12
i=i+1;
rho(i)=1.2*exp(-h/10.4);
cl(i) = clalpha * alpha * pi/180;
V(i) = sqrt(2*W/rho(i)/S/cl(i));
L(i) = 0.5 * rho(i) * V(i) * V(i) * S * cl(i);
cd(i) = cd0 + k * cl(i) * cl(i);
D(i) = 0.5 * rho(i) * V(i) * V(i) * S * cd(i);
clcd(i) = cl(i)/cd(i);
p(i) = D(i)*V(i);
ang(i) = alpha;
T(i)=Tsl*(rho(i)/1.2).^0.75;
end
Vm(k1,:) = V; % Save Results
Dm(k1,:) = D; % Save Results
Tm(k1,:) = T; % Save Results
xix = find(diff(sign(0.5*V.*V.*rho.*S.*cd-T))); % Approximate Indices Of Intersections
if ~isempty(xix) % Skip If Empty
for k3 = 1:numel(xix)
idxrng = max(xix(k3)-1,1) : min(numel(V),xix(k3)+1);
xv(k3) = interp1(D(idxrng)-T(idxrng), V(idxrng), 0);
yv(k3) = interp1(V(idxrng),D(idxrng), xv(k3));
Xm{k1,k3} = xv(k3); % Intersection X-Matrix
Ym{k1,k3} = yv(k3); % Intersection Y-Matrix
end
end
figure(1)
plot(V,D, 'DisplayName',sprintf('D_{%d}',k1))
hold on
plot(V,T, 'DisplayName',sprintf('T_{%d}',k1))
plot(xv, yv, 's', 'DisplayName',sprintf('Intersections %d',k1)) % Plot Intersections (Optional)
end
legend('Location','eastoutside','NumColumns',2)
Intersections = cell2table(cat(2,Xm,Ym), 'VariableNames',{'X1','X2','Y1','Y2'})
Intersections = 9×4 table
X1 X2 Y1 Y2 ______ ____________ ______ ____________ 55.445 {0×0 double} 800 {0×0 double} 55.654 {0×0 double} 744.34 {0×0 double} 55.885 {0×0 double} 692.55 {0×0 double} 55.914 {[ 21.1789]} 644.37 {[644.3652]} 55.866 {[ 23.3517]} 599.53 {[599.5325]} 55.586 {[ 25.8488]} 557.82 {[557.8191]} 54.967 {[ 28.7803]} 519.01 {[519.0081]} 53.781 {[ 32.3659]} 482.9 {[482.8973]} 51.479 {[ 37.2235]} 449.3 {[449.2990]}
% Vm
% Dm
% Tm
% fun = @(V) 0.5*V.*V.*rho.*S.*cd-T;
% x=fzero(fun,1);
See the documentation on interp1 to understand the intersections calculations.
.
  10 Kommentare
Alina Abdikadyr
Alina Abdikadyr am 27 Jan. 2023
Sorry, but in the table, there are Vmax and Vmin, how could I find a corresponding values of P for calculated Vmax and Vmin
The table gives Pmax and Pmin from a plot, but how to get P values for Vmax and Vmin?
Xmax Xmin Y1 Y2 Pmax Pmin Vmax Vmin
______ ____________ _____ ______________ _____ ______ ______ ______
56.013 {0×0 double} 25000 {0×0 double } 30859 9283.1 61.008 15.752
55.822 {0×0 double} 22198 {0×0 double } 33404 10049 66.04 17.051
55.235 {0×0 double} 19710 {0×0 double } 36159 10878 71.486 18.458
53.92 {0×0 double} 17501 {0×0 double } 39141 11775 77.382 19.98
51.119 {[ 23.8756]} 15539 {[1.5539e+04]} 42369 12746 83.764 21.628
39.638 {[ 39.6377]} 13798 {[1.3798e+04]} 45863 13797 90.673 23.412
Star Strider
Star Strider am 27 Jan. 2023
Use the returned indices to refer to them.
I calculated them as:
PmaxV(k1,:) = P(Vxidx);
PminV(k1,:) = P(Vnidx);
where ‘Vxidx’ and ‘Vnidx’ are returned from:
[maxV,Vxidx] = max(V);
[minV,Vnidx] = min(V);
And added them to the ‘Intersections’ table as well —
W = 10000;
S = 40;
AR = 7;
cd0 = 0.005;
k = 1 / pi / AR;
clalpha = 2*pi;
Psl=25000;
figure(1);hold on; xlabel('V');ylabel('P')
hv=0:1.6484:8.242;
for k1 = 1:numel(hv)
h = hv(k1);
i=0;
for alpha = 1:0.25:15
i=i+1;
rho(i)=1.225*exp(-h/10.4);
cl(i) = clalpha * alpha * pi/180;
V(i) = sqrt(2*W/rho(i)/S/cl(i));
L(i) = 0.5 * rho(i)* V(i) * V(i) * S * cl(i);
cd(i) = cd0 + k * cl(i) * cl(i);
D(i) = 0.5 * rho(i) * V(i) * V(i) * S * cd(i);
clcd(i) = cl(i)/cd(i);
P(i) = D(i)*V(i);
ang(i) = alpha;
Ph(i)=Psl*(rho(i)/1.225).^0.75;
end
Vm(k1,:) = V;
Pm(k1,:) = P;
Pmh(k1,:) = Ph;
xix = find(diff(sign(D.*V-Ph)));
if ~isempty(xix)
for k3 = 1:numel(xix)
idxrng = max(xix(k3)-1,1) : min(numel(V),xix(k3)+1);
xv(k3) = interp1(P(idxrng)-Ph(idxrng), V(idxrng), 0);
yv(k3) = interp1(V(idxrng),P(idxrng), xv(k3));
Xm{k1,k3} = xv(k3); % Intersection X-Matrix
Ym{k1,k3} = yv(k3); % Intersection Y-Matrix
end
end
figure(1)
plot(V,P)
hold on
plot(V,Ph)
xlim([0 90])
ylim([0 45000])
plot(xv, yv, 's', 'DisplayName',sprintf('Intersections %d',k1)) % Plot Intersections (Optional)
[maxP,idx] = max(P);
[minP,idx] = min(P);
Pmax(k1,:) = maxP;
Pmin(k1,:) = minP;
[maxV,Vxidx] = max(V);
[minV,Vnidx] = min(V);
Vmax(k1,:) = maxV;
Vmin(k1,:) = minV;
PmaxV(k1,:) = P(Vxidx);
PminV(k1,:) = P(Vnidx);
end
%legend('Location','eastoutside','NumColumns',2)
Intersections = cell2table(cat(2,Xm,Ym), 'VariableNames',{'Xmax','Xmin','Y1','Y2'});
Intersections = addvars(Intersections, Pmax,Pmin,Vmax,Vmin,PmaxV,PminV, 'After','Y2')
Intersections = 6×10 table
Xmax Xmin Y1 Y2 Pmax Pmin Vmax Vmin PmaxV PminV ______ ____________ _____ ______________ _____ ______ ______ ______ _____ _____ 56.013 {0×0 double} 25000 {0×0 double } 30859 9283.1 61.008 15.752 30859 12261 55.822 {0×0 double} 22198 {0×0 double } 33404 10049 66.04 17.051 33404 13273 55.235 {0×0 double} 19710 {0×0 double } 36159 10878 71.486 18.458 36159 14367 53.92 {0×0 double} 17501 {0×0 double } 39141 11775 77.382 19.98 39141 15552 51.119 {[ 23.8756]} 15539 {[1.5539e+04]} 42369 12746 83.764 21.628 42369 16835 39.638 {[ 39.6377]} 13798 {[1.3798e+04]} 45863 13797 90.673 23.412 45863 18223
% Xm1=cell2mat(Xm);
.

Melden Sie sich an, um zu kommentieren.

Weitere Antworten (1)

John D'Errico
John D'Errico am 23 Jan. 2023
Bearbeitet: John D'Errico am 23 Jan. 2023
You cannot use fzero on a vector. A vector of elements is NOT a function. It is just a list of numbers. And while you may think of it as a function, it is not. Essentially, a vector of numbers is just a picture of a function.
Suppose you were not feeling well today, so you decided to visit your doctor. But rather than appearing for an examination, you decided to just send your doctor a picture of yourself. Would that be sufficient to solve your problem? Rather than bring your car in for an oil change, would you just send a picture of your car to the JiffyLube place? Fzero needs a function, not just a list of numbers.
Instead, you will need to use tools to identify where the curves cross, perhaps using methods of interpolation. Or, better yet, you could use the actual functions themselves, now as true functions that fzero can evaluate at any point.

Kategorien

Mehr zu MATLAB finden Sie in Help Center und File Exchange

Produkte


Version

R2022b

Community Treasure Hunt

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

Start Hunting!

Translated by