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

Nothing to add to the answer of @drhagen but I just want to thank you for this beautifully written question! — Toghe, Jun 22 '17 at 15:05

score 4 · Accepted Answer · answered Jun 22 '17 at 14:53

4

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.

answered Jun 22 '17 at 14:53

drhagen

8,331
8
53
82

@drhagen 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! – Experience111 Jun 23 '17 at 07:09
@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. – drhagen Jun 23 '17 at 09:51
@drhagen Wow I didn't know about automatic linear indexing in matlab, I'd better look it up. Thanks for your help ! – Experience111 Jun 26 '17 at 06:12

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

1 Answers1