Fundamental issue with trapz() function (or not?)

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,300);
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 = trapz(w,Sbb)./(2.*pi);
%AA(ii)=A;

figure(1)
plot(w,Sbb);
set(gca,'FontSize',13)
xlabel('\omega')
ylabel('S_{bb}')
%xlim([0,10000])

Upon compiling the code above, one would find a curve with a very narrow width centered about 0. Furthermore, it can be seen that the area under the curve was computed using the trapz() function and stored in A (well it's divided by 2pi but that is just scaling). A = 0.0442 in this case as observed from the workspace and all is well.

However as I decrease the range of my w. Say:

w = linspace(-0.1,0.1,300);

My area would give A = 2.1705 in this case (it increased). This doesn't make sense given that the tail ends on both sides of the curve barely contributes to the area under the curve (it's a spectral density curve) since they are so small. It is implying that as my range (w) decreases, my area increases which is physically absurd. Furthermore, upon some investigation, it seems that the area would peak roughly when the range is at:

w = linspace(-8e-03,8e-03,300)

Which gives A = 9.3343. Further reducing the range would then decrease the area. Does anyone know what seems to be the issue?

Thanks in advance!

Answer 1

[David Goodmanson]

Hi Alvin,

--- Modified Answer ----

The width of the peak is around 1e-4. So for trapz to work, you need the spacing in the w array to be significantly less than that. Using linspace changes the spacing, and for both the (-5,5) cases and the (-.1,1) cases the spacing is not fine enough. If you use a small, constant spacing then

w = -5:1e-6:5
A = 9.4247

w = -.1:1e-6:.1
A = 9.4218

[Alvin]

Thanks! For some reason I assumed that trapz() takes care of the non-unit spacing. Confirmed it myself again that this was not the case after looking it up!

[Alvin]

Edit: I double checked and the issue doesn't seem to be resolved. According to the Matlab documentation on trapz, if I have a uniform spacing on my x-axis of length delw then the Area would be given as:

A = delw*trapz(Sbb)./(2.*pi)

Or just

 A = trapz(w,Sbb)./(2.*pi)

These two are equivalent. However

 A = delw*trapz(w,Sbb)./(2.*pi)

is not the same as the previous two (you can check it yourself). The third equation is simply the scaling of the second equation by a factor delw. As such I think I evaluated my A correctly and the issue is still unresolved.

Apologies for not checking more thoroughly.

[David Goodmanson]

Hi Alvin,

you're right, putting w into the argument of trapz should take care of delw, so there is something else happening. See the modified answer above.

[Alvin]

I've checked thoroughly and this seems to be what I'm looking for. Thanks for the correction!