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

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]

I confirm your findings on my system.

If you do

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;
lowerbounds = [-1./(1:1000), 0];

C = integral(integrand,lowerbound,0, 'waypoints', lowerbounds);
Cnw = integral(integrand,lowerbound,0);


for i = 1:length(lowerbounds)-1
    Cp(i) = integral(integrand,lowerbounds(i),lowerbounds(i+1), 'waypoints', lowerbounds);
    C(i) = integral(integrand,lowerbounds(i),0, 'waypoints', lowerbounds);
    Cnw(i) = integral(integrand,lowerbounds(i),0);
end
subplot(1,3,1);
plot(C)
title('integral with waypoints')
subplot(1,3,2);
plot(cumsum(Cp(1:end-1)))
title('cumsum of integral between waypoints')
subplot(1,3,3)
plot(Cnw)
title('integral without waypoints')

then you will see that the partial sums (subplot 2) are all exactly 0 as far as integral() is concerned, and you will also see that whether you use waypoints or not makes a difference.

If you use waypoints, then the spike near 200 occurs at -1/151 .

>> G = simplify(sym(integrand)); integral(integrand,-1/151,0), integral(integrand,-1/151,0,'waypoints',lowerbounds), vpa(int(G, -1/151, 0)), vpaintegral(G, -1/151, 0)

ans =

      9.33637318816886e-08

ans =

      4.37729538998267e-09

ans =
 
0.000000093363731881684444589887123401249
 
 
ans =
 
1.7008e-15669
>> vpa(integrand(sym(-1/151)))
 
ans =
 
1.6757183152487926083618414343337e-31756827

Notice that even vpa(int()) vs vpaintegral() shows a major difference

Numerically it looks likely that double precision integration should be finding essentially 0 for most of the range.

I suspect we are seeing failures of approximation when the values are so close to 0.

[Nahman]

I suspected approximation errors but was hoping there would be a systematic way to resolve them, oh well. Thanks a lot for your time!

[Walter Roberson]

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]

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]

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.

Answer 1

[Matt J]

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')