How to plot a graph correctly?

I have a function of two variables F(t,d) and am building a plane, and I want to plot a function where F(t,d) - C_2=0, I use function 'contour', but I get a few curves and if I take contour(t,d,F) or contour(d,t,F) only the location of the coordinate axes changes, but the graph remains the same. In conclusion, you need to get a graph t(d) for F() - C_2=0.

My code:

%% initial conditions
global k0 h_bar ksi m E C_2
Ef = 2.77*10^3; 
Kb = physconst('boltzmann'); % 1.38*10^(-23)
m = 9.1093837*10^(-31);
Tc = 1.2;
ksi = 10^(-9);
h_bar = (1.0545726*10^(-34));
k0 = (ksi/h_bar)*sqrt(2.*m.*pi.*Kb.*Tc);
C_2 = 6.81;
t = linspace(0.1, 2);
d = linspace(10, 1000);
%[TT,DD] = meshgrid(t,d);
%contour(TT, DD, @f_calc);
for i=1:numel(d)
    for j = 1:numel(t)
        F(i,j) = f_calc(t(j),d(i));
    end
end
%plot(t,F)
contour(d,t,F)
%xlabel('d') 
%ylabel('t')
%surf(t,d,F)
function F = f_calc(t, d) 
global k0 h_bar ksi m C_2
nD = floor(375/(2*pi.*t*1.2) - 0.5);
F = 0;
for k = 0:nD
    F = F + 1/(2*k+1).*imag(f_lg(k,t,d)+1i.*d.*k0.*((f_p1(k,t)-f_p2(k,t))./2)+1i*f_arg_1(k,t,d)-1i*f_arg_2(k,t,d));
end
F = -(1/d).*F - 7.826922968141167;
end
function p1 = f_p1(n,t)
p1 = ((1+1i)./sqrt(2)).*sqrt(t.*(2.*n+1));
end
function p2 = f_p2(n,t)
p2 = sqrt(3601+1i.*t.*(2.*n+1));
end 
function arg_1 = f_arg_1(n,t,d)
global k0
arg_1 = angle(1+exp(-1i.*d*k0.*f_p1(n,t)));
end
function arg_2 = f_arg_2(n,t,d)
global k0
arg_2 = angle(1+exp(-1i.*d*k0.*f_p2(n,t)));
end
function n_lg = f_lg(n,t,d)
global k0;
arg_of_lg = (1+exp(-1i.*d.*k0.*f_p1(n,t)))/(1+exp(-1i.*d.*k0.*f_p2(n,t)));
n_lg = log(abs(arg_of_lg));
end

Answer 1

[Star Strider]

I am not certain that I understand what you want to do.

If you only want the contour at ‘F=0’, supply that as the fourth, ‘Levels’ argument to contour:

contour(d,t,F,[0 0])

Try this —

%% initial conditions
global k0 h_bar ksi m E C_2
Ef = 2.77*10^3; 
Kb = physconst('boltzmann'); % 1.38*10^(-23)
m = 9.1093837*10^(-31);
Tc = 1.2;
ksi = 10^(-9);
h_bar = (1.0545726*10^(-34));
k0 = (ksi/h_bar)*sqrt(2.*m.*pi.*Kb.*Tc);
C_2 = 6.81;
t = linspace(0.1, 2);
d = linspace(10, 1000);
%[TT,DD] = meshgrid(t,d);
%contour(TT, DD, @f_calc);
for i=1:numel(d)
  for j = 1:numel(t)
    F(i,j) = f_calc(t(j),d(i));
  end
end
%plot(t,F)
contour(d,t,F,[0 0])
xlabel('d') 
ylabel('t')

%surf(t,d,F)
function F = f_calc(t, d) 
global k0 h_bar ksi m C_2
nD = floor(375/(2*pi.*t*1.2) - 0.5);
F = 0;
for k = 0:nD
    F = F + 1/(2*k+1).*imag(f_lg(k,t,d)+1i.*d.*k0.*((f_p1(k,t)-f_p2(k,t))./2)+1i*f_arg_1(k,t,d)-1i*f_arg_2(k,t,d));
