I have an error where it says unrecognized function or variable x when using fsolve.

2022 1 19
1 个回答
回答已采纳
更新时间：2022 1 22
12 次查看（30 天）
Ran in:

0 个投票

Ran in:
I have a variable that is called x(2) and I need to solve this using fsolve. However, variables like rho and speedsound are dependent on this x(2) variable to be solved hence I have another function (atmos) below. Does anyone know how I can fix this error in my script?
clear all;
clc;
W = 1248.5*9.81;
S = 17.1;
list_aoa = linspace(0.0438914335842468,0.316137577406202,100);
h_V_aoa_eledefl = zeros(100,4);
% CL = (6.44*aoa + 0.355*eledefl);
% CD = (0.03 + 0.05*(6.44*aoa + 0.355*eledefl)^2);
% Cm = (0.05 - 0.683*aoa - 0.923*eledefl);
% thrust = (3 * ( (7 + V/speedsound )*200/3 + h*(2*(V/speedsound) - 11) ));
% V = x(1);
% height = x(2);
% eledefl = x(3);
% CL = (6.44*x(2) + 0.355*x(3));
% CD = (0.03 + 0.05*(6.44*x(2) + 0.355*x(3))^2);
% thrust = ( 3 * ( 200/3*( 7 + x(1)/speedsound ) + height/1000*( 2*x(1)/speedsound - 11) ) );
for i=1:1:100
func_1 = @(x) ( 0.5 * rho * x(1)^2 * S * (6.44*list_aoa(i) + 0.355*x(3)) + ( 3 * ( 200/3*( 7 + x(1)/speedsound ) + x(2)/1000*( 2*x(1)/speedsound - 11) ) )*sin(list_aoa(i)) - W ) ;
func_2 = @(x) ( 0.5 * rho * x(1)^2 * S * (0.03 + 0.05*(6.44*list_aoa(i) + 0.355*x(3))^2) - ( 3 * ( 200/3*( 7 + x(1)/speedsound ) + x(2)/1000*( 2*x(1)/speedsound - 11) ) )*cos(list_aoa(i)) );
func_3 = @(x) ( 0.05 - 0.683*list_aoa(i) - 0.923*x(3) );
[T, p, rho, speedsound] = atmos(x(2));
FUNC = @(x) [func_1(x); func_2(x); func_3(x)];
X = fsolve(@(x) FUNC(x), [10 0 -30/180*pi]);
h_V_aoa_eledefl(a,1) = X(2); % h
h_V_aoa_eledefl(a,2) = X(1); %V
h_V_aoa_eledefl(a,3) = list_aoa(i); %aoa
h_V_aoa_eledefl(a,4) = X(3); %eledefl
end
Unrecognized function or variable 'x'.
%% Atmos function
% Calculating T, p ,rho and speedsound for every altitude in the ISA atmosphere
function [T, p, rho, speedsound] = atmos(h)
h1 = 11000; % Height of tropopause
h2 = 20000; % End height of table
g = 9.81;
R = 287;
c = 6.51e-3; % temperature lapse dt/dh = - c = -6.51 degcelcius/km
T0 = 15+273.15; % Temperature sea level
p0 = 101325; % pressure sealevel
rho0 = 101325/R/T0; % density sealevel = pressure / R*T, R=287, T = 15 degcelcius
T1 = T0 - c*h1; % Temperature at 11km
p1 = p0 * (T1/T0)^5.2506; % Pressure at 11km
rho1 = rho0 * (T1/T0)^4.2506; % Density at 11km
T2 = T1; % Temperature at 20km
p2 = p1 * exp(-g/(R*T2)*(h2-h1)); % Pressure at 20km
rho2 = rho1 * exp(-g/(R*T2)*(h2-h1)); % Density at 20km
if h <= h1
% disp('Troposphere');
T = T0 - c*h;
p = p0 * (T/T0)^5.2506;
rho = rho0 * (T/T0)^4.2506;
speedsound = (1.4*R*T)^0.5;
elseif h <= h2
% disp('Tropopause');
T = T1;
p = p1 * exp(-g/(R*T)*(h-h1));
rho = rho1 * exp(-g/(R*T)*(h-h1));
speedsound = (1.4*R*T)^0.5;
end
return
end

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

 采纳的回答

Torsten
Torsten 2022-1-19

0 个投票

