how to plot a certain level in a contour?

I 'm using contour option in Matlab R2014a to plot many curves: contour(x,y,c,cilevels); with cilevels = [0.0,.005,0.019];

but I need to select one curve (level) from the contour: the option peaks is not working? and I have tried : contour(x,y,c,[1 1]); but it was error also?

could you please help me?

Answer 1

[Star Strider]

Extracting the contours for a specific level is straightforward.

Try this —

[x,y,c] = peaks(20);
cilevels = [0.0,.005,0.019];

figure
[c,h] = contour(x,y,c,cilevels);

axlim = axis;

Level = 0.005;                                          % Choose Level
startIdx = find(c(1,:) == Level);
vlen = c(2,startIdx);

for k = 1:numel(startIdx)
    % Q = [startIdx(k); vlen(k); startIdx(k)+vlen(k)]
    idxrng = startIdx(k)+1 : startIdx(k)+vlen(k);
    xc{k} = c(1,idxrng);
    yc{k} = c(2,idxrng);
end

figure
hold on
for k = 1:numel(startIdx)
    plot(xc{k}, yc{k})
end
hold off
axis(axlim)

.

[Mohamed Hajjaj]

I'm using R2014a , may be this make the error:

Error in csorrsommerfeldstabilitytemp (line 156)

[t_Res,t_alphas,imag(t_c)] = peaks(20);

[Walter Roberson]

[t_Res,t_alphas,imag(t_c)] = peaks(20);

That tries to assign the third output to an array named imag indexed by t_c

Chances are that t_c does not happen entirely consist of positive integers, if it exists at all. Chances are that t_c does not happen to have 400 (20 * 20) elements.

Note that the code does not assign to the imaginary component of t_c. Note too that if the assignment did somehow suceed, that afterwards you would be unable to call the imag function.

[Mohamed Hajjaj]

t_c already exist and has a large number of values, it is a complex number, so its imaginary part is a real number. If imag(t_c) is a real and cilevels = [0.0,.005, 0.0019];

so does the suggested solution fail?

[Star Strider]

@Walter Roberson — Thank you.

@Mohamed Hajjaj — I no longer have access to R2014a. (That was the last release prior to a major change that also introduced ‘Handle Graphics 2’).

I would have to have whatever ‘t_c’ is in order to se if I can simulate your code. It might be possible to create the result you want using scatteredInterpolant (introduced in R2013a), however I would have to have your data to determine that.

I also do not entirely understand the reason you want to use the peaks function to create that contour. Please descrribe what you actually want to do, and what problem you are solving.

.

[Mohamed Hajjaj]

It works with

p = peaks(2);

contour(t_Res,t_alphas,imag(t_c),[p p]);

Thank you very much for your help.

If I need to save the data (points of Res, alpha, imag c) of selected curve (level) in a file, what the option I use?

[Star Strider]

My pleasure!

Without actually having yur data to work with, craeazting it as a cell array is the best option, and then using the writecell function. Another option, of course, is to save it as a .mat file. It depends on what you want to do with it later.

I still do not understand what you are doing. If I did, I could probably provide a more exact reply.

If my Answer helped you solve your problem, please Accept it!

.

[Mohamed Hajjaj]

I'm working in the field of computational fluid dynamics, the stability of fluid flow. Plotting the curves is important and also reporting certain data (critical values), so I need the numeric values for certain curves to validate my results.

[Star Strider]

I need the actual data and requirements, or some sort of analogous problem that can be directly applied to what you want to do here, in order to proceed with this. I cannot guess what you want to do.

[Mohamed Hajjaj]

Ok, the analogous problem in the below link where the m file is available:

https://www.cdsimpson.net/2015/04/temporal-and-spatial-stability-analysis.html

[Star Strider]

I edited your last comment to make the link live.

I’m still not sure what you want to do. If you want to create the contours in this plot image:

figure
imshow(imread('https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiA3lvtI3-S69hEb-g9XeGhmjcbv4GChuCG7__Yg4VaUxYXlpZJXRvDLembtZkARscC3EbNlL5Yvsx6yZwOkZV7sYONYShxHro147aqPT0CcxcLQFuI7jDH0ECBODvbcFDxpkzxUB9-VcLW/s1600/temporal_stability.png'))

you will have to calculate those contours. You can integrate the differential equations using any of the appropriate functions for that available in core MATLAB. I am not certain which one would be best, however the usual pne to begin with is ode45. If the parameters vary over a large range of magnitudes, creating a ‘stiff’ system, then ode15s would probably be best to start with. The derivatives apparently are with respect to time, and the contours are most likely different integrated results plotted against each other for different values of a particular parameter. (Comnputational fluid dynamics is not among my areas of expertise.)

I copied the code, and while it appears to run correctly (I substituted ode15s for ‘cs_rk4’ because the provided function was taking too long), it’s throwing an error that I don’t understand:

Index exceeds the number of array elements. Index must not exceed 3.

I added a plot call to see the results of the integration that you can remove.

It’s not obvious to me what is being counted, how to reset that counter, or how to reset the iteration limiits of the loop is throwing that error to match the actual number of whatever is being created earlier, whatever it is.

I leave troubleshooting this to you —

% Orr-Sommerfeld Stability Analysis
% 
% This script computes the Blasius velocity profile using a 4th order
% Runge-Kutta shooting method.  A finite difference method is then used to
% discretize the Orr-Sommerfeld equation, and the eigenvalue problem is
% solved using a real alpha value.  The process is repeated with a complex
% alpha, and a map is created to relate alpha and omega. A solution of the
% zeros of omega in the complex plane is then found to plot the spatial
% stability curves. There are some differences in the plots output by this
% code and those found in the report. The plots were labeled, scaled, and
% made monochromatic by hand in post-processing.
% 
% Programmed by:
% Christopher Simpson
% 11 March 2015
% cdsimpson.net

