MATLAB Answers

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);
which produces this
seems to confirm my suspicions. Would anyone know why this would be the case?


Show 2 older comments
Walter Roberson
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
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.

Sign in to comment.

Accepted Answer

Matt J
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].
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;
C(i) = integral(integrand,lowerbound,0,'ArrayValued',true)
plot([C-Ctrapz1;C-Ctrapz2].'); legend('C-Ctrapz1','C-Ctrapz2')


Sign in to comment.

More Answers (0)

Translated by