function main
W = 1248.5*9.81;
S = 17.1;
list_aoa = linspace(0.0438914335842468,0.316137577406202,100);
h_V_aoa_eledefl = zeros(100,4);
% CL = (6.44*aoa + 0.355*eledefl);
% CD = (0.03 + 0.05*(6.44*aoa + 0.355*eledefl)^2);
% Cm = (0.05 - 0.683*aoa - 0.923*eledefl);
% thrust = (3 * ( (7 + V/speedsound )*200/3 + h*(2*(V/speedsound) - 11) ));
% V = x(1);
% height = x(2);
% eledefl = x(3);
% CL = (6.44*x(2) + 0.355*x(3));
% CD = (0.03 + 0.05*(6.44*x(2) + 0.355*x(3))^2);
% thrust = ( 3 * ( 200/3*( 7 + x(1)/speedsound ) + height/1000*( 2*x(1)/speedsound - 11) ) );
for i=1:1:100
X = fsolve(@(x) FUNC(x,W,S,list_aoa,i), [10 0 -30/180*pi]);
h_V_aoa_eledefl(a,1) = X(2); % h
h_V_aoa_eledefl(a,2) = X(1); %V
h_V_aoa_eledefl(a,3) = list_aoa(i); %aoa
h_V_aoa_eledefl(a,4) = X(3); %eledefl
end
end
function RES = FUNC(x,W,S,list_aoa,i)
[T, p, rho, speedsound] = atmos(x(2));
RES(1) = ( 0.5 * rho * x(1)^2 * S * (6.44*list_aoa(i) + 0.355*x(3)) + ( 3 * ( 200/3*( 7 + x(1)/speedsound ) + x(2)/1000*( 2*x(1)/speedsound - 11) ) )*sin(list_aoa(i)) - W ) ;
RES(2) = ( 0.5 * rho * x(1)^2 * S * (0.03 + 0.05*(6.44*list_aoa(i) + 0.355*x(3))^2) - ( 3 * ( 200/3*( 7 + x(1)/speedsound ) + x(2)/1000*( 2*x(1)/speedsound - 11) ) )*cos(list_aoa(i)) );
RES(3) = ( 0.05 - 0.683*list_aoa(i) - 0.923*x(3) );
end
%% Atmos function
% Calculating T, p ,rho and speedsound for every altitude in the ISA atmosphere
function [T, p, rho, speedsound] = atmos(h)
h1 = 11000; % Height of tropopause
h2 = 20000; % End height of table
g = 9.81;
R = 287;
c = 6.51e-3; % temperature lapse dt/dh = - c = -6.51 degcelcius/km
T0 = 15+273.15; % Temperature sea level
p0 = 101325; % pressure sealevel
rho0 = 101325/R/T0; % density sealevel = pressure / R*T, R=287, T = 15 degcelcius
T1 = T0 - c*h1; % Temperature at 11km
p1 = p0 * (T1/T0)^5.2506; % Pressure at 11km
rho1 = rho0 * (T1/T0)^4.2506; % Density at 11km
T2 = T1; % Temperature at 20km
p2 = p1 * exp(-g/(R*T2)*(h2-h1)); % Pressure at 20km
rho2 = rho1 * exp(-g/(R*T2)*(h2-h1)); % Density at 20km
if h > h2
disp('Error occured.')
end
if h <= h1
% disp('Troposphere');
T = T0 - c*h;
p = p0 * (T/T0)^5.2506;
rho = rho0 * (T/T0)^4.2506;
speedsound = (1.4*R*T)^0.5;
elseif h <= h2
% disp('Tropopause');
T = T1;
p = p1 * exp(-g/(R*T)*(h-h1));
rho = rho1 * exp(-g/(R*T)*(h-h1));
speedsound = (1.4*R*T)^0.5;
end
return
end

7 个评论
显示 5更早的评论 隐藏 5更早的评论