end
F = -(1/d).*F - 7.826922968141167;
end
function p1 = f_p1(n,t)
p1 = ((1+1i)./sqrt(2)).*sqrt(t.*(2.*n+1));
end
function p2 = f_p2(n,t)
p2 = sqrt(3601+1i.*t.*(2.*n+1));
end 
function arg_1 = f_arg_1(n,t,d)
global k0
arg_1 = angle(1+exp(-1i.*d*k0.*f_p1(n,t)));
end
function arg_2 = f_arg_2(n,t,d)
global k0
arg_2 = angle(1+exp(-1i.*d*k0.*f_p2(n,t)));
end
function n_lg = f_lg(n,t,d)
    global k0;
    arg_of_lg = (1+exp(-1i.*d.*k0.*f_p1(n,t)))/(1+exp(-1i.*d.*k0.*f_p2(n,t)));
    n_lg = log(abs(arg_of_lg));
end

.

[Dmitry]

This graph should be reversed. Why, if you plot t(d) or d(t), do they not change?

[Star Strider]

Reversing it is easy enough —

%% initial conditions
global k0 h_bar ksi m E C_2
Ef = 2.77*10^3; 
Kb = physconst('boltzmann'); % 1.38*10^(-23)
m = 9.1093837*10^(-31);
Tc = 1.2;
ksi = 10^(-9);
h_bar = (1.0545726*10^(-34));
k0 = (ksi/h_bar)*sqrt(2.*m.*pi.*Kb.*Tc);
C_2 = 6.81;
t = linspace(0.1, 2);
d = linspace(10, 1000);
%[TT,DD] = meshgrid(t,d);
%contour(TT, DD, @f_calc);
for i=1:numel(d)
  for j = 1:numel(t)
    F(i,j) = f_calc(t(j),d(i));
  end
end
%plot(t,F)
figure
[c,h] = contour(d,t,F,[0 0]);
xlabel('d') 
ylabel('t')

tv = c(2, 2:end)

tv = 1×191

0.1000 0.1086 0.1162 0.1155 0.1192 0.1227 0.1238 0.1249 0.1251 0.1262 0.1274 0.1284 0.1295 0.1307 0.1317 0.1327 0.1337 0.1349 0.1357 0.1366 0.1375 0.1384 0.1407 0.1434 0.1460 0.1485 0.1517 0.1543 0.1575 0.1576

dv = c(1, 2:end)

dv = 1×191

596.2918 590.0000 580.0000 570.0000 568.4323 560.0000 550.0000 540.0000 530.0000 520.0000 510.0000 500.0000 490.0000 480.0000 470.0000 460.0000 450.0000 440.0000 430.0000 420.0000 410.0000 408.2938 410.0000 420.0000 430.0000 440.0000 450.0000 460.0000 470.0000 470.1067

figure
plot(tv, dv)
xlabel('t')
ylabel('d')
grid

%surf(t,d,F)
function F = f_calc(t, d) 
global k0 h_bar ksi m C_2
nD = floor(375/(2*pi.*t*1.2) - 0.5);
F = 0;
for k = 0:nD
    F = F + 1/(2*k+1).*imag(f_lg(k,t,d)+1i.*d.*k0.*((f_p1(k,t)-f_p2(k,t))./2)+1i*f_arg_1(k,t,d)-1i*f_arg_2(k,t,d));
end
F = -(1/d).*F - 7.826922968141167;
end
function p1 = f_p1(n,t)
p1 = ((1+1i)./sqrt(2)).*sqrt(t.*(2.*n+1));
end
function p2 = f_p2(n,t)
p2 = sqrt(3601+1i.*t.*(2.*n+1));
end 
function arg_1 = f_arg_1(n,t,d)
global k0
arg_1 = angle(1+exp(-1i.*d*k0.*f_p1(n,t)));
end
function arg_2 = f_arg_2(n,t,d)
global k0
arg_2 = angle(1+exp(-1i.*d*k0.*f_p2(n,t)));
end
function n_lg = f_lg(n,t,d)
    global k0;
    arg_of_lg = (1+exp(-1i.*d.*k0.*f_p1(n,t)))/(1+exp(-1i.*d.*k0.*f_p2(n,t)));
    n_lg = log(abs(arg_of_lg));
end

With only one continuous contour, getting the data from the contour function ‘c’ output is straightforward. The first column are the level (here 0), and the number of columns of data (here 191). The first row (2:end) are the ‘x’ data and the second row (2:end) are the ‘y’ data. Reversing the plot simply means reversing those vectors in the plot argument.

What else would you like to do with the ‘c’ (contour) data, ‘tv’ and ‘dv’?

.

EDIT — (26 Jan 2023 at 20:48)

‘I changed the constant that we subtract from F and now there are solutions there.’

Does that affect the contour plot in any way, or is it doing what you want it to do as currently written?

