Area Mach Number Relation

The following code shows me how to get converged Mach number solutions for both subsonic and supersonic given Area ratio (ARatio), however, how can i input a range of ARatio instead of just one value so that the solutions are an array of converged mach numbers?

clear;

clc;

%% INPUTS

% Define some paramters

g = 1.4;

gm1 = g-1;

gp1 = g+1;

% Define anonymous function with two inputs (M and ARatio)

% - Will be used in the methods below

% - Pass M and ARatio as arguments to AM_EQN to get function value

% - funVal = AM_EQN(M,ARatio)

AM_EQN = @(M,ARatio) sqrt((1/M^2)*(((2+gm1*M^2)/gp1)^...

(gp1/gm1)))-ARatio;

% Solve for Msub and Msup using this area ratio (A/A*)

ARatio = 1.5;

% Error tolerance

errTol = 1e-4;

% Flags for printing iterations to screen

verboseBisection = 0;

verboseIncremental = 0;

%% SUBSONIC INCREMENTAL SEARCH

% Initial values

dM = 0.1; % Initial M step size

M = 1e-6; % Initial M value

iConvSub = 0; % Initial converge index

if (verboseIncremental == 1)

fprintf('Incremental Search Method: Subsonic\n');

fprintf('-----------------------------------\n');

end

% Iterate to solve for root

iterMax = 100; % Maximum iterations

stepMax = 100; % Maximum step iterations

for i = 1:1:iterMax

for j = 1:1:stepMax

% Evaluate function at j and j+1

fj = AM_EQN(M,ARatio);

fjp1 = AM_EQN(M+dM,ARatio);

% Print iterations to command window

if (verboseIncremental == 1)

fprintf('fj | fjp1: %3.4f\t%3.4f\n',fj,fjp1);

end

% Update M depending on sign change or not

% - If no sign change, then we are not bounding root yet

% - If sign change, then we are bounding the root, update dM

if (fj*fjp1 > 0)

M = M + dM; % Update M

elseif (fj*fjp1 < 0)

dM = dM*0.1; % Refine the M increment

break; % Break out of j loop

end

end % END: j Loop

% Check for convergence

if (abs(fj-fjp1) <= errTol) % If converged

iConvSub = i; % Set converged index

break; % Exit loop

end

end % END: i Loop

% Set subsonic Mach number to final M from iterations

Msub = M;

%% SUPERSONIC INCREMENTAL SEARCH

% Initial values

dM = 1; % Initial M step size

M = 1+1e-6; % Initial M value

iConvSup = 0; % Initial converge index

if (verboseIncremental == 1)

fprintf('\nIncremental Search Method: Supersonic\n');

fprintf('-------------------------------------\n');

end

% Iterate to solve for root

iterMax = 100; % Maximum iterations

stepMax = 100; % Maximum step iterations

for i = 1:1:iterMax

for j = 1:1:stepMax

% Evaluate function at j and j+1

fj = AM_EQN(M,ARatio);

fjp1 = AM_EQN(M+dM,ARatio);

% Print iterations to command window

if (verboseIncremental == 1)

fprintf('fj | fjp1: %3.4f\t%3.4f\n',fj,fjp1);

end

% Update M depending on sign change or not

% - If no sign change, then we are not bounding root yet

% - If sign change, then we are bounding the root, update dM

if (fj*fjp1 > 0)

M = M + dM; % Update M

elseif (fj*fjp1 < 0)

dM = dM*0.1; % Refine the M increment

break; % Break out of j loop

end

end % END: j Loop

% Check for convergence

if (abs(fj-fjp1) <= errTol) % If converged

iConvSup = i; % Set converged index

break; % Exit loop

end

end % END: i Loop

% Set supersonic Mach number to final M from iterations

Msup = M;

% Print solutions to command window

fprintf('==== INCREMENTAL SEARCH SOLVER ====\n');

fprintf('Msub: %3.4f after %i iterations\n',Msub,iConvSub);

fprintf('Msup: %3.4f after %i iterations\n',Msup,iConvSup);

fprintf('===================================\n\n');

Answer 1

[Chris]

You can format posts to look llike ML text and make them easier to read. Without reading any of your code it sounds like you can write a wraper function if you can not eaisly vectorize the code.

[Steven Castrillon]

Hopefully this is better formatting. Here my input ARatio is: ARatio = 1.5

However if i input ARatio as: ARatio = (0.1:0.1:10) i get the following error

Error in AREAMACH2 (line 62)

if (fj*fjp1 > 0)

