How to plot a numerical integral using the trapezoidal method

I want to plot a figure based on an equation. The equation in mind is composed of two terms. The first term is a simple equation, consisting of different parameters. But, the second equation with the same parameters, is an integral which has to be taken into account numerically. Of the numerical methods, I used the trapezoidal method. I'm working with the following code:

 clc
clear all
close all
e=1.60217662d-19;       
kB=1.3806488d-23;       
h=6.62607004d-34;       
hbar=h/(2*pi);          
T=300;                  
R=kB*T;
tau=0.0658d-12;
gamma=1/tau; 
vf=1d6;                 
omega=2*pi*1d12*(0:0.1:10);  
uc=e*[0.01,0.015,0.02];
for  k=1:length(omega)
    for  n=1:length(uc)
         D=(e^2)/(4*hbar);
         eps=e*(0:0.01:10);
         W1=1./(1+exp((eps-uc(n))/(R)));
         W2=1./(1+exp((-eps-uc(n))/(R)));
         W=W2-W1;
         E=(1i*(e^2))/(pi*hbar^2);
         sigma_inter(k,n)=trapz(eps,(E*(omega(k)+1i*gamma)*W)./((omega(k)+1i*gamma)^2-4*(eps/hbar).^2));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
         A=(uc(n)/R)+2*log(1+exp(-uc(n)/R));
         B=(1i*(e^2)*R)/(pi*hbar^2);
         C=omega(k)+1i*gamma;
        sigma_intra(k,n)=(A*B)./C;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        sigma(k,n)=(sigma_intra(k,n)+sigma_inter(k,n))/D;
    end
yyaxis left
plot(omega/(1d12*2*pi),imag(sigma(k,n)));
ylabel('Im[\sigma]');
yyaxis right
plot(omega/(1d12*2*pi),real(sigma(k,n)));
ylabel('Re[\sigma]');
xlabel('Frequency (THz)');
title('Graphene Surface Conductivity');
end

When plotting the imaginary and real parts of the equation, the figures are weird. As it seems, the dimensions of some parameters are not in agreement with each other. so I think, it's all about the dimensionality. But I have no idea how to resolve this problem.

[Mostafa alfarq]

wow what this huge code

[Walter Roberson]

We do not know what the graph "should" look like, so we have no basis to know if the output figures are weird or normal.

I do not see anything obviously wrong when I look at the output.

Answer 1

[David Goodmanson]

Hi Mahdi,

I believe your units are correct, although dividing by D gives normalized conductivity and not absolute conductivity.

I think the problems are caused by doing plotting within the for loops. Occasionally the autoplotting produces some funky choices when given free reign.

If you move the last 'end' to the location just above the plot commands and remove the one-at-a-time indexing [ sigma(k,n) ] in the plots, then the lines after the first end statement become

end
figure(1)
yyaxis left
plot(omega/(1d12*2*pi),imag(sigma));
ylabel('Im[\sigma]');
yyaxis right
plot(omega/(1d12*2*pi),real(sigma));
ylabel('Re[\sigma]');
xlabel('Frequency (THz)');
title('Graphene Surface Conductivity');

and you get what appear to be good results.