Integrate a function that doesn't directly operate element-wise

Question

I'm working on a Matlab script that involves fourth order tensors calculations volume integrals. Let H(r,theta,phi) be the function I want to integrate. Assume that H cannot be obtained by simple operations on r, theta and phi.

My problem is that in Matlab as in any other code I know :

All input functions must accept arrays and operate elementwise. The function FUN(X,Y,Z)
must accept arrays X, Y, Z of the same size and return an array of corresponding values.

This is from the implementation of the integral3 function from Matlab.

If I try with this simple function :

fun = @(X,Y,Z) X.*Y.*Z

There is no problem at all and if I integrate it over [0,1] x [0,1] x [0,1], I get the right result :

integral3(fun,0,1,0,1,0,1)

returns 0.125 which is correct.

The problem is that as I said, I can't perform simple calculations with the vectors to obtain H and I am forced to do things more or less this way :

function [result] = fun(x,y,z)
    sz = length(x);
    result  = zeros(1,sz);
    for i=1:sz
        result(i) = x(i)*y(i)*z(i);
    end
end

This function works on its own and returns exactly the same results as the other one I introduced earlier. However, when I try to use integral3 I get this error :

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

But it can clearly be seen from the definition of my function that I specifically made it the size of the input.

I don't understand what's wrong and I'm not sure I have any other solution to compute this function than using this kind of syntax.

Thanks a lot for your time and help :)

Polemics

Nothing to add to the answer of @drhagen but I just want to thank you for this beautifully written question! – Bod

Answer 1

You are on the right track, but you are building an array with the correct length but the wrong size. It is a subtle difference in Matlab. I'm guessing that integral3 is passing in a column vector, but your function always returns a row vector. The column and row vectors have the same "length" sz, but different "sizes": the column vector is [sz,1] and the row vector is [1,sz]. The code below does what you want because it uses size to ensure that all the dimensions of the output match the input and numel to loop over the individual elements:

function result = fun(x,y,z)
    sz = size(x);
    result = zeros(sz);
    for i = 1:numel(x)
        result(i) = x(i)*y(i)*z(i);
    end
end

A good rule-of-thumb is to only use size and numel and never use length, which is a contender for worst function in Matlab.

Ventris

@Ventris It works when I use size and a double loop over size(1) and size(2) just in case the function integral3 uses an array instead of a column or a row vector, thank you ! Also, thank you for the length advice! – Polemics

@JeremyDiallo Actually, due to linear indexing, my version will work if supplied with an array, even if its an array that's 2D, 3D, 4D, etc. – Ventris

@Ventris Wow I didn't know about automatic linear indexing in matlab, I'd better look it up. Thanks for your help ! – Polemics