MATLAB Answers

Nahman
0

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

Asked by Nahman
on 16 Aug 2019
Latest activity Edited by 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
integral_failure.jpg
seems to confirm my suspicions. Would anyone know why this would be the case?

  5 Comments

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
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
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.

1 Answer

Answer by Matt J
on 16 Aug 2019
Edited by Matt J
on 16 Aug 2019
 Accepted Answer

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')
untitled.png

  0 Comments

Sign in to comment.