Trouble when using dblquad for product of functions

5 Nov 2013
1 Answer
Answer Accepted
Updated 8 Nov 2013
2 Views (30 days)

0 votes

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
Show 3 older comments Hide 3 older comments

Mike Hosea
Mike Hosea on 6 Nov 2013
Edited: Mike Hosea on 6 Nov 2013
You've provided enough information for me to compare the integrand when var = 1, i = 1, j = 2 to var = 1, i = 2, j = 1, but I can't see what the corresponding var = 2 versions look like.
Note that DBLQUAD is obsolete. The new routine is called INTEGRAL2. To use it you will need to reformulate your call sites a little bit. Instead of
dtemp(i,j) = dblquad('integrand',(0+(m)*h),(0+(m+1)*h),(0+(m)*h),(0+(m+1)*h),[],@QUADL,var,i,j,h);
you would write
dtemp(i,j) = integral2(@(x,y)integrand(x,y,var,i,j,h),(0+(m)*h),(0+(m+1)*h),(0+(m)*h),(0+(m+1)*h));
Jonathan
Jonathan on 6 Nov 2013
Edited: Jonathan on 6 Nov 2013
I typed help integral2 in my command line and integral2 is not found.
VAR = 1 means that the first function is a function of x
VAR = 2 means that the first function is a function of y
here is the entire integrand
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
if (i == 3)
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 == 4)
if (j == 1)
g = (-2*x/h + 3*x.^2/h^2).*( 6*x/h^2 + 6*x.^2/h^3).*H;
end
if (j == 2)
g = (-2*x/h + 3*x.^2/h^2).*(1 - 4*x/h + 3*x.^2/h^2).*H;
end
if (j == 3)
g = (-2*x/h + 3*x.^2/h^2).*(6*x/h^2 - 6*x.^2/h^3).*H;
end
if (j == 4)
g = (-2*x/h + 3*x.^2/h^2).*(-2*x/h + 3*x.^2/h^2).*H;
end
end
end
if (var == 2)
if (i == 1)
if (j == 1)
g = (- 6*y/h^2 + 6*y.^2/h^3).*( 6*x/h^2 + 6*x.^2/h^3).*H;
end
if (j == 2)
g = (- 6*y/h^2 + 6*y.^2/h^3).*(1 - 4*x/h + 3*x.^2/h^2).*H;
end
if (j == 3)
g = (- 6*y/h^2 + 6*y.^2/h^3).*(6*x/h^2 - 6*x.^2/h^3).*H;
end
if (j == 4)
g = (- 6*y/h^2 + 6*y.^2/h^3).*(-2*x/h + 3*x.^2/h^2).*H;
end
end
if (i == 2)
if (j == 1)
g = (1 - 4*y/h + 3*y.^2/h^2).*( 6*x/h^2 + 6*x.^2/h^3).*H;
end
if (j == 2)
g = (1 - 4*y/h + 3*y.^2/h^2).*(1 - 4*x/h + 3*x.^2/h^2).*H;
end
if (j == 3)
g = (1 - 4*y/h + 3*y.^2/h^2).*(6*x/h^2 - 6*x.^2/h^3).*H;
end
if (j == 4)
g = (1 - 4*y/h + 3*y.^2/h^2).*(-2*x/h + 3*x.^2/h^2).*H;
end
end
if (i ==3)
if (j == 1)
g = (6*y/h^2 - 6*y.^2/h^3).*( 6*x/h^2 + 6*x.^2/h^3).*H;
end
if (j == 2)
g = (6*y/h^2 - 6*y.^2/h^3).*(1 - 4*x/h + 3*x.^2/h^2).*H;
end
if (j == 3)
g = (6*y/h^2 - 6*y.^2/h^3).*(6*x/h^2 - 6*x.^2/h^3).*H;
end
if (j == 4)
g = (6*y/h^2 - 6*y.^2/h^3).*(-2*x/h + 3*x.^2/h^2).*H;
end
end
if (i == 4)
if (j == 1)
g = (-2*y/h + 3*y.^2/h^2).*( 6*x/h^2 + 6*x.^2/h^3).*H;
end
if (j == 2)
g = (-2*y/h + 3*y.^2/h^2).*(1 - 4*x/h + 3*x.^2/h^2).*H;
end
if (j == 3)
g = (-2*y/h + 3*y.^2/h^2).*(6*x/h^2 - 6*x.^2/h^3).*H;
end
if (j == 4)
g = (-2*y/h + 3*y.^2/h^2).*(-2*x/h + 3*x.^2/h^2).*H;
end
end
end
end
Mike Hosea
Mike Hosea on 6 Nov 2013
Edited: Mike Hosea on 6 Nov 2013
Unless your version is very old, you can probably use QUAD2D instead of INTEGRAL2 for this problem.
for j = 1:4
for i = 1:4
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));
dtemp(i,j) = dtemp(i,j) - quad2d(@(x,y)integrand(x,y,2,i,j,h),(0+(m)*h),(0+(m+1)*h),(0+(m)*h),(0+(m+1)*h));
end
end
Jonathan
Jonathan on 8 Nov 2013
What is the most accurate I can make these integration schemes?
Mike Hosea
Mike Hosea on 8 Nov 2013
Edited: 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.

Sign in to answer this question.

 Accepted Answer

Mike Hosea
Mike Hosea on 6 Nov 2013

0 votes

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)

Categories

Find more on Data Distribution Plots in Help Center and File Exchange

Tags

Asked:

Jonathan
Jonathan
on 5 Nov 2013

Edited:

Mike Hosea
Mike Hosea
on 8 Nov 2013

Community Treasure Hunt

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

Start Hunting!