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

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

Answer 1

[Torsten]

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

[Rachel Ong]

Hi, thanks for your help but I think there is still an error and it cant run. What does RES do?

[Torsten]

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]

But for this, it shouldnt exceed the 20km mark of the x(2) > h2. How should I edit this to make it work?

[Torsten]

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]

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]

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.

[Rachel Ong]

I see. Thank you for your help!