.

[Dmitry]

tv and dv are these vectors of values 't' and 'd' by which the graph is plotted for F = 0?

[Star Strider]

Yes.

[Dmitry]

Do I understand correctly that this is a graph t(d) for F=0?

[Star Strider]

Yes, because:

[c,h] = contour(d,t,F,[0 0]);

plots the contour for ‘F=0’ and the ‘c’ result (in this instance) contains only the contour for that level.

Fortunately, it is a continuous contour, so it can be directly accessed as a single vector for each (x,y) (or here (d,t)) coordinate. (If it were not a single contour, it would still be possible to retrieve it, however that would be more complicated to code.)

[Dmitry]

In conclusion, we have a plane F(t,d) that intersects the plane X Y (I mean d t) by functuon "contour" we make a slice for F=0, look from above and see the graph t(d). am I right?

[Torsten]

We have a surface z = F(t,d) over a rectangle [tmin,tmax] x [dmin,dmax], cut this surface by the plane z=0, look from above and see the curve of tuples (t,d) for which F(t,d) = 0.

The same as fimplicit does.

I suspect you didn't choose the search region for d big enough. Contour shows that d must be about 400 - 600 for that the cut with the plane z=0 and the surface is not empty.

[Star Strider]

@Dmitry — As I understand the original problem, yes.

[Dmitry]

Star Strider, Sir, I set the vector 'd' so that the step between the values is 0.1, but if we look at the vector 'c' mf we see that the step for 'd' is much smaller there, how can this be explained?

And why do I have values of t=0 on the final graph for 't' if I count 't{0.1:2}' ?

Code:

%% initial conditions

global k0 h_bar ksi m E C_2

Ef = 2.77*10^3;

Kb = physconst('boltzmann'); % 1.38*10^(-23)

m = 9.1093837*10^(-31);

Tc = 1.2;

ksi = 10^(-9);

h_bar = (1.0545726*10^(-34));

k0 = (ksi/h_bar)*sqrt(2.*m.*pi.*Kb.*Tc);

C_2 = 6.81;

t = linspace(0.1, 2,901);

d = linspace(10, 100,901);

%[TT,DD] = meshgrid(t,d);

%contour(TT, DD, @f_calc);

for i=1:numel(d)

for j = 1:numel(t)

F(i,j) = f_calc(t(j),d(i));

end

end

%plot(t,F)

figure

[c,h] = contour(t,d,F,[0 0]);

xlabel('d')

ylabel('t')

tv = c(2, 2:end)

dv = c(1, 2:end)

figure

plot(tv, dv)

xlabel('d')

ylabel('t')

grid

%surf(t,d,F)

function F = f_calc(t, d)

global k0 h_bar ksi m C_2

nD = floor(375/(2*pi.*t*1.2) - 0.5);

F = 0;

for k = 0:nD

F = F + 1/(2*k+1).*imag(f_lg(k,t,d)+1i.*d.*k0.*((f_p1(k,t)-f_p2(k,t))./2)+1i*f_arg_1(k,t,d)-1i*f_arg_2(k,t,d));

end

F = -(1/d).*F - 7.826922968141167;

end

function p1 = f_p1(n,t)

p1 = ((1+1i)./sqrt(2)).*sqrt(t.*(2.*n+1));

end

function p2 = f_p2(n,t)

p2 = sqrt(3601+1i.*t.*(2.*n+1));

end

function arg_1 = f_arg_1(n,t,d)

global k0

arg_1 = angle(1+exp(-1i.*d*k0.*f_p1(n,t)));

end

function arg_2 = f_arg_2(n,t,d)

global k0

arg_2 = angle(1+exp(-1i.*d*k0.*f_p2(n,t)));

end

function n_lg = f_lg(n,t,d)

global k0;

arg_of_lg = (1+exp(-1i.*d.*k0.*f_p1(n,t)))/(1+exp(-1i.*d.*k0.*f_p2(n,t)));

n_lg = log(abs(arg_of_lg));

end

[Dmitry]

Star Strider, if I want to delete from the vector 'c' points at which 't'=0, can I use such an insertion?:

for i=1:numel(c)

for j = 1:numel(c)

if c[1;2:end] = 0;

delete(c[i,j])

end

end

after block:

[c,h] = contour(t,d,F,[0 0]);

xlabel('d')

ylabel('t')

[Star Strider]

