Trouble when using dblquad for product of functions

1 view (last 30 days)
Jonathan
Jonathan on 5 Nov 2013
Edited: Mike Hosea on 8 Nov 2013
Hello, I am using dblquad to integrate a function based on some parameters. Here is the call inside of a script file
for j = 1:4
for i = 1:4
dtemp(i,j) = dblquad('integrand',(0+(m)*h),(0+(m+1)*h),(0+(m)*h),(0+(m+1)*h),[],@QUADL,1,i,j,h);
dtemp(i,j) = dtemp(i,j) - dblquad('integrand',(0+(m)*h),(0+(m+1)*h),(0+(m)*h),(0+(m+1)*h),[],@QUADL,2,i,j,h);
end
end
The first couple of possible integrands from the 'integrand' function are
function [ g ] = integrand(x, y, var,i, j, h)
a = 1.040394;
b = 0.064362;
H = a/(b*sqrt(2*pi))*exp(-(x-y).^2/(2*b^2));
if (var == 1)
if (i == 1)
if ( j == 1)
g = ( 6*x/h^2 + 6*x.^2/h^3).*( 6*x/h^2 + 6*x.^2/h^3).*H;
end
if ( j == 2)
g = ( 6*x/h^2 + 6*x.^2/h^3).*(1 - 4*x/h + 3*x.^2/h^2).*H;
end
if ( j == 3)
g = ( 6*x/h^2 + 6*x.^2/h^3).*(6*x/h^2 - 6*x.^2/h^3).*H;
end
if ( j == 4)
g = ( 6*x/h^2 + 6*x.^2/h^3).*(-2*x/h + 3*x.^2/h^2).*H;
end
end
if (i == 2)
if (j == 1)
g = (1 - 4*x/h + 3*x.^2/h^2).*( 6*x/h^2 + 6*x.^2/h^3).*H;
end
if (j == 2)
g = (1 - 4*x/h + 3*x.^2/h^2).*(1 - 4*x/h + 3*x.^2/h^2).*H;
end
if (j == 3)
g = (1 - 4*x/h + 3*x.^2/h^2).*(6*x/h^2 - 6*x.^2/h^3).*H;
end
if (j == 4)
g = (1 - 4*x/h + 3*x.^2/h^2).*(-2*x/h + 3*x.^2/h^2).*H;
end
end
.
.
.
.
So the output matrix dtemp should ultimately be symmetric, however it is not, any ideas? Thank you for your time and help.
-Jonathan
  5 Comments
Mike Hosea
Mike Hosea on 8 Nov 2013
Well, QUAD2D is fast enough that you can probably afford to crank down the tolerances. Do this like
dtemp(i,j) = quad2d(@(x,y)integrand(x,y,1,i,j,h),(0+(m)*h),(0+(m+1)*h),(0+(m)*h),(0+(m+1)*h),'AbsTol',100*eps,'RelTol',100*eps);
Generally you want to set AbsTol to something that is small enough to be considered essentially zero. A reltol of 100*eps = 2.2e-14 should give you about 13 significant digits, give or take, unless the answer gets close to zero, in which case reltol doesn't matter, only abstol. (Conversely, if the answer is not close to zero, abstol doesn't matter very much, only reltol).

Sign in to comment.

Accepted Answer

Mike Hosea
Mike Hosea on 6 Nov 2013
Your integrand function is not defined to be symmetric. We have
dtemp(1,2) = dblquad('integrand',(0+(m)*h),(0+(m+1)*h),(0+(m)*h),(0+(m+1)*h),[],@QUADL,1,1,2,h) ...
- dblquad('integrand',(0+(m)*h),(0+(m+1)*h),(0+(m)*h),(0+(m+1)*h),[],@QUADL,2,1,2,h);
Let me declutter this by writing DOUBLEINTEGRAL to represent the integration with the stated limits. It turns out that
dtemp(1,2) = DOUBLEINTEGRAL(( 6*x/h^2 + 6*x.^2/h^3).*(1 - 4*x/h + 3*x.^2/h^2).*H) - ...
DOUBLEINTEGRAL((- 6*y/h^2 + 6*y.^2/h^3).*(1 - 4*x/h + 3*x.^2/h^2).*H)
and
dtemp(2,1) = DOUBLEINTEGRAL((1 - 4*x/h + 3*x.^2/h^2).*( 6*x/h^2 + 6*x.^2/h^3).*H) - ...
DOUBLEINTEGRAL((1 - 4*y/h + 3*y.^2/h^2).*( 6*x/h^2 + 6*x.^2/h^3).*H)
Why should these be the same? If we re-order the terms to make them more similar in form, we get
dtemp(1,2) = DOUBLEINTEGRAL(( 6*x/h^2 + 6*x.^2/h^3).*(1 - 4*x/h + 3*x.^2/h^2).*H) - ...
DOUBLEINTEGRAL((- 6*y/h^2 + 6*y.^2/h^3).*(1 - 4*x/h + 3*x.^2/h^2).*H)
dtemp(2,1) = DOUBLEINTEGRAL(( 6*x/h^2 + 6*x.^2/h^3).*(1 - 4*x/h + 3*x.^2/h^2).*H) - ...
DOUBLEINTEGRAL(( 6*x/h^2 + 6*x.^2/h^3).*(1 - 4*y/h + 3*y.^2/h^2).*H)
The second integral of each pair is different.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!