% clear all;
% close all;
% clc;

%% INPUTS
h = 0.1;   % eta step size - use this to control accuracy of calculations

% Velocity Profile Inputs

% % Use this for Poisseuille Flow, modify for other types
% xi = linspace(0,1,50);
% U = 1 - xi.^2;
% Upp = -2*ones(1,length(xi));
use_Blasius = 1; % use Blasius velocity profile below
transform_u = 1; % transform velocity to xi coordinates

% Temporal Stability Inputs
run_temporal_analysis = 1 ; % 1 = run temporal analysis, 0 = skip
t_Re_min = 1e3; % minimum Re
t_Re_max = 1e5; % maximum Re
t_N_Re = 40;    % number of Re
t_Res = logspace(log10(t_Re_min), log10(t_Re_max), t_N_Re);

t_alpha_min = 0.001;    % minimum alpha for temporal analysis
t_alpha_max = 1.3;      % maximum alpha
t_N_alpha = 40;         % number of alpha to use
t_alphas = linspace(t_alpha_min,t_alpha_max,t_N_alpha);

% Spatial Stability Inputs
run_spatial_analysis  = 1 ; % 1 = run spatial analysis, 0 = skip
s_Re_min = 1e3; % minimum Re
s_Re_max = 1e5; % maximum Re
s_N_Re = 100;   % number of Re - this needs to be large for resolution
s_Res = logspace(log10(s_Re_min), log10(s_Re_max), s_N_Re);

s_alphar_min = 0.000000001; % min real component of alpha
s_alphar_max = 2;   % max real component of alpha
s_N_alphar = 200;   % number of real alphas - large for good resolution
s_alphars = linspace(s_alphar_min,s_alphar_max,s_N_alphar);

% desired contour levels
s_aibyRelevels = [ 0, -1e-6, -2e-6, -4e-6, -6e-6, -7e-6 ];
s_alphais = zeros(1,length(s_aibyRelevels));
s_alphas = zeros(length(s_alphais),length(s_alphars));
for m1 = 1:length(s_alphais)
    for m2 = 1:length(s_alphars)
        alpha = s_alphars(m2)+1i*s_alphais(m1);
        s_alphas(m1,m2) = alpha;
    end
end


%% DEFINED FIGURE NUMBERS (do not correspond to report)
FIG_VEL_PROFILE = 1;
FIG_TEM_EIG = 2;
FIG_TEM_CONT = 3;
FIG_SPA_EIG = 4;
FIG_SPA_MAP = 5;
FIG_SPA_CONT = 6;

