using QUADGK vs QUADL numerical integration algorithm

2011 12 21
4 个回答
回答已采纳
5 次查看（30 天）

0 个投票

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?

请先登录，再回答此问题。

 采纳的回答

Mike Hosea
Mike Hosea 2011-12-23

1 个投票

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

更多回答（3 个）

the cyclist
the cyclist 2011-12-21

0 个投票

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.
Nicholas
Nicholas 2011-12-21

0 个投票

Your are definitely right! In your proposed case there is no problem, in fact my function is a bit complicated and the problem occurs when some of its coefficients are zero and the function becomes linear.
Here's the code (I tried to write it as short as possible):
function [q] = testquad()
Elm = 70e9; Elc = 168e9; nim = 0.3; nic = 0.3; h = 0.1;
a = 1; b = 0; c = 0; p = 0;
q(1) = quadgk(@(z) funz(z,h,a,b,c,p,Elm,Elc,nim,nic,1,1),-.5*h,.5*h);
q(2) = quadgk(@(z) funz(z,h,a,b,c,p,Elm,Elc,nim,nic,2,1),-.5*h,.5*h);
q(3) = quadgk(@(z) funz(z,h,a,b,c,p,Elm,Elc,nim,nic,3,1),-.5*h,.5*h);
end
function f1 = funz(z,h,a,b,c,p,Elm,Elc,nim,nic,count,n)
f0 = ((1 - a.*(.5 + z./h) + b.*((.5 + z./h).^c)).^p); % fgm law
ni = ((nic - nim).*f0 + nim);
El = ((Elc - Elm).*f0 + Elm);
switch count
case 1 % Q11
Q = El./(1-ni.^2);
case 2 % Q12
Q = (ni.*El)./(1-ni.^2);
case 3 % Q66
Q = El./(2.*(1+ni));
end
f1 = Q.*z.^n;
end
Copy that in a new m-file and run it, you'll see this kind of error, one for each QUADGK function call
Warning: Reached the limit on the maximum number of intervals in use.
Approximate bound on error is 1.1e-09. The integral may not exist, or
it may be difficult to approximate numerically. Increase MaxIntervalCount
to 918 to enable QUADGK to continue for another iteration.
Thanks

2 个评论
显示 无 隐藏 无

the cyclist
the cyclist 2011-12-21
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
Mike Hosea 2011-12-23

0 个投票

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.

类别

在 帮助中心 和 File Exchange 中查找有关 Loops and Conditional Statements 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by