MATLAB: Area Mach Number Relation

mach numberMATLABnumerical integration

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?

Here my input ARatio is: ARatio = 1.5

and i get two solutions: Msub & Msup

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)

I need to get Msub & Msup for an array of ARatios from 0.1 to 10, in order to then plot ARatio and the Mach number solutions

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');

Best Answer

The easiest way is to just enclose the essence of your calculations in a loop:

% Solve for Msub and Msup using this area ratio (A/A*)
ARatio = 0:0.25:2.5;
% 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;
for k = 1:numel(ARatio)
    % 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(k));
            fjp1 = AM_EQN(M+dM,ARatio(k));
            
            % 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(k));
            fjp1 = AM_EQN(M+dM,ARatio(k));
            
            % 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');
end

producing:

==== INCREMENTAL SEARCH SOLVER ====
Msub: 1000.0000 after 0 iterations
Msup: 10001.0000 after 0 iterations
===================================
==== INCREMENTAL SEARCH SOLVER ====
Msub: 1000.0000 after 0 iterations
Msup: 10001.0000 after 0 iterations
===================================
==== INCREMENTAL SEARCH SOLVER ====
Msub: 1000.0000 after 0 iterations
Msup: 10001.0000 after 0 iterations
===================================
==== INCREMENTAL SEARCH SOLVER ====
Msub: 1000.0000 after 0 iterations
Msup: 10001.0000 after 0 iterations
===================================
==== INCREMENTAL SEARCH SOLVER ====
Msub: 1000.0000 after 0 iterations
Msup: 10001.0000 after 0 iterations
===================================
==== INCREMENTAL SEARCH SOLVER ====
Msub: 0.5533 after 5 iterations
Msup: 1.5997 after 5 iterations
===================================
==== INCREMENTAL SEARCH SOLVER ====
Msub: 0.4303 after 5 iterations
Msup: 1.8541 after 6 iterations
===================================
==== INCREMENTAL SEARCH SOLVER ====
Msub: 0.3566 after 5 iterations
Msup: 2.0433 after 6 iterations
===================================
==== INCREMENTAL SEARCH SOLVER ====
Msub: 0.3059 after 5 iterations
Msup: 2.1972 after 6 iterations
===================================
==== INCREMENTAL SEARCH SOLVER ====
Msub: 0.2685 after 5 iterations
Msup: 2.3282 after 6 iterations
===================================
==== INCREMENTAL SEARCH SOLVER ====
Msub: 0.2395 after 5 iterations
Msup: 2.4428 after 6 iterations
===================================

It would likely be easier to create a function from your code, then call it with different inputs, producing the desired output.

Make necessary changes to create the result you want.

Best Answer

Related Solutions

MATLAB: Stopping a function after 30 iterations

Related Question