# Increasing the bounds of integration on a strictly positive function decreases the result

1 view (last 30 days)
Nahman on 16 Aug 2019
Edited: Matt J on 16 Aug 2019
Here is the code that I am running
lambda = 40.8382-1i*6.4807e-15;
L0 = 4.5866e-8;
integrand = @(z) airy(0,-z./L0-lambda).*conj(airy(0,-z./L0-lambda));
lowerbound = -1;
C = integral(integrand,lowerbound,0,'ArrayValued',true)
As you can see I am integrating the function "integrand", which is a complex valued airy function times its conjugate, meaning it is strictly positive. C in this case is 9.3364e-08. When changing the bound of lower integration to -0.01, C becomes 1.9183e-08, which is understandably lower. However, setting lowerbound = -0.001 gives 9.3364e-08 again. Now I know most of the weight of this function is between 0 and about 1e-5 so I suspect that 9.3364e-08 is the right answer, but that there are specific points where it does not evaluate properly. Running this loop:
for i = 1:1000
lowerbound = -1/(i);
C(i) = integral(integrand,lowerbound,0,'ArrayValued',true);
end
figure
plot(C)
which produces this seems to confirm my suspicions. Would anyone know why this would be the case?

Walter Roberson on 16 Aug 2019
My tests suggest that the partial results are numerically insignificant below about -2.4e-6 (-1/-416666) where the value rises to about 9.3e-19 . The peak appears to be about -1.825e-6 where it reaches about 0.28
Walter Roberson on 16 Aug 2019
Specifying 'abstol' 1e-12 or below helps a lot. You will still get 4 inexplicable peaks, but the absolute differences are in the range 1e-17
Nahman on 16 Aug 2019
You are right, I did not know about setting the AbsTol but it seems that I can keep the errors at a reasonable for this specific usecase.

Matt J on 16 Aug 2019
Edited: Matt J on 16 Aug 2019
A plot of the function reveals highly oscillatory behavior in the interval [-2e-6,0] which then flattens out to zero in [-1,-2e-6].
fplot(integrand,[-2.5e-6,0],'MeshDensity',1e2) To discretize the integral reliably, you need to narrow the calculation to the sub-interval [-2e-6,0] and ensure a discretization interval that is finer than the oscillations. Below, I compare the result of integral() against trapz() with two levels of discretization fineness,
for i=1:10
lowerbound = -2.5e-6/i;
zs=linspace(lowerbound,0,1e5);
Ctrapz1(i)=trapz(zs,integrand(zs))
zs=linspace(lowerbound,0,1e6);
Ctrapz2(i)=trapz(zs,integrand(zs))
C(i) = integral(integrand,lowerbound,0,'ArrayValued',true)
end
plot([C-Ctrapz1;C-Ctrapz2].'); legend('C-Ctrapz1','C-Ctrapz2') 