Increasing the bounds of integration on a strictly positive function decreases the result
2 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
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?
5 Kommentare
Walter Roberson
am 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
Akzeptierte Antwort
Matt J
am 16 Aug. 2019
Bearbeitet: Matt J
am 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')
0 Kommentare
Weitere Antworten (0)
Siehe auch
Kategorien
Mehr zu Numerical Integration and Differentiation finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!