The new contour result turned out to be something of a challenge. The problem is that unlike the first version, this version produces several (13) different contours, and the problem was stitching them together correctly. Also, it was necessary to sort them by the ‘yv’ vector rather than the ‘xv’ vector (perhaps the sorting is not strictly necessary, however I wanted to be certain because weird lines in a plot usually indicate unsorted data).

This takes too long to run here, so I can’t demonstrate it, however it does run successfully —

%% initial conditions
% global k0 h_bar ksi m E C_2
Ef = 2.77*10^3; 
Kb = physconst('boltzmann'); % 1.38*10^(-23)
m = 9.1093837*10^(-31);
Tc = 1.2;
ksi = 10^(-9);
h_bar = (1.0545726*10^(-34));
k0 = (ksi/h_bar)*sqrt(2.*m.*pi.*Kb.*Tc);
C_2 = 6.81;
t = linspace(0.1, 2,901);
d = linspace(10, 100,901);
%[TT,DD] = meshgrid(t,d);
%contour(TT, DD, @f_calc);
for i=1:numel(d)
  for j = 1:numel(t)
    F(i,j) = f_calc(t(j),d(i), k0, h_bar, ksi, m, C_2);
  end
end
%plot(t,F)
figure
[c,h] = contour(t,d,F,[0 0]);
xlabel('d') 
ylabel('t')
% tv = c(2, 2:end);
% dv = c(1, 2:end);
[tv,dv] = getContour0(c,h);
Q1 = [size(c); size(tv)]
figure
plot(tv, dv)
xlabel('d')
ylabel('t')
grid
%surf(t,d,F)
function F = f_calc(t, d, k0, h_bar, ksi, m, C_2) 
% global k0 h_bar ksi m C_2
nD = floor(375/(2*pi.*t*1.2) - 0.5);
F = 0;
for k = 0:nD
    F = F + 1/(2*k+1).*imag(f_lg(k,t,d,k0)+1i.*d.*k0.*((f_p1(k,t)-f_p2(k,t))./2)+1i*f_arg_1(k,t,d,k0)-1i*f_arg_2(k,t,d,k0));
end
F = -(1/d).*F - 7.826922968141167;
end
function p1 = f_p1(n,t)
p1 = ((1+1i)./sqrt(2)).*sqrt(t.*(2.*n+1));
end
function p2 = f_p2(n,t)
p2 = sqrt(3601+1i.*t.*(2.*n+1));
end 
function arg_1 = f_arg_1(n,t,d,k0)
% global k0
arg_1 = angle(1+exp(-1i.*d*k0.*f_p1(n,t)));
end
function arg_2 = f_arg_2(n,t,d,k0)
% global k0
arg_2 = angle(1+exp(-1i.*d*k0.*f_p2(n,t)));
end
function n_lg = f_lg(n,t,d,k0)
%     global k0;
    arg_of_lg = (1+exp(-1i.*d.*k0.*f_p1(n,t)))/(1+exp(-1i.*d.*k0.*f_p2(n,t)));
    n_lg = log(abs(arg_of_lg));
end

function [xv,yv] = getContour0(M,C)
Levels = C.LevelList
    for k = 1:numel(Levels)
        idx = find(M(1,:) == 0);
        ValidV = rem(M(2,idx),1) == 0;
        StartIdx{k,:} = idx(ValidV);
        VLen{k,:} = M(2,StartIdx{k});
%         Q3 = numel(StartIdx{k})
            for k1 = 1:numel(StartIdx{k})
%                 Q4 = StartIdx{k}(k1)
                idxv = StartIdx{k}(k1)+1 : StartIdx{k}(k1)+VLen{k}(k1);          % Index For Contour 'k1'
                xc{k1} = M(1,idxv);
                yc{k1} = M(2,idxv);
            end
    end
    xy = [cell2mat(xc); cell2mat(yc)].';
    xy = sortrows(xy,2);
    xv = xy(:,1).';
    yv = xy(:,2).';

%     figure
%     plot(xv, yv, '-r')
%     grid
%     title('Test Plot')
end

The ‘getContour0’ function gets the components of the zero-contour, stitches them together, and sorts them. I kept a few of the debugging lines in and commented out.

I also removed the global calls and added the appropriate arguments to the functions that use them. Avoid global variables like the plague they are unless they are absolutely necessary. (I have never found them to be necessary at all, other than in early PC versions of MATLAB in the mid-1990s when passing extra parameters was not possible.) See the documentation section on Passing Extra Parameters for an extended discussion on them.

.

