using QUADGK vs QUADL numerical integration algorithm
2 visualizzazioni (ultimi 30 giorni)
Mostra commenti meno recenti
Hi all, I have some issues while using quadgk with handle functions. Sometimes it gives me some warnings, where quadl function does not give me any warning at all! I can not understand why. My handle function is like: y(x) = x The integration is from -0.5 to 0.5
Why does it happen that?
0 Commenti
Risposta accettata
Mike Hosea
il 23 Dic 2011
Yes, you're on the right track. The problem here is that you are computing in finite precision and you don't have enough of it for the accuracy you are requesting. Let's suppose we integrate f(x) = x*2e11 from -0.05 to 0.05. The exact result is 0, so let's pretend we obtain this result by calculating (0.05 - -0.05)*(f(-0.05)+f(0.05))/2, i.e. 0.1*(-1e10+ 1e10)/2 = 0. Further suppose we have just one bit of roundoff error instead, so that we end up computing 0.1(-1e10+eps(-1e10) + 1e10)/2 = 0.05*eps(1e10) = 9.536743164062501e-08 . The default tolerances for QUADGK are reltol = 1e-6 and abstol = 1e-10. A single bit of roundoff error results here in an absolute error of 9.5e-8 and of course infinite relative error. This is numerically hopeless. You're essentially asking for the exact answer without so much as one bit of roundoff error anywhere. If you loosen the absolute tolerance to 1e-7 (quad and quadl have even a looser default tolerance than that), there is no warning:
q(1) = quadgk(@(z) funz(z,h,a,b,c,p,Elm,Elc,nim,nic,1,1),-.5*h,.5*h,'abstol',1e-7);
q(2) = quadgk(@(z) funz(z,h,a,b,c,p,Elm,Elc,nim,nic,2,1),-.5*h,.5*h,'abstol',1e-7);
q(3) = quadgk(@(z) funz(z,h,a,b,c,p,Elm,Elc,nim,nic,3,1),-.5*h,.5*h,'abstol',1e-7);
>> testquad
ans =
1.0e-07 *
-0.111758708953857 -0.158324837684631 -0.018626451492310
0 Commenti
Più risposte (3)
the cyclist
il 21 Dic 2011
Neither of these bits of code cause a problem for me. Did you do something differently? Maybe you could post the code you ran and the warnings you received.
% quadgk
f1 = @(x) x;
Q1 = quadgk(f1,-0.5,0.5)
% quadl
f2 = @(x) x;
Q2 = quadl(f2,-0.5,0.5)
EDIT IN RESPONSE TO ADDITIONAL INFO:
Here is a drastically simplified version of your code that exhibits the warning you saw.
% Function to integrate
funz = @(z) 1.e11 * z;
% Plot variable over range that is going to be integrated
z = -0.05:0.01:0.05;
plot(z,funz(z),'.-')
% Integrate
q = quadgk(@(z) funz(z),min(z),max(z))
I think the essence is that you are taking an integral of a function that is of order 10^11, but has an integral that is zero. In the course of MATLAB estimating the error in that integral (in order to determine whether it should stop making smaller intervals, and return a value), you are getting values that are tiny in comparison to the function itself.
That might not be exactly correct, but it is definitely the gist of the warning.
I have not gone back and looked at the ramifications for your actual function. However, my guess is that you could pull out some large scaling factor while doing the integral, and put it back in after the integration.
0 Commenti
Nicholas
il 21 Dic 2011
2 Commenti
the cyclist
il 21 Dic 2011
I'll try to take a look, but it might not be soon. In the meantime, a couple suggestions:
One is that if you have not done so, trying plotting what your functions look like over the interval that you are trying to integrate. Sometimes that exposes some unexpected weirdness. The other, related to the use of this forum, is that you might want to delete this "answer" (which is not really an answer) and incorporate the code into your question itself.
Mike Hosea
il 23 Dic 2011
I'm not sure what notifications the poster gets. Since I added an "I'll take a look" answer first and later edited it, maybe I should add an "I took a look" answer so that the OP will get notified.
0 Commenti
Vedere anche
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!