correct function handle for integral2 in matlab

Question

I created a function in matlab that returns a vector like

function w = W_1D(x,pos,h)
w=zeros(1,length(x));
        if (h~=0)
            xmpos = x-pos; 
            inds1 = (-h <= xmpos) & (xmpos < 0);
            w(inds1) = xmpos(inds1)./h + 1;
            inds2 = (0 <= xmpos) & (xmpos <= h);
            w(inds2) = -xmpos(inds2)./h + 1;
        else
            error('h shouldn't be 0')
        end
 end

Thus, in the end, there is a vector w of size length(x). Now i created a second function like

function f = W_2D(x,y,pos_1,pos_2,h)
  w_x = W_1D(x,pos_1,h);
  w_y = W_1D(y,pos_2,h);
  f = w_x'*w_y;
end

where length(x)=length(y). Thus, the function W_2D obviously returns a matrix. But when I now try to evaluate the integral over a rectangular domain like e.g.

V = integral2(@(x,y) W_2D(x,y,2,3,h),0,10,0,10);

matlab returns some errors:

Error using integral2Calc>integral2t/tensor (line 242)
Integrand output size does not match the input size.

Error in integral2Calc>integral2t (line 56)
[Qsub,esub] = tensor(thetaL,thetaR,phiB,phiT);

Error in integral2Calc (line 10)
[q,errbnd] = integral2t(fun,xmin,xmax,ymin,ymax,optionstruct);

Error in integral2 (line 107)
Q = integral2Calc(fun,xmin,xmax,yminfun,ymaxfun,opstruct);

I also tried to vary something in the W_2D-function: instead of f = w_x'*w_y; I tried f = w_x.'*w_y; or w_y = transpose(w_y); f = kron(w_x,w_y);, but there is always this error with the Integrand output size-stuff. Can anyone explain, where my fault is?

EDIT: After Werner's hint with the keyboard debugging method, I can tell you the following. The first step returns w_x of type <1x154 double>, w_y is <1x192 double>, x and y are both <14x14 double>. In the next step, f appears with a value of <154x192 double>. Then everything disappears, except x and y and the matlab-function integral2Calc.m appears in the editor and it jumps to the Function Call Stack integral2t/tensor and after some more steps, the error occurs here

    Z = FUN(X,Y);  NFE = NFE + 1;
    if FIRSTFUNEVAL
        if ~isfloat(Z)
            error(message('MATLAB:integral2:UnsupportedClass',class(Z)));
        end
        % Check that FUN is properly vectorized. This is important here
        % because we (otherwise) always pass in square matrices, which
        % reduces the probability of the user generating an error by
        % using matrix functions instead of elementwise functions.
        Z1 = FUN(X(VTSTIDX),Y(VTSTIDX)); NFE = NFE + 1;
        if ~isequal(size(Z),size(X)) || ~isequal(size(Z1),size(VTSTIDX))
            % Example:
            % integral2(@(x,y)1,0,1,0,1)
            error(message('MATLAB:integral2:funSizeMismatch'));
        end

Hope that information is detailed enough...I have no idea what happenes, because my example is exact as it is given on the mathworks site about integral2, isn't it?

Maybe I should precise a bit more, what I wanna do: since W_2D gives me a surface w(x,y) of a compactly supported 2-dimensional hat-function, stored in a matrix w, I want to calculate the volume between the (x,y)-plane and the surface z=w(x,y)...

EDIT2: I still do not understand how to handle the problem, that integral2 creates matrices as inputs for my W_1D-functions, which are called in W_2D and intended to have a <1xn double>-valued input and return a <1xn double> output, but at least I can simply use the following to solve the integration over the tensor product by using two one-dimensional integral-calls, that is

V = integral(@(x)integral(@(y)W_1D(y,3,h),0,10).*W_1D(x,2,h),0,10);

Answer 1

This first function is quite wrong. You are not indexing the array positions while you are doing w = x inside for.

Besides, if that would work, you are returning a line vector, that is, size 1xlength(x) and when you do w_x'*w_y you are doing length(x)x1 times 1xlength(y), which would give you a matrix length(x)*length(y).

Consider correcting your function:

function w = W_1D(x,pos)
   w = zeros(length(x),1); % Allocate w as column vector, so that the product gives a scalar (as I suppose that it is what you want.
   for ii=1:length(x) % Here, so that is indexes w and x elements as you need
       w(ii)=x(ii) - pos; % I changed your code to something that makes sense, but I don't know if that is what you want to do, you have to adapt it to work correctly.
   end
end

You may also want to debug your functions, consider adding keyboard before your operations and check what they are returning using dbstep. I.e:

function f = W_2D(x,y,pos_1,pos_2)
  w_x = W_1D(x,pos_1);
  w_y = W_1D(y,pos_2);
  keyboard
  f = w_x'*w_y;
end

Execution will stop at keyboard, then you can check w_x size, w_y size, and do dbstep to go after f = w_x'*w_y and see what it returned. After you finish debug, you can do dbcont so that it will continue execution.

This answer is a draft as it is quite difficult to help you with the information you have provided. But I think you can start working the things out with this. If you have more doubts feel free to ask.