[Dmitry]

Star Strider, Why do we have much more values in the vector 'c' than we set for 'd' 't' by linespace()?

And we get 2 graphs at the output, how do they differ from each other?

[Star Strider]

With respect to the differences in ‘c’, and ‘tv’ or ‘dv’, it is because there are 13 columns in ‘c’ that indicate the level and the number of values in the corresponding vector segments. (In this instance, there are 13 contour segments, not 1 contour as in the earlier example.) The difference in lengths accounts for deleting those columns (that are simply there to provide information about the contours). That accounts for the difference between the lengths of ‘c’, and ‘tv’ or ‘dv’. The vectors returned from ‘getContour0’ each have 18344 values.

The contours themselves are going to be different from the values in the linspace calls because the linspace calls define the matrices, not the contours. I doubt that there is any specific way of relating the matrix sizes to the lenghts of the contour vectors, other than that they are likely proportional only (and that being due to the resolution of the matrix on which contour operates), and do not have any specific functional relationship.

.

[Dmitry]

Star Strider, are these all contours for the zero level?

[Star Strider]

As far as I can tell, yes.

For whatever reason, contour occasionally breaks up a contour at one level into several different contours. Sometimes this is obvious, for example taking the contours at 0 of the peaks function, although here it is less obvious, at least to me, since I do not understand what you are doing. There are apparently non-zero regions between the zero segments here, and that breaks the 0 contour into disjointed segments.

You would likely have to look at a small section of the contour plot near the 0 region and plot all the available contours in that region to see what it is doing. (It may be necessary to specify those contours as described in the documentation section on levels. Use the xlim function to restrict the region of the contour plot, for example [400 600].)

I leave this to you because I do not understand what your code calculates, so I cannot interpret that result.

Answer 2

[Torsten]

Try "fimplicit".

[Dmitry]

When I use this function I have an empty graph

[Torsten]

Sure. At least for the ranges of t and d where we searched, there was no solution pair for the equation F(t,d) - C_2 = 0. Did you find one in the meantime ?

[Dmitry]

If we look at this solution, we can see that there is a line but if you plot contour(d,t,F,[0 0]) orcontour(d,t,F,[0 0]), do they not change and that's weird. I changed the constant that we subtract from F and now there are solutions there.

%% initial conditions
global k0 h_bar ksi m E C_2
Ef = 2.77*10^3; 
Kb = physconst('boltzmann'); % 1.38*10^(-23)
m = 9.1093837*10^(-31);
Tc = 1.2;
ksi = 10^(-9);
h_bar = (1.0545726*10^(-34));
k0 = (ksi/h_bar)*sqrt(2.*m.*pi.*Kb.*Tc);
C_2 = 6.81;
t = linspace(0.1, 2);
d = linspace(10, 1000);
%[TT,DD] = meshgrid(t,d);
%contour(TT, DD, @f_calc);
for i=1:numel(d)
    for j = 1:numel(t)
        F(i,j) = f_calc(t(j),d(i));
    end
end
%plot(t,F)
contour(d,t,F,[0 0])
xlabel('d') 
ylabel('t')

%surf(t,d,F)
function F = f_calc(t, d) 
global k0 h_bar ksi m C_2
nD = floor(375/(2*pi.*t*1.2) - 0.5);
F = 0;
for k = 0:nD
    F = F + 1/(2*k+1).*imag(f_lg(k,t,d)+1i.*d.*k0.*((f_p1(k,t)-f_p2(k,t))./2)+1i*f_arg_1(k,t,d)-1i*f_arg_2(k,t,d));
end
F = -(1/d).*F - 7.826922968141167;
end
function p1 = f_p1(n,t)
p1 = ((1+1i)./sqrt(2)).*sqrt(t.*(2.*n+1));
end
function p2 = f_p2(n,t)
p2 = sqrt(3601+1i.*t.*(2.*n+1));
end 
function arg_1 = f_arg_1(n,t,d)
global k0
arg_1 = angle(1+exp(-1i.*d*k0.*f_p1(n,t)));
end
function arg_2 = f_arg_2(n,t,d)
global k0
arg_2 = angle(1+exp(-1i.*d*k0.*f_p2(n,t)));
end
function n_lg = f_lg(n,t,d)
global k0;
arg_of_lg = (1+exp(-1i.*d.*k0.*f_p1(n,t)))/(1+exp(-1i.*d.*k0.*f_p2(n,t)));
n_lg = log(abs(arg_of_lg));
end