Increasing the bounds of integration on a strictly positive function decreases the result
2 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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 commentaires
Walter Roberson
le 16 Août 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
Réponse acceptée
Matt J
le 16 Août 2019
Modifié(e) : Matt J
le 16 Août 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 commentaires
Plus de réponses (0)
Voir également
Catégories
En savoir plus sur Numerical Integration and Differentiation dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!