Trouble when using dblquad for product of functions

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

[Mike Hosea]

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]

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]

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]

What is the most accurate I can make these integration schemes?

[Mike Hosea]

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).

Answer 1

[Mike Hosea]

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.