How to plot a numerical integral using the trapezoidal method

11 次查看（过去 30 天）
Mahdi Yavarian
Mahdi Yavarian 2017-11-22
编辑： David Goodmanson 2019-3-25
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.
  2 个评论
Walter Roberson
Walter Roberson 2019-3-24
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.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

David Goodmanson
David Goodmanson 2019-3-25
编辑：David Goodmanson 2019-3-25
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.

类别

在 Help Center 和 File Exchange 中查找有关 Creating, Deleting, and Querying Graphics Objects 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by