How to plot a numerical integral using the trapezoidal method
조회 수: 5 (최근 30일)
이전 댓글 표시
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
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
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.
댓글 수: 0
참고 항목
카테고리
Help Center 및 File Exchange에서 Annotations에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!