%% COMPUTE BLASIUS VELOCITY PROFILE
if ( use_Blasius )
    fprintf('Solving for Blasius Velocity Profile ... ');
    tic; % start counting elapsed time

    % Blasius ODE Definition:  f''' + ff'' = 0  ->  f''' = -ff''
    fun = @(eta,f) [ ...
        f(2); ...
        f(3); ...
        -f(1)*f(3); ...
        ];

    eta0 = 0;   % Initial eta station (@ wall )
    etaf = 10;  % Final eta station (@ "freestream")
    target = 1; % Freestream BC: f'(inf) = U(inf) = 1
    tol = 1e-6; % Shooting method tolerance

    f0 = [ 0; 0; .45 ]; % Initial condition for f, f', f'' at wall

    go = 1;
    while ( go == 1 ) % loop until converged
        % Standard 4th-order Runge-Kutta method
        % [ eta, f ] = cs_rk4( fun, eta0, etaf, f0, h );
        [eta, f] = ode15s(fun, [eta0 etaf], f0);
        
        figure                                          % <— PLOT CALL ADDED
        plot(eta, f)
        grid
        xlabel('\eta')
        ylabel('f')

        % compute error, adjust f''(0)
        err = f(2,end) - target;
        if ( abs(err) < tol )
            go = 0;
        else
            f0(3,1) = f0(3,1) * ( target-err );
        end
    end
    U = f(2,:);             % U   = f'(eta)
    Up = f(3,:);            % U'  = f''(eta)
    Upp = -f(1,:).*f(3,:);  % U'' = f'''(eta)
    elapsed_time = toc;     % compute time elapsed
    fprintf('Done (%f sec)\n',elapsed_time);
end

Solving for Blasius Velocity Profile ... 

Done (0.443772 sec)

%% TRANSFORM ETA TO XI
if ( transform_u )
    fprintf('Transforming Coordinates ... ');
    tic; % start counting elapsed time

    % find delta in terms of eta
    minval = 1; % initialize minimization variable
    index = 0;  % index of the edge of the Boundary Layer
    delta_vel = 0.999;  % percent of freestream velocity

    % find edge of Boundary layer, where U == delta_vel
    for n = 1:length(U)
        if ( abs(U(n)-delta_vel) < minval )
            index = n;
            minval = abs(U(n)-delta_vel);
        end
    end

    delta = eta(index); % Boundary layer thickness
    xi = eta/delta;     % define xi s.t. xi = 1 where eta = delta

    % transform to xi coordinate system
    % U = U;            % U(xi)   = U(eta)
    Up = Up*delta;      % U'(xi)  = U'(eta)*delta
    Upp = Upp*delta^2;  % U''(xi) = U''(eta)*delta^2

    elapsed_time = toc; % compute elapsed time
    fprintf('Done (%f sec)\n',elapsed_time);
end

Transforming Coordinates ... 

Done (0.001105 sec)

NMAX = length(xi);  % look at all eigenvalues
% NMAX = index;       % only look at eigenvalues in B.L.

%% TEMPORAL STABILITY ANALYSIS
if ( run_temporal_analysis )
    fprintf('Temporal Stability Analysis ...\n');
    tic; % start counting elapsed time

    h = xi(2)-xi(1);    % xi step size
    for l = 1:length(t_Res)
        Re = t_Res(l);
        for m = 1:length(t_alphas)
            alpha = t_alphas(m);
            fprintf('\ta = %.2f, Re = %d ... ',alpha,Re);

            % generate the A matrix
            A = zeros(NMAX); % initialize A

            n = 2; % start at n = 2, go to n = nmax-1

            % coefficients as defined in my paper
            a1 = -U(n)*alpha^3 - Upp(n)*alpha + 1i*alpha^4/Re;
            a2 = U(n)*alpha - 2*1i*alpha^2/Re;
            a3 = 1i/Re;

            % 2nd order Central Difference Method
            A(n,n-1) = a2 - 4*a3/h^2;
            A(n,n)   = a1*h^2 - 2*a2 + 6*a3/h^2;
            A(n,n+1) = a2 - 4*a3/h^2;
            A(n,n+2) = a3/h^2;
            % Values used in reference, not really central difference
            %   A(n,n-1) = a2-a3;
            %   A(n,n)   = a1*h^2 - 2*a2 + 3*a3;
            %   A(n,n+1) = a2-a3;
            %   A(n,n+2) = a3;

            for n = 3:NMAX-2
                % coefficients as defined in my paper
                a1 = -U(n)*alpha^3 - Upp(n)*alpha + 1i*alpha^4/Re;
                a2 = U(n)*alpha - 2*1i*alpha^2/Re;
                a3 = 1i/Re;
                % 2nd order Central Difference Method
                A(n,n-2) = a3/h^2;
                A(n,n-1) = a2 - 4*a3/h^2;
                A(n,n)   = a1*h^2 - 2*a2 + 6*a3/h^2;
                A(n,n+1) = a2 - 4*a3/h^2;
                A(n,n+2) = a3/h^2;
            end

            n = NMAX-1;
            a1 = -U(n)*alpha^3 - Upp(n)*alpha + 1i*alpha^4/Re;
            a2 = U(n)*alpha - 2*1i*alpha^2/Re;
            a3 = 1i/Re;
            % 2nd order Central Difference Method
            A(n,n-2) = a3/h^2;
            A(n,n-1) = a2 - 4*a3/h^2;
            A(n,n)   = a1*h^2 - 2*a2 + 6*a3/h^2;
            A(n,n+1) = a2 - 4*a3/h^2;

            A = A(2:end-1,2:end-1); % remove first/last rows/cols of A (0 BCs)

            % generate B
            B = zeros(NMAX);    % initialize B
            for n = 2:NMAX-1
                % coefficients as defined in my paper
                b1 = -alpha^3;
                b2 = alpha;
                % 2nd order Central Difference Method
                B(n,n-1) = b2;
                B(n,n)   = b1*h^2 - 2*b2;
                B(n,n+1) = b2;
            end
            B = B(2:end-1,2:end-1); % remove first/last rows/cols of B (0 BCs)


            [V,e] = eig(B\A);   % invert B, compute eigenvalues of LHS
            e = diag(e);    % turn diagonal eigenvalue matrix into vector

            % store most unstable eigenvalue
            [maxe,cindex] = max(imag(e));
            fprintf('c = %f + %fi\n',real(e(cindex)),imag(e(cindex)));
            t_c(m,l) = e(cindex);

            % Plot eigenvalues of the discrete system
            figure(FIG_TEM_EIG)
            plot(real(e),imag(e),'.','MarkerSize',6)
            xlabel('c_r')
            ylabel('c_i')
            title(sprintf('alpha = %.2f, Re = %d',alpha,Re));
            axis tight

        end
    end
    elapsed_time = toc; % compute elapsed time
    fprintf('Done (%f sec)\n',elapsed_time);
end

Temporal Stability Analysis ...

	a = 0.00, Re = 1000 ... 

Index exceeds the number of array elements. Index must not exceed 3.

%% SPATIAL STABILITY ANALYSIS
if ( run_spatial_analysis )
    fprintf('Spatial Stability Analysis ...\n');
    tic

    h = xi(2)-xi(1);    % xi step size
    
    num_pts = zeros(1,length(s_alphais));
    re_plot = zeros(length(s_alphais),1);
    om_plot = re_plot;
    for l = 1:length(s_Res)
        Re = s_Res(l);
        oldzerocount = 0;
        for m1 = 1:length(s_alphais)
            for m2 = 1:length(s_alphars)
                s_alphas(m1,m2) = s_alphars(m2) + 1i*s_aibyRelevels(m1)*Re;
                alpha = s_alphas(m1,m2);
                fprintf('\ta = %.2f+%.2fi, Re = %d ... ',real(alpha),imag(alpha),Re);

                % generate the A matrix
                A = zeros(NMAX); % initialize A

                n = 2; % start at n = 2, go to n = nmax-1

                % coefficients as defined in my paper
                a1 = -U(n)*alpha^3 - Upp(n)*alpha + 1i*alpha^4/Re;
                a2 = U(n)*alpha - 2*1i*alpha^2/Re;
                a3 = 1i/Re;

                % 2nd order Central Difference Method
                A(n,n-1) = a2 - 4*a3/h^2;
                A(n,n)   = a1*h^2 - 2*a2 + 6*a3/h^2;
                A(n,n+1) = a2 - 4*a3/h^2;
                A(n,n+2) = a3/h^2;

                for n = 3:NMAX-2
                    % coefficients as defined in my paper
                    a1 = -U(n)*alpha^3 - Upp(n)*alpha + 1i*alpha^4/Re;
                    a2 = U(n)*alpha - 2*1i*alpha^2/Re;
                    a3 = 1i/Re;
                    % 2nd order Central Difference Method
                    A(n,n-2) = a3/h^2;
                    A(n,n-1) = a2 - 4*a3/h^2;
                    A(n,n)   = a1*h^2 - 2*a2 + 6*a3/h^2;
                    A(n,n+1) = a2 - 4*a3/h^2;
                    A(n,n+2) = a3/h^2;
                end

                n = NMAX-1;
                a1 = -U(n)*alpha^3 - Upp(n)*alpha + 1i*alpha^4/Re;
                a2 = U(n)*alpha - 2*1i*alpha^2/Re;
                a3 = 1i/Re;
                % 2nd order Central Difference Method
                A(n,n-2) = a3/h^2;
                A(n,n-1) = a2 - 4*a3/h^2;
                A(n,n)   = a1*h^2 - 2*a2 + 6*a3/h^2;
                A(n,n+1) = a2 - 4*a3/h^2;

                A = A(2:end-1,2:end-1); % remove first/last rows/cols of A (0 BCs)

                % generate B
                B = zeros(NMAX);    % initialize B
                for n = 2:NMAX-1
                    % coefficients as defined in my paper
                    b1 = -alpha^2;
                    b2 = 1;
                    % 2nd order Central Difference Method
                    B(n,n-1) = b2;
                    B(n,n)   = b1*h^2 - 2*b2;
                    B(n,n+1) = b2;
                end
                B = B(2:end-1,2:end-1); % remove first/last rows/cols of B (0 BCs)


                [V,e] = eig(B\A);   % invert B, compute eigenvalues of LHS
                e = diag(e);    % turn diagonal eigenvalue matrix into vector
                
                % store most unstable eigenvalue
                [maxe,oindex] = max(imag(e));
                fprintf('c = %f + %fi\n',real(e(oindex)),imag(e(oindex)));
                s_o(l,m1,m2) = e(oindex);

                
                % Plot eigenvalues of the discrete system
                figure(FIG_SPA_EIG)
                plot(real(e),imag(e),'.','MarkerSize',6)
                axis tight
                xlabel('\omega_r')
                ylabel('\omega_i')
                title(sprintf('alpha = %.2f+%.2fi, Re = %d',real(alpha),imag(alpha),Re));
            end
            
            % store eigenvalues
            for n = 1:length(s_o(l,m1,:))
                omi(n) = imag(s_o(l,m1,n));
                omr(n) = real(s_o(l,m1,n));
            end
            
            zerocount = 0;
            % need to find where omega_i = 0
            for n = 1:length(s_o(l,m1,:))-1
                % if there is a change in sign
                if ( sign(omi(n+1))~=sign(omi(n)) )
                    % estimate the zero by linear interpolation
                    omega = interp1(omi(n:n+1),omr(n:n+1),0);
                    yy = spline(omr,omi); %define a cubic spline
                    myfun = @(x) ppval(yy,x);
                    omega = fzero(myfun,omega); % compute more accurate zero
                    if ( omega > 0 ) % if omega_i > 0, store it for plot
                        zerocount = zerocount+1;
                        num_pts(m1) = num_pts(m1) + 1;
                        re_plot(m1,num_pts(m1)) = Re;
                        om_plot(m1,num_pts(m1)) = omega/Re;
                    end
                end
            end
            fprintf('Omega crosses the real axis %d times\n',zerocount);
            
            % plot alpha-omega map
            figure(FIG_SPA_MAP)
            subplot(1,2,1)
            grid on
            plot(real(s_alphas(m1,:)),imag(s_alphas(m1,:)));
            xlabel('\alpha_r');
            ylabel('\alpha_i');
            subplot(1,2,2)
            grid on
            plot(omr,omi);
            xlabel('\omega_r');
            ylabel('\omega_i');
            if ( oldzerocount > 0 && zerocount == 0 )
                break
            else
                oldzerocount = zerocount;
            end
        end
    end
    elapsed_time = toc; % compute elapsed time
    fprintf('Done (%f sec)\n',elapsed_time);
end

%% OUTPUT RESULTS
if ( run_temporal_analysis )
    cilevels = [0,.005,.01,.015,.0019];
    figure(FIG_TEM_CONT);
    contour(t_Res,t_alphas,imag(t_c),cilevels);
    set(gca,'XScale','log');
    xlabel('$Re_\delta$','Interpreter','latex','Fontsize',14);
    ylabel('$\bar{\alpha}$','Interpreter', 'latex','Fontsize',14);
    colorbar
    grid on;
end

if ( run_spatial_analysis )
    figure(FIG_SPA_CONT)
    prop = ['b.';'r.';'g.';'k.';'c.';'m.'];
    for n = 1:length(s_alphais)
        if ( num_pts(n) > 0 )
            hold on
            plot(re_plot(n,1:num_pts(n)),om_plot(n,1:num_pts(n)),prop(n,:),'MarkerSize',6);
            hold off
        end
    end
    legend('0','-1e-6','-2e-6','-4e-6','-6e-6','-7e-6','location','best')
    
    set(gca,'Xscale','log');
    set(gca,'Yscale','log');
    title('Spatial Stability Curves');
    xlabel('Re_\delta');
    ylabel('\omega / Re_\delta');
end

% uprop = ['k- ';'k: ';'k--'];
uprop = ['b';'r';'g'];
figure(FIG_VEL_PROFILE);
plot(U,xi,uprop(1,:),Up,xi,uprop(2,:),Upp,xi,uprop(3,:))
title('Blasius Velocity Profile')
xlabel('U, U'', U''''')
ylabel('\xi = y/\delta')%,'Interpreter','latex','FontSize',16)
legend('U','U''','U''''','location','best')

%CS

function [ t, w ] = cs_rk4( fun, t0, tf, y0, h )
% cs_rk4 Runge-Kutta 4th Order Numerical Integration
%   [ t, y ] = cs_rk4( fun, t0, tf, y0, h )
%
%   input   description
%   -----   -----------
%   fun     function handle for ydot function
%   t0      initial time
%   tf      final time
%   y0      starting value(s) (column vector)
%   h       time step (delta time)
%
%   output  description
%   ------  -----------
%   t       row vector containing solution times
%   y       row vector(s) containing solution values
%
% Example:
%   [ t, y ] = cs_rk4( @fun, 1, 3, 0, .2 )
%
%   function yp = fun( t, y )
%       yp = 1 + y/t + (y/t)^2;;
%   end
%
%
% Programmed by:
%   Christopher Simpson
%   27 January 2015
%   www.cdsimpson.net

t = [t0:h:tf];
w = zeros(length(y0),length(t)); % initialize w array
w(:,1) = y0;
% fprintf('\nstp\teqn\tti\t\t\twi\t\t\tk1\t\t\tk2\t\t\tk3\t\t\tk4\t\t\twi+1\n');
% fprintf('----------------------------------------------------------------------------------------\n');
for i = 1:(length(t)-1);
    ti = t(i);
    wi = w(:,i);
    k1 = h*fun(ti,wi);
    k2 = h*fun(ti+h/2,wi+k1/2);
    k3 = h*fun(ti+h/2,wi+k2/2);
    k4 = h*fun(ti+h,wi+k3);
    
    w(:,i+1) = wi + 1/6*(k1+2*k2+2*k3+k4);
%     for j = 1:length(w(:,i))
%         fprintf('%d\t%d\t%.6f\t%.6f\t%.6f\t%.6f\t%.6f\t%.6f\t%.6f\n',i,j,ti,wi(j),k1(j),k2(j),k3(j),k4(j),w(j,i+1));
%     end
%     fprintf('----------------------------------------------------------------------------------------\n');
end

end

Substituting ode15s works, so it’s best to not change that.

.

[Mohamed Hajjaj]

I have no problem about the excution and understanding the code itself, but in the first figure above, I need to extract the numerical values of certaion curve (level ci =0 for example) for t_Res, t_alphas that corespond to this curve?

cilevels = [0,.005,.01,.015,.0019];

figure(FIG_TEM_CONT);

contour(t_Res,t_alphas,imag(t_c),cilevels);

[Star Strider]

I suggest adapting the code I already wrote to the problem of extracting the contours. The author is using R2014b or R2015a at the latest, considering that the code is dated March 2015.

If you want to extract more than one contour, the code changes to this:

[x,y,c] = peaks(20);
cilevels = [0.0,.005,0.019,0.25];

figure
[c,h] = contour(x,y,c,cilevels);

axlim = axis;

Lvls = h.LevelList

Lvls = 1×4

0 0.0050 0.0190 0.2500

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

for k1 = 1:numel(Lvls)
    Level = Lvls(k1)                                        % Choose Level
    startIdx = find(c(1,:) == Level);
    vlen = c(2,startIdx);
    for k2 = 1:numel(startIdx)
        % Q = [startIdx(k); vlen(k); startIdx(k)+vlen(k)]
        idxrng = startIdx(k2)+1 : startIdx(k2)+vlen(k2);
        xc{k1,k2} = c(1,idxrng);
        yc{k1,k2} = c(2,idxrng);
    end
end

Level = 0

Level = 0.0050

Level = 0.0190

Level = 0.2500

% xc

figure
hold on
for k1 = 1:size(xc,1)
    for k2 = 1:size(xc,2)
        hp{k1,k2} = plot(xc{k1,k2}, yc{k1,k2}, DisplayName=["("+k1+","+k2+")"]);
    end
end
hold off
axis(axlim)
legend([hp{:}], Location='best')

Your contours will probably not be disjointed as these are, although they could be.

You may need to write the author to find out what the problem is with using ode15s rather than the provided integrator, with the code then throwing the indexing error. The code is not internally documented, wo it’s difficult for me to figure out how it works.

.

EDIT — (4 Mar 2025 at 12:25)

NOTE: My code pulls out every part of every contour (they are discontinuous) and then plots it. However, I suspect from looking at the code that creates this, that the contours may not be contours at all and instead are different trajectories of a specific variable for specific values of a particular parameter. I cannot get the code to work (why they didn’t use the MATLAB ODE solvers in code written in 2015 mystifies me), so I cannot verify that.

.

[Star Strider]

‘Ok, If I need to save x and y data into a data file to print or to read it, what is the option here?’

In my code:

save('YourFileName.mat', 'xc', 'yc')

Since these are cell arrays (required because the sizes vary), no other format will work, even using writecell (and readcell).

My complete code with that addition —

[x,y,c] = peaks(20);
cilevels = [0.0,.005,0.019,0.25];

figure
[c,h] = contour(x,y,c,cilevels);

axlim = axis;

Lvls = h.LevelList

Lvls = 1×4

0 0.0050 0.0190 0.2500

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

for k1 = 1:numel(Lvls)
    Level = Lvls(k1)                                        % Choose Level
    startIdx = find(c(1,:) == Level);
    vlen = c(2,startIdx);
    for k2 = 1:numel(startIdx)
        % Q = [startIdx(k); vlen(k); startIdx(k)+vlen(k)]
        idxrng = startIdx(k2)+1 : startIdx(k2)+vlen(k2);
        xc{k1,k2} = c(1,idxrng);
        yc{k1,k2} = c(2,idxrng);
    end
end

Level = 0

Level = 0.0050

Level = 0.0190

Level = 0.2500

xc

xc = 4x3 cell array
    {1x30 double}    {1x27 double}    {0x0  double}
    {1x29 double}    {1x9  double}    {1x36 double}
    {1x13 double}    {1x59 double}    {0x0  double}
    {1x66 double}    {1x9  double}    {0x0  double}

yc 

yc = 4x3 cell array
    {1x30 double}    {1x27 double}    {0x0  double}
    {1x29 double}    {1x9  double}    {1x36 double}
    {1x13 double}    {1x59 double}    {0x0  double}
    {1x66 double}    {1x9  double}    {0x0  double}

save('TestWrite.mat', 'xc','yc')                    % 'save' Arrays To '.mat' File

xcyc = load('TestWrite.mat')                        % 'load' Arrays From '.mat' File

xcyc = struct with fields:
    xc: {4x3 cell}
    yc: {4x3 cell}

xc = xcyc.xc

xc = 4x3 cell array
    {1x30 double}    {1x27 double}    {0x0  double}
    {1x29 double}    {1x9  double}    {1x36 double}
    {1x13 double}    {1x59 double}    {0x0  double}
    {1x66 double}    {1x9  double}    {0x0  double}

yc = xcyc.yc

yc = 4x3 cell array
    {1x30 double}    {1x27 double}    {0x0  double}
    {1x29 double}    {1x9  double}    {1x36 double}
    {1x13 double}    {1x59 double}    {0x0  double}
    {1x66 double}    {1x9  double}    {0x0  double}

figure
hold on
for k1 = 1:size(xc,1)
    for k2 = 1:size(xc,2)
        hp{k1,k2} = plot(xc{k1,k2}, yc{k1,k2}, DisplayName=["("+k1+","+k2+")"]);
    end
end
hold off
% axis(axlim)
legend([hp{:}], Location='best')

The save function will not work with cell arrays and the '-ascii' tag.

Your options are quite limited with these data.

.

[Star Strider]

I ran the code on MATLAB Online, saving the results in a .mat file. (It took a bit more than 30 minutes, so it is not possible to run it here, even using the faster and more efficient MATLAB ODE integrators. It iterates too many times.)

This code loads the Orr_Somerfield_Results.mat file, uses the original code to plot the contours, and then uses my code to extract the contours, printing the (x,y) coordinates of the requested contours, taken from the code and the matrix of results. You can copy my code here, download the Orr_Somerfield_Results.mat file from this Comment, and run it locally to get the results you want.

Try this —

load('Orr_Somerfield_Results.mat')

    cilevels = [0,.005,.01,.015,.0019];
    figure
    [c,h] = contour(t_Res,t_alphas,imag(t_c),cilevels, ShowText=1);
    set(gca,'XScale','log');
    xlabel('$Re_\delta$','Interpreter','latex','Fontsize',14);
    ylabel('$\bar{\alpha}$','Interpreter', 'latex','Fontsize',14);
    colormap(turbo)
    colorbar
    grid on;

figure
surf(t_Res,t_alphas,imag(t_c), FaceAlpha=0.5, EdgeAlpha=0.5)
hold on
contour3(t_Res,t_alphas,imag(t_c),cilevels, 'g', LineWidth=2, ShowText=1);
hold off
set(gca,'XScale','log');
xlabel('$Re_\delta$','Interpreter','latex','Fontsize',14);
ylabel('$\bar{\alpha}$','Interpreter', 'latex','Fontsize',14);
colormap(turbo)
colorbar
grid on;
zlim([-0.01 0.025])
title('Surface Plot With Contours')

    
    
Lvls = h.LevelList

Lvls = 1×5

0 0.0019 0.0050 0.0100 0.0150

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

for k1 = 1:numel(Lvls)
    Level = Lvls(k1)                                        % Choose Level
    startIdx = find(c(1,:) == Level);
    vlen = c(2,startIdx);
    for k2 = 1:numel(startIdx)
        % Q = [startIdx(k); vlen(k); startIdx(k)+vlen(k)]
        idxrng = startIdx(k2)+1 : startIdx(k2)+vlen(k2);
        xc{k1,k2} = c(1,idxrng);
        yc{k1,k2} = c(2,idxrng);
    end
end

Level = 0

Level = 0.0019

Level = 0.0050

Level = 0.0100

Level = 0.0150

xc

xc = 5x1 cell array
    {[100000 9.7305e+04 8.8862e+04 8.4810e+04 7.8965e+04 7.0170e+04 6.9952e+04 6.2355e+04 6.0509e+04 5.5410e+04 5.2992e+04 4.9239e+04 4.5477e+04 4.3755e+04 ... ] (1x127 double)}
    {[100000 8.8862e+04 8.7972e+04 7.8965e+04 7.5216e+04 7.0170e+04 6.4421e+04 6.2355e+04 5.5410e+04 5.5285e+04 4.9239e+04 4.7671e+04 4.3755e+04 4.1180e+04 ... ] (1x121 double)}
    {[100000 8.8862e+04 8.8777e+04 7.8965e+04 7.5976e+04 7.0170e+04 6.5187e+04 6.2355e+04 5.6054e+04 5.5410e+04 4.9239e+04 4.8347e+04 4.3755e+04 4.1797e+04 ... ] (1x116 double)}
    {[100000 9.3919e+04 8.8862e+04 7.9905e+04 7.8965e+04 7.0170e+04 6.8381e+04 6.2355e+04 5.8758e+04 5.5410e+04 5.0631e+04 4.9239e+04 4.3755e+04 4.3723e+04 ... ] (1x103 double)}
    {[ 100000 9.4434e+04 8.8862e+04 7.9261e+04 7.8965e+04 7.0170e+04 6.7192e+04 6.2355e+04 5.7387e+04 5.5410e+04 4.9264e+04 4.9239e+04 4.3755e+04 4.2404e+04 ... ] (1x87 double)}

yc 

yc = 5x1 cell array
    {[0.3950 0.4007 0.4257 0.4340 0.4534 0.4666 0.4673 0.4938 0.5006 0.5248 0.5339 0.5538 0.5672 0.5773 0.5969 0.6005 0.6254 0.6338 0.6565 0.6672 0.6868 0.7005 ... ] (1x127 double)}
    {[0.3746 0.3984 0.4007 0.4231 0.4340 0.4483 0.4673 0.4745 0.5001 0.5006 0.5262 0.5339 0.5528 0.5672 0.5798 0.6005 0.6077 0.6338 0.6357 0.6632 0.6672 0.6913 ... ] (1x121 double)}
    {[0.3425 0.3672 0.3674 0.3919 0.4007 0.4173 0.4340 0.4432 0.4673 0.4697 0.4962 0.5006 0.5229 0.5339 0.5500 0.5672 0.5775 0.6005 0.6054 0.6336 0.6338 0.6614 ... ] (1x116 double)}
    {[0.2877 0.3008 0.3116 0.3341 0.3364 0.3615 0.3674 0.3870 0.4007 0.4132 0.4340 0.4399 0.4671 0.4673 0.4939 0.5006 0.5211 0.5339 0.5485 0.5672 0.5762 0.6005 ... ] (1x103 double)}
    {[ 0.2233 0.2342 0.2450 0.2675 0.2681 0.2914 0.3008 0.3158 0.3341 0.3412 0.3674 0.3675 0.3933 0.4007 0.4196 0.4340 0.4464 0.4673 0.4735 0.5006 0.5008 0.5270 ... ] (1x87 double)}

save('TestWrite.mat', 'xc','yc')                    % 'save' Arrays To '.mat' File

xcyc = load('TestWrite.mat')                        % 'load' Arrays From '.mat' File

xcyc = struct with fields:
    xc: {5x1 cell}
    yc: {5x1 cell}

xc = xcyc.xc

xc = 5x1 cell array
    {[100000 9.7305e+04 8.8862e+04 8.4810e+04 7.8965e+04 7.0170e+04 6.9952e+04 6.2355e+04 6.0509e+04 5.5410e+04 5.2992e+04 4.9239e+04 4.5477e+04 4.3755e+04 ... ] (1x127 double)}
    {[100000 8.8862e+04 8.7972e+04 7.8965e+04 7.5216e+04 7.0170e+04 6.4421e+04 6.2355e+04 5.5410e+04 5.5285e+04 4.9239e+04 4.7671e+04 4.3755e+04 4.1180e+04 ... ] (1x121 double)}
    {[100000 8.8862e+04 8.8777e+04 7.8965e+04 7.5976e+04 7.0170e+04 6.5187e+04 6.2355e+04 5.6054e+04 5.5410e+04 4.9239e+04 4.8347e+04 4.3755e+04 4.1797e+04 ... ] (1x116 double)}
    {[100000 9.3919e+04 8.8862e+04 7.9905e+04 7.8965e+04 7.0170e+04 6.8381e+04 6.2355e+04 5.8758e+04 5.5410e+04 5.0631e+04 4.9239e+04 4.3755e+04 4.3723e+04 ... ] (1x103 double)}
    {[ 100000 9.4434e+04 8.8862e+04 7.9261e+04 7.8965e+04 7.0170e+04 6.7192e+04 6.2355e+04 5.7387e+04 5.5410e+04 4.9264e+04 4.9239e+04 4.3755e+04 4.2404e+04 ... ] (1x87 double)}

yc = xcyc.yc

yc = 5x1 cell array
    {[0.3950 0.4007 0.4257 0.4340 0.4534 0.4666 0.4673 0.4938 0.5006 0.5248 0.5339 0.5538 0.5672 0.5773 0.5969 0.6005 0.6254 0.6338 0.6565 0.6672 0.6868 0.7005 ... ] (1x127 double)}
    {[0.3746 0.3984 0.4007 0.4231 0.4340 0.4483 0.4673 0.4745 0.5001 0.5006 0.5262 0.5339 0.5528 0.5672 0.5798 0.6005 0.6077 0.6338 0.6357 0.6632 0.6672 0.6913 ... ] (1x121 double)}
    {[0.3425 0.3672 0.3674 0.3919 0.4007 0.4173 0.4340 0.4432 0.4673 0.4697 0.4962 0.5006 0.5229 0.5339 0.5500 0.5672 0.5775 0.6005 0.6054 0.6336 0.6338 0.6614 ... ] (1x116 double)}
    {[0.2877 0.3008 0.3116 0.3341 0.3364 0.3615 0.3674 0.3870 0.4007 0.4132 0.4340 0.4399 0.4671 0.4673 0.4939 0.5006 0.5211 0.5339 0.5485 0.5672 0.5762 0.6005 ... ] (1x103 double)}
    {[ 0.2233 0.2342 0.2450 0.2675 0.2681 0.2914 0.3008 0.3158 0.3341 0.3412 0.3674 0.3675 0.3933 0.4007 0.4196 0.4340 0.4464 0.4673 0.4735 0.5006 0.5008 0.5270 ... ] (1x87 double)}

figure
hold on
for k1 = 1:size(xc,1)
    for k2 = 1:size(xc,2)
        hp{k1,k2} = plot(xc{k1,k2}, yc{k1,k2}, DisplayName=string(Lvls(k1)));
    end
end
hold off
grid on                                    % <— OPTIONAL
Ax = gca;
Ax.XScale='log';
% axis(axlim)
% axis('square')
lgd = legend([hp{:}], Location='best');
title(lgd, 'Contour Levels')
xlabel('$\Re_\delta$','Interpreter','latex','Fontsize',14);
ylabel('$\bar{\alpha}$','Interpreter', 'latex','Fontsize',14);

.

[Mohamed Hajjaj]

It is geart and I'm grateful to you for all your help and your effort.

[Star Strider]

My pleasure!

If my Answer helped you solve your problem, please Accept it!

.

Answer 2

[Cris LaPierre]

I think you want the (x,y) coordinates of your levels. That can be obtained using this syntax: M = contour(___)

However, be sure to read the documentation to understand how the 2xN output array is formatted:

Contour matrix, returned as a two-row matrix of following form.

Z1, x1,1, x1,2, ..., x1,N1, Z2, x2,1, x2,2, ..., x2,N2, Z3, ...

N1, y1,1, y1,2, ..., y1,N1, N2, y2,1, y2,2, ..., y2,N2, N3, ...

The columns of the matrix define the contour lines. Each contour line starts with a column containing Z and N values:

• Zi — The height of the ith contour line
• Ni — The number of vertices in the ith contour line
• (xij, yij) — The coordinates of the vertices for the ith contour line, where j ranges from 1 to Ni

Conside this example:

p = peaks(25);
m = contour(p)

m = 2×259

-6.0000 14.4744 14.0000 13.0767 13.4954 14.0000 14.5773 14.4744 -4.0000 15.3227 15.0000 14.0000 13.0000 12.1330 12.0000 11.5885 11.9601 12.0000 13.0000 13.3965 14.0000 15.0000 15.3397 16.0000 16.1878 16.1552 16.0000 15.3227 -2.0000 16.0658 7.0000 6.0000 5.7662 6.0000 7.0000 7.1890 7.0000 6.0000 19.0000 5.0000 4.7869 4.4921 4.5628 5.0000 5.1433 6.0000 7.0000 7.0338 7.8094 8.0000 8.2328 8.1836 8.0000 7.2888 7.0000 6.0000 5.7510 5.0000 27.0000 4.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Inspecting the output m, I see that the first 2 entries correspopnd to the contour levels for -6 and -4. There are 7 (x,y) coordinates designating the -6 contour, and 15 for -4. You won't know the number of points a priori, so you'll probably need a loop to extract the data. Here's one way.

i = 1;
data = [];
while i < length(m)
    tmp.clvl = m(1,i);
    num = m(2,i);

    tmp.x = m(1,i+1:i+num);
    tmp.y = m(2,i+1:i+num);
    i = i+num+1;
    
    data = [data,tmp];
end
data(1)

ans = struct with fields:
    clvl: -6
       x: [14.4744 14 13.0767 13.4954 14 14.5773 14.4744]
       y: [6 5.7662 6 7 7.1890 7 6]

Now that I have that data, I can plot it for comparison.

figure
data(1).clvl

ans = -6

plot(data(1).x,data(1).y)

Note that some contour levels have multiple independent contours (0, for example here). Each will appear as their own item in the output array. Said another way, there are 2 entries for 0 here; entries 5 and 6.

data(5)

ans = struct with fields:
    clvl: 0
       x: [7.4524 7.1551 7.0405 7.0568 7.2600 8 8.0478 9 9.1791 10 10.8642 11 12 12.5792 13 13.9966 14 15 15.6699 16 17 17.3197 18 18.8752 19 20 20.6058 21 ... ] (1x35 double)
       y: [1 2 3 4 5 5.8471 6 6.7570 7 7.4351 8 8.0587 8.6493 9 9.3414 10 10.0068 10.5148 10 9.8749 9.1884 9 8.4892 8 7.8475 7.1987 7 6.5074 6 5.9504 5.1578 5 4.2651 4 3.1712]

data(6)

ans = struct with fields:
    clvl: 0
       x: [1 2 3 3.0329 4 5 5.0995 6 7 7.0792 8 8.7235 9 10 10.0422 10.6609 10.4867 10 9.9630 9.3154 9 8.4278 8 7.3619 7 6.1244 6 5 4.2841 4 3 2 1.7134 1]
       y: [20.9823 20.6637 20.0768 20 19.6921 19.1473 19 18.6728 18.0895 18 17.5227 17 16.8355 16.0295 16 15 14 13.0895 13 12 11.5847 11 10.4536 10 9.5173 9 ... ] (1x34 double)

So ensure you set a cilevel for the line whose (x,y) data you want.

[Mohamed Hajjaj]

Thank you very much for your response, It seems that is near to my requisite . I have tried your example, it was some error:

clvl: 0

x: [1x35 double]

y: [1x35 double]

for data(5),

clvl: 0

x: [1x34 double]

y: [1x34 double]

for data(6)

[Cris LaPierre]

I don't see an error. The output you pasted is showing the size of the vectors stored in the X and Y fields rather than the actual values, but that's not an error. There are 35 values for data(5) and 34 for data(6)

[Mohamed Hajjaj]

Ok, If I need to save x and y data into a data file to print or to read it, what is the option here?

[Cris LaPierre]

The easiest way to view the values is display them on the screen like I did (red box below). Another option is to view the variable in the Variable Editor. Double-click on the variable in the Workspace (red arrow) to open it in the Variable Editor (orange box below)

If a field is too big to show, its size is displayed instead. To see the contents, double click on the hyperlink size in a cell (green arrow above) to open that in the Variable Editor (image below).

If you do want to save to file, there are many options. I would combine them into a table and use writetable to save them to a file.

% select data to save
clvl = data(5).clvl
X = data(5).x;
Y = data(5).y;

% save to a txt file
tbl = table(X,Y);
writetable(tbl,"contour_" + clvl + ".txt")

[Mohamed Hajjaj]

It helps me a lot, thank you very much.