If this can be solved by using a wraper function could you please help?

clear;
clc;
%% INPUTS
% Define some paramters
g   = 1.4;
gm1 = g-1;
gp1 = g+1;
% Define anonymous function with two inputs (M and ARatio)
% - Will be used in the methods below
% - Pass M and ARatio as arguments to AM_EQN to get function value
% -     funVal = AM_EQN(M,ARatio)
AM_EQN = @(M,ARatio) sqrt((1/M^2)*(((2+gm1*M^2)/gp1)^...
                            (gp1/gm1)))-ARatio;
% Solve for Msub and Msup using this area ratio (A/A*)
ARatio = 1.5;
% Error tolerance
errTol = 1e-4;
% Flags for printing iterations to screen
verboseBisection   = 0;
verboseIncremental = 0;
%% SUBSONIC INCREMENTAL SEARCH
% Initial values
dM       = 0.1;                                                             % Initial M step size
M        = 1e-6;                                                            % Initial M value
iConvSub = 0;                                                               % Initial converge index
if (verboseIncremental == 1)
    fprintf('Incremental Search Method: Subsonic\n');
    fprintf('-----------------------------------\n');
end
% Iterate to solve for root
iterMax = 100;                                                              % Maximum iterations
stepMax = 100;                                                              % Maximum step iterations
for i = 1:1:iterMax
    for j = 1:1:stepMax
        
        % Evaluate function at j and j+1
        fj   = AM_EQN(M,ARatio);
        fjp1 = AM_EQN(M+dM,ARatio);
        
        % Print iterations to command window
        if (verboseIncremental == 1)
            fprintf('fj | fjp1: %3.4f\t%3.4f\n',fj,fjp1);
        end
        
        % Update M depending on sign change or not
        % - If no sign change, then we are not bounding root yet
        % - If sign change, then we are bounding the root, update dM
        if (fj*fjp1 > 0)
            M = M + dM;                                                     % Update M
        elseif (fj*fjp1 < 0)
            dM = dM*0.1;                                                    % Refine the M increment
            break;                                                          % Break out of j loop
        end
        
    end % END: j Loop
    
    % Check for convergence
    if (abs(fj-fjp1) <= errTol)                                             % If converged
        iConvSub = i;                                                       % Set converged index
        break;                                                              % Exit loop
    end
    
end % END: i Loop
% Set subsonic Mach number to final M from iterations
Msub = M;
%% SUPERSONIC INCREMENTAL SEARCH
% Initial values
dM       = 1;                                                               % Initial M step size
M        = 1+1e-6;                                                          % Initial M value
iConvSup = 0;                                                               % Initial converge index
if (verboseIncremental == 1)
    fprintf('\nIncremental Search Method: Supersonic\n');
    fprintf('-------------------------------------\n');
end
% Iterate to solve for root
iterMax = 100;                                                              % Maximum iterations
stepMax = 100;                                                              % Maximum step iterations
for i = 1:1:iterMax
    for j = 1:1:stepMax
        
        % Evaluate function at j and j+1
        fj   = AM_EQN(M,ARatio);
        fjp1 = AM_EQN(M+dM,ARatio);
        
        % Print iterations to command window
        if (verboseIncremental == 1)
            fprintf('fj | fjp1: %3.4f\t%3.4f\n',fj,fjp1);
        end
        
        % Update M depending on sign change or not
        % - If no sign change, then we are not bounding root yet
        % - If sign change, then we are bounding the root, update dM
        if (fj*fjp1 > 0)
            M = M + dM;                                                     % Update M
        elseif (fj*fjp1 < 0)
            dM = dM*0.1;                                                    % Refine the M increment
            break;                                                          % Break out of j loop
        end
        
    end % END: j Loop
    
    % Check for convergence
    if (abs(fj-fjp1) <= errTol)                                             % If converged
        iConvSup = i;                                                       % Set converged index
        break;                                                              % Exit loop
    end
    
end % END: i Loop
% Set supersonic Mach number to final M from iterations
Msup = M;
% Print solutions to command window
fprintf('==== INCREMENTAL SEARCH SOLVER ====\n');
fprintf('Msub: %3.4f after %i iterations\n',Msub,iConvSub);
fprintf('Msup: %3.4f after %i iterations\n',Msup,iConvSup);
fprintf('===================================\n\n');

Answer 2

[Chris]

Put all that code info a function with inputs/outputs; call it for each ARatio you want. https://en.wikipedia.org/wiki/Wrapper_function

Control statements with vectors are commonly inapropreate.