Rachel Ong
Rachel Ong 2022-1-20
Hi, thanks for your help but I think there is still an error and it cant run. What does RES do?
Torsten
Torsten 2022-1-20
编辑：Torsten 2022-1-20
The error is that you don't supply a formula for T, p, rho, speedsound if x(2) > h2.
And this seems to be the case during computation.
RES(1),...,RES(3) is the residual that you return for FSOLVE. It's the same as if you write
function RES = FUNC(x,W,S,list_aoa,i)
[T, p, rho, speedsound] = atmos(x(2));
func_1 = @(x) ( 0.5 * rho * x(1)^2 * S * (6.44*list_aoa(i) + 0.355*x(3)) + ( 3 * ( 200/3*( 7 + x(1)/speedsound ) + x(2)/1000*( 2*x(1)/speedsound - 11) ) )*sin(list_aoa(i)) - W ) ;
func_2 = @(x) ( 0.5 * rho * x(1)^2 * S * (0.03 + 0.05*(6.44*list_aoa(i) + 0.355*x(3))^2) - ( 3 * ( 200/3*( 7 + x(1)/speedsound ) + x(2)/1000*( 2*x(1)/speedsound - 11) ) )*cos(list_aoa(i)) );
func_3 = @(x) ( 0.05 - 0.683*list_aoa(i) - 0.923*x(3) );
RES = [func_1(x); func_2(x); func_3(x)];
end
But it's not necessary to use functions in this context.
Rachel Ong
Rachel Ong 2022-1-20
But for this, it shouldnt exceed the 20km mark of the x(2) > h2. How should I edit this to make it work?
Torsten
Torsten 2022-1-21
编辑：Torsten 2022-1-21
I don't know whether it's correct, but I changed the unit of x(2) in your equations (from km to m ?) and this gives a result.
function main
W = 1248.5*9.81;
S = 17.1;
list_aoa = linspace(0.0438914335842468,0.316137577406202,100);
h_V_aoa_eledefl = zeros(100,4);
% CL = (6.44*aoa + 0.355*eledefl);
% CD = (0.03 + 0.05*(6.44*aoa + 0.355*eledefl)^2);
% Cm = (0.05 - 0.683*aoa - 0.923*eledefl);
% thrust = (3 * ( (7 + V/speedsound )*200/3 + h*(2*(V/speedsound) - 11) ));
% V = x(1);
% height = x(2);
% eledefl = x(3);
% CL = (6.44*x(2) + 0.355*x(3));
% CD = (0.03 + 0.05*(6.44*x(2) + 0.355*x(3))^2);
% thrust = ( 3 * ( 200/3*( 7 + x(1)/speedsound ) + height/1000*( 2*x(1)/speedsound - 11) ) );
X0 = [60 2 0.02];
for i=1:1:100
X = fsolve(@(x) FUNC(x,W,S,list_aoa,i), X0)
h_V_aoa_eledefl(i,1) = X(2); % h
h_V_aoa_eledefl(i,2) = X(1); %V
h_V_aoa_eledefl(i,3) = list_aoa(i); %aoa
h_V_aoa_eledefl(i,4) = X(3); %eledefl
X0 = X;
end
figure(1)
plot(h_V_aoa_eledefl(:,3),h_V_aoa_eledefl(:,1))
figure(2)
plot(h_V_aoa_eledefl(:,3),h_V_aoa_eledefl(:,2))
figure(3)
plot(h_V_aoa_eledefl(:,3),h_V_aoa_eledefl(:,4))
end
function RES = FUNC(x,W,S,list_aoa,i)
[T, p, rho, speedsound] = atmos(x(2));
RES(1) = ( 0.5 * rho * x(1)^2 * S * (6.44*list_aoa(i) + 0.355*x(3)) + ( 3 * ( 200/3*( 7 + x(1)/speedsound ) + x(2)*( 2*x(1)/speedsound - 11) ) )*sin(list_aoa(i)) - W ) ;
RES(2) = ( 0.5 * rho * x(1)^2 * S * (0.03 + 0.05*(6.44*list_aoa(i) + 0.355*x(3))^2) - ( 3 * ( 200/3*( 7 + x(1)/speedsound ) + x(2)*( 2*x(1)/speedsound - 11) ) )*cos(list_aoa(i)) );
RES(3) = ( 0.05 - 0.683*list_aoa(i) - 0.923*x(3) );
end
%% Atmos function
% Calculating T, p ,rho and speedsound for every altitude in the ISA atmosphere
function [T, p, rho, speedsound] = atmos(h)
h1 = 11000; % Height of tropopause
h2 = 20000; % End height of table
g = 9.81;
R = 287;
c = 6.51e-3; % temperature lapse dt/dh = - c = -6.51 degcelcius/km
T0 = 15+273.15; % Temperature sea level
p0 = 101325; % pressure sealevel
rho0 = 101325/R/T0; % density sealevel = pressure / R*T, R=287, T = 15 degcelcius
T1 = T0 - c*h1; % Temperature at 11km
p1 = p0 * (T1/T0)^5.2506; % Pressure at 11km
rho1 = rho0 * (T1/T0)^4.2506; % Density at 11km
T2 = T1; % Temperature at 20km
p2 = p1 * exp(-g/(R*T2)*(h2-h1)); % Pressure at 20km
rho2 = rho1 * exp(-g/(R*T2)*(h2-h1)); % Density at 20km
if h > h2
disp('Error occured.')
end
if h <= h1
% disp('Troposphere');
T = T0 - c*h;
p = p0 * (T/T0)^5.2506;
rho = rho0 * (T/T0)^4.2506;
speedsound = (1.4*R*T)^0.5;
elseif h <= h2
% disp('Tropopause');
T = T1;
p = p1 * exp(-g/(R*T)*(h-h1));
rho = rho1 * exp(-g/(R*T)*(h-h1));
speedsound = (1.4*R*T)^0.5;
end
return
end
Rachel Ong
Rachel Ong 2022-1-22
I see, the results are not right but I think at least I get how to solve it now. I will try to edit it. However, can I ask why the workspace is empty after running? I thought that it would be saved under h_V_aoa_eledefl but the workspace is just empty. Thanks for your help though really appreciate it!
Torsten
Torsten 2022-1-22
However, can I ask why the workspace is empty after running?
I think you will have to remove the first line
function main
and the
end
at the end of the function.

请先登录，再进行评论。

更多回答（0 个）

类别

在 帮助中心 和 File Exchange 中查找有关 Annotations 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by