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

1 view (last 30 days)
idermann on 27 Mar 2023
Commented: VBBV on 27 Mar 2023
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 on 27 Mar 2023
Edited: idermann on 27 Mar 2023
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.

VBBV on 27 Mar 2023
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 on 27 Mar 2023
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 on 27 Mar 2023
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

Dyuman Joshi on 27 Mar 2023
Edited: Dyuman Joshi on 27 Mar 2023
"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 on 27 Mar 2023
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 on 27 Mar 2023
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.

Categories

Find more on Parallel Computing Fundamentals in Help Center and File Exchange

R2022b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!