Area difference in using trapz() through a loop and without a loop.

조회 수: 5 (최근 30일)
Alvin
Alvin 2017년 11월 4일
다시 열림: Walter Roberson 2017년 11월 5일
clear
C1 = 0;
k1 = 1;
k2 = 0;
n1 = 0;
n2 = 0;
nm = 1;
D = 0;
O = 0;
gam = 1e-04.*k1;
alph = 0;
w = linspace(-5,5,1e+06);
delw = w(2)-w(1);
TA = (w+D)./(k1./2);
TB = (w-O)./(gam/2);
xa = 1-1i*TA;
xb = 1-1i*TB;
C = (2.*1i.*sqrt(C1))./(sqrt(gam).*((xa.*xb)+C1));
D = (2.*xa)./(sqrt(gam).*((xa.*xb)+C1));
MC = (abs(C)).^2;
MD = (abs(D)).^2;
Sbb = 2.*pi*(MC.*(n1+0.5)+MD.*(nm+0.5));
A = delw.*trapz(Sbb)./(2.*pi);
figure(1)
plot(w,Sbb);
set(gca,'FontSize',13)
Upon compiling the code above, one can observe that my area under the curve of Sbb versus w as given by A = delw.*trapz(Sbb)./(2.*pi) gives A = 9.4247. This setting is at C1 = 0. If I change my C1 to C1 = 1, my new A value gives A = 6.2834 (try it!). All is good.
Now, I know that A is dependent on C1, so I tried to loop C1 over a range of values like the following (it looks long but the code is almost exactly the same but with minor tweaks) to observe the trend:
clear
%C1 = 0:1:1;
C1 = linspace(0,500,1000);
k1 = 1;
k2 = 0;
n1 = 0;
n2 = 0;
nm = 1;
D = 0;
O = 0;
gam = 1e-04.*k1;
alph = 0;
w = linspace(-5,5,1e+06);
delw = w(2)-w(1);
for ii=1:length(C1)
TA = (w+D)./(k1./2);
TB = (w-O)./(gam/2);
xa = 1-1i*TA;
xb = 1-1i*TB;
C = (2.*1i.*sqrt(C1(ii)))./(sqrt(gam).*((xa.*xb)+C1(ii)));
D = (2.*xa)./(sqrt(gam).*((xa.*xb)+C1(ii)));
MC = (abs(C)).^2;
MD = (abs(D)).^2;
Sbb = 2.*pi.*(MC.*(n1+0.5)+MD.*(nm+0.5));
A = delw.*trapz(Sbb)./(2.*pi);
AA(ii)=A;
end
figure(2)
plot(C1,AA);
set(gca,'FontSize',13)
What I get is a bunch of zagged lines which does not seem to fit the bill. My first motivation was to switch the range of my C1 so that they are in increments of one unit and just up to 1 and recompile it:
C1 = 0:1:1;
%C1 = linspace(0,500,1000);
If I look over at my workspace, at C1 = 1, my AA(2) = 9.6485 which is clearly not 6.2834 as computed without the loop (but still at C1=1). I am absolutely convinced that 6.2834 is the correct value and 9.6485 isn't (Since A should decrease with increasing C1). What is the problem here? I would greatly appreciate any help with this.
Thanks in advance!
  댓글 수: 2
John D'Errico
John D'Errico 2017년 11월 4일
편집: John D'Errico 2017년 11월 4일
You accepted the answer the last time you asked virtually the identical question. What was wrong with that answer?
In fact, that answer was correct to your last question. That you do not understand the answer is no reason to ask the same question multiple times.
Alvin
Alvin 2017년 11월 4일
I did not say that anything is wrong with the answer to my last question. In fact the answer made sense. This is a different question looping over a particular parameter albeit in the same code structure as before.
I suggest you re-read my last post. If you've read it with proper scrutiny, you'd noticed the difference.

댓글을 달려면 로그인하십시오.

답변 (0개)

카테고리

Help CenterFile Exchange에서 Numerical Integration and Differentiation에 대해 자세히 알아보기

태그

Community Treasure Hunt

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

Start Hunting!

Translated by