# Plot a multi-variable equation in (discontinuous) regions

The goal is to display the relation between K and E in the following equation:

$$\cos(K) = 5 \text{sinc}( 5.12\sqrt{E} ) - \cos(5.12\sqrt{E} )$$

The intended result is:

However, using the following code:

k = linspace(0, 3*pi, 99);  % 0 < k < 3pi , crystal momentum.
E = linspace(0.1,  3.5, 99);% Energy.
[X, Y] = meshgrid(k, E);

P = 5;                      % 2mU_0b/hbar^2.
Z = cos(X)  - ( P * sinc( 5.12 * sqrt(Y) ) )  - cos( 5.12 * sqrt(Y) );

contour(X, Y, Z);

xlabel('k')
ylabel('E(k)')
legend('cos(K) - P * sin(E) - cos(E)')
title('Dispersion relation for periodic potential in extended zone.')


what I currently get unfortunately bears no resemblance at all:

I just don't see how to separate the graph into three (discontinuous) regions, $K = [0, \pi], [\pi, 2\pi], [2\pi, 3\pi]$, and what value of $E$ should be used for each of these regions.

How to plot the equation the same way as in the first graph?

The first plot is from here, where they do some modular arithmetic trick by adding $i*floor(k ./ \pi)$ on the l.h.s. and $i*floor(5.12*\sqrt{E} ./ \pi)$ on the r.h.s. There is something similar here.

• I tried to find the solution using bisection and am now convinced that the plot and the closed-form dispersion relation do not correspond to each other. Highlight the zero level-set in your contour plot and you'll see why, e.g., v = [0,0]; contour(X,Y,Z,v). – Biswajit Banerjee Sep 14 '18 at 22:59

First, you originally wrote that your equation is $\cos(K)=5 \text{ sinc}(5.12\sqrt{E})- \cos(5.12\sqrt{E})$, but you clearly meant $\cos(K)=5 \text{ sinc}(5.12\sqrt{E}) + \cos(5.12\sqrt{E})$

Second, the $\text{sinc}()$ function has differing conventions. The other question is based on Mathematica, which uses the convention $\text{sinc}(x) \equiv \frac{\sin x}{x}$, while Matlab/Octave use the convention $\text{sinc}(x) \equiv \frac{\sin (\pi x)}{\pi x}$. (To be precise of course we need to separately define $\text{sinc}(0)\equiv 1$.) This is why your first contour plot shows "no resemblance at all".

Compensating for these differing conventions and altering your code slightly so that it only shows the contours for Z=0 (as suggested in a comment by Biswajit Banerjee) yields

k = linspace(0, 3*pi, 99);  % 0 < k < 3pi , crystal momentum.
E = linspace(0.1,  3.5, 99);% Energy.
[X, Y] = meshgrid(k, E);

P = 5;                      % 2mU_0b/hbar^2.
Z = cos(X)  - ( P * pi * sinc( 5.12 * sqrt(Y)/pi ) )  - cos( 5.12 * sqrt(Y) );
contour(X, Y, Z,[0,0],'linewidth',2,'linecolor','k');
for i=1:3
line(i*[pi,pi],[0,3.5],'color','k','linestyle','--')
end

xlabel('k')
ylabel('E(k)')
legend('cos(K) - P * sin(E) - cos(E)')
title('Dispersion relation for periodic potential in extended zone.')


which produces the plot

The best way I could come up with to isolate sections of these curves is to simply run contour for the separate regions and plot them together. The following code accomplishes this and generalizes it to N sections.

P = 5;                      % 2mU_0b/hbar^2.
a=5.12;
N=3;

hold on;

for i=1:N
k = linspace((i-1)*pi, i*pi, 1000);  % 0 < k < 3pi , crystal momentum.
E = linspace(((i-1)*pi/a)^2,  (i*pi/a)^2, 1000);% Energy.
[X, Y] = meshgrid(k, E);

Z = cos(X)  - ( P * pi * sinc( a * sqrt(Y)/pi ) )  - cos( a * sqrt(Y) );
contour(X, Y, Z,[0;0],'linewidth',2,'linecolor','k');

line(i*[pi,pi],[0,(N*pi/a)^2+0.2],'color','k','linestyle','--')
end

xlabel('k')
xlim([0,N*pi+0.2])
ylim([0,(N*pi/a)^2+0.2])
ylabel('E(k)')


Here is something near to the intended goal:

But the code is ugly:

k = linspace(0,  3*pi, 100);

hbar = 1;
m = 1;
E =  ( (hbar^2) .* (k.^2) ) ./  (2*m);

E(33:66) = E(33:66) + 3;
E(66:100) = E(66:100) + 6;

E(33) = E(32);
E(66) = E(65);

E(34) = NaN;
E(67) = NaN;

hold on
plot(k,  E, [pi pi], [0 60], '--k', [2*pi 2*pi], [0 60], '--k');
plot(k(34), E(33), 'ob', 'MarkerFaceColor','b',...
k(67), E(66), 'ob',  'MarkerFaceColor', 'b',...
k(34), E(35), 'ob', 'MarkerFaceColor', 'b',...
k(67), E(68), 'ob' , 'MarkerFaceColor', 'b');
hold off

xlabel('k')
ylabel('E(k)')
legend('E = hbar^2 k^2 / 2m', 'x = pi', 'x = 2pi')
title('Dispersion relation for periodic potential in extended zone.')

txt = '<- Band Gap';
text(k(35), E(33) +2, txt)