MATLAB: Trouble when using dblquad for product of functions

dblqaudintregrationquadl

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

Best Answer

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.
Related Question