Trouble when using dblquad for product of functions

5 Nov. 2013
1 Respuesta
Respuesta aceptada
Actualizado a las 8 Nov. 2013
13 Visualizaciones (30 días)

0 votos

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 comentarios
Mostrar 3 comentarios más antiguos Ocultar 3 comentarios más antiguos

Mike Hosea
Mike Hosea el 6 de Nov. de 2013
Editada: Mike Hosea el 6 de Nov. de 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 el 6 de Nov. de 2013
Editada: Jonathan el 6 de Nov. de 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 el 6 de Nov. de 2013
Editada: Mike Hosea el 6 de Nov. de 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 el 8 de Nov. de 2013
What is the most accurate I can make these integration schemes?
Mike Hosea
Mike Hosea el 8 de Nov. de 2013
Editada: Mike Hosea el 8 de Nov. de 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).

Iniciar sesión para comentar.

Iniciar sesión para responder a esta pregunta.

 Respuesta aceptada

Mike Hosea
Mike Hosea el 6 de Nov. de 2013

0 votos

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.

Más respuestas (0)

Categorías

Más información sobre Data Distribution Plots en Centro de ayuda y File Exchange.

Etiquetas

Preguntada:

Jonathan
Jonathan
el 5 de Nov. de 2013

Editada:

Mike Hosea
Mike Hosea
el 8 de Nov. de 2013

Community Treasure Hunt

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

Start Hunting!

Translated by