using QUADGK vs QUADL numerical integration algorithm

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?

Answer 1

[Mike Hosea]

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

Answer 2

[the cyclist]

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.

Answer 3

[Nicholas]

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

[the cyclist]

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.

[Nicholas]

How can I do that?

Answer 4

[Mike Hosea]

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.