How to pass on arguments in the form of two grids and return a matrix the elements of which involve conditional statements?

Consider the following function of the two variables r and t which returns 5 if and 8 if :

function a = f(r, t)

a = 5 * (t < r) + 8 * (r < t);

end

I would like to evaluate this function on a grid. The following implementation does the job:

N = 4;

DT = 1/N;

t = linspace(DT/2, 1-DT/2, N)';

r = t + 0.3 * DT;

[rr, tt] = ndgrid(r, t);

a = f(rr, tt)

Now suppose that 5 and 8 are actually very expensive operations that should only be evaluated at the corresponding condition.

I tried the following naive code but I only obtained a scalar and not the desired matrix.

function a = f(r, t)

if (t < r)

a = 5;

else

a = 8;

end

end

I was wondering what could be missing in my implementation. Any help is highly appreciated. Thank you!

[idermann]

It would of course be possible in the Matlab script to call the function times using 2 for loops but this is obviously very slow.

Answer 1

[VBBV]

N = 4;
DT = 1/N;
t = linspace(DT/2, 1-DT/2, N)';
r = t + 0.3 * DT;

[rr, tt] = ndgrid(r, t);
a = f(rr, tt)

a = 1×16

5 5 5 5 5 5 5 5 5 5 8 8 8 8 8 8

 
function a = f(r, t)
    
    idx1 = (t < r);  % 
    
    a1 = 5*ones(numel(idx1(idx1==1)),1);
    
    idx2 = r < t;    
    a2 = 8*ones(numel(idx2(idx2==1)),1);    
    
    a = [a1.' a2.'];

end

[idermann]

Thanks for the answer. What one should get is rather

5     8     8     8
5     5     8     8
5     5     5     8
5     5     5     5

Reshaping what we get from the function you wrote gives something different. But I should be able to figure out the rest myself.

[VBBV]

This is easier method than before, and something that you actually want

N = 4;
DT = 1/N;
t = linspace(DT/2, 1-DT/2, N)';
r = t + 0.3 * DT;

[rr, tt] = ndgrid(r, t);
a = f(rr, tt)

a = 4×4

5 8 8 8 5 5 8 8 5 5 5 8 5 5 5 5

 
 
function a = f(r, t)
    
    idx1 = (t < r);  % 
    A = idx1*5;    % this is easier than before
    
    idx2 = r < t;   
    B = idx2*8;    
    
    a = A+B;

end

Answer 2

[Dyuman Joshi]

"I tried the following naive code but I only obtained a scalar and not the desired matrix"

You obtained a scalar because you assigned a scalar.

And you are using non-scalar value as a condition to if-else statements, which evaluates all the values to a condition and proceeds accordingly.

if [0 1; 1 2]
    disp('if')
elseif 3*ones(5)
    disp('elseif')
else
    disp('else')
end

elseif

It depends on the operation but this should work for a good amount of cases. If this doesn't work, please specify the 'expensive operation' you are trying to implement.

function a = f(r, t)

id1 = r>t;
%do expensive operation on these indices

id2 = r<t; %or id2 = ~id1
%do expensive operation on these indices

end

[idermann]

Thanks for the answer. How to do expensive operations for these indices? Suppose that for r > t the function returns some complicated function of r and t that we denote g1(r,t). We would like to calculate g1(r,t) only when we are sure that r > t.

Also, for r < t, the function returns some complicated function of r and t that we denote g2(r,t). How to specify that already in the definition of the function f so as to avoid computing both g1 and g2 for all the values of r and t? I can simply make f return g1 .* (r > t) + g2 .* (r < t) but this would mean that the function f would have to compute g1 when r < t and g2 when r > t (very expensive operations). Any help is very much appreciated. Thank you,

[Dyuman Joshi]

Could you provide an example for g1 and g2? Or are you asking for a general case?

[idermann]

g1 and g2 involve infinite sums of terms involving hypergeom functions. They are really very complicated and there is no real interest in providing their expressions explicitly here I guess. If we could have a solution for general g1 and g2, that would be a dream. Thanks Dyuman.

[idermann]

I realized that using parfor for parallel computation leads to a quick evaluation of the kernel functions without the need for vectorizing the implementation. Sometimes simple approaches save both time and energy. Thanks for your valuable help anyway.