How to code a certain part of this equation?

I am trying to write a code for solving population based equations (PBEs) using cell average technique (CAT) by Kumar et al. (2006). Link to the paper: PAPER

However, I am stucking with interpreting a certain part of the whole formulation (mathematical notation) which I have circled in red in the figure. Can anyone help me in correcting me or explaining it? My attempt is somewhat as follows:

clear 
close all;
clc;
% ======== User handle starts ==================
%% Initializing
v_max    =  100.0;   % Maximum size/volume
v_min    =  1e-5;    % Minimum size/volume
r        =  2;       % Geometric ratio
% v0       =  0.01;    % Initial average size
% N0       =  1;       % Initial total number
% tspan    =  [0 20];  % time span for integration

%% Generating discretized grid
disp('Creating the geometric grid..')
m     = ceil(log(v_max/v_min)/log(r)); % Number of pivots
mv    = m + 1; % Number of volume boundaries
vi    = zeros(mv,1);
vi(1) = v_min;

for i = 2:mv
    vi(i) = vi(i-1)*r; %cell edges
end

x = zeros(m,1); % Creating pivots
w = zeros(m,1); % Bin width

for i = 1:m
    x(i) = (vi(i) + vi(i+1))/2; %cell centers/pivots
    w(i) = vi(i+1) - vi(i); % width
end

p = 1; %fixed pivot
for j = 1: m
    for k = 1:j
        v = x(j) + x(k) ;
        for i = 1:m-1
            if v < x(i+1) && v >= x(i)
                itvl(p,:) = [x(j), x(k)] ;
                Z(i,j,k) = 1;
                p = p +1;
                break;
            end
        end
    end
end

p = 1; %CAT
for j = 1: m
    for k = 1:j
        v = x(j) + x(k) ;
        for i = 1:mv-1
            if v < vi(i+1) && v >= vi(i)
                itvl(p,:) = [x(j), x(k)] ;
                Z(i,j,k) = 1;
                p = p +1;
            end
        end
    end
end
return
%%%%%%%% i-th equation %%%%%%%%%%%%
for i = 1:m
    for j = 1:m
        for k = 1:j
            B(i) = (1 - (1/2)*Z(i,j,k)*dirac_delta(i,j))*beta(j,k)*N(j)*N(k) ;
        end
    end
end

Answer 1

[Torsten]

for i = 1:m-1
    xmhalf = x(i);
    xphalf = x(i+1);
    B_agg(i) = 0.0;
    for k=1:m
        for j=k:m
            if vi(j) + vi(k) < xphalf && vi(j) + vi(k) >= xmhalf
                B_agg(i) = B_agg(i) + (1-0.5*dirac_delta(j,k))*beta(j,k)*N(j)*N(k);
            end
        end
    end
end

Not sure about the loop limits - recheck !

[Angshuman Podder]

Thank you for your solution. Really appreciate it! I think it makes sense now. I modified your code a little and tried to implement it in the overall code. However, I am still running into a solution not matching the published results. I have outlined the entire algorithm here in the figure:

Can you please have a look at my attempt and point out errors which is as follows. Thanks in advance.

function [dNdt] = aggregation_CA(t,y,rcenter,redge)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Function to calculate aggregation rate 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N = y;
% Initialize
n    = length(rcenter) ; %no. of cell centers %one less than edge 
B    = zeros(n,1) ;
M    = zeros(n,1) ;
abar = zeros(n,1) ;
D    = zeros(n,1) ;
                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                    
x = rcenter;%(4/3)*pi*r.^3;
vi = redge;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Step 1: Calculate birth rates 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% a. Calculation of kernels
for i = 1: n
    for j = 1:n
     beta(i,j) = 1; %constant kernel %test problem 1
    end
end
 
% b. Calculation of the net rate of addition of particles to cell i by
% coagulation of particles in lower cells
p=1;
for i = 1:n
    xmhalf = vi(i);
    xphalf = vi(i+1);
    B_agg(i) = 0.0;
    for k=1:n
        for j=k:n
            if x(j) + x(k) < xphalf && x(j) + x(k) >= xmhalf
                B_agg(i) = B_agg(i) + (1-0.5*dirac_delta(j,k))*beta(j,k)*N(j)*N(k);
                %interval(p,:) = [x(k), x(j)] ;
                %del(p,:) = x(j) - x(k);
                p = p +1;
            end
        end
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Step 2: Calculate volume averages
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% a. Calculation of fluxes into cell i as a result of coagulations

for i = 1:n
    xmhalf = vi(i);
    xphalf = vi(i+1);
    M(i) = 0.0;
    for k=1:n
        for j=k:n
            if x(j) + x(k) < xphalf && x(j) + x(k) >= xmhalf
                M(i) = M(i) + (1-0.5*dirac_delta(j,k))*beta(j,k)*N(j)*N(k)*(x(j)+x(k));   
            end
        end
    end
end

% b. Average volume of newborn particles in the ith cell
for i = 1:n
    abar(i) = M(i)/B(i);    
end

for i = 1:n
    if B(i) == 0 % if no particles is there in a bin, avoid NaN error
        abar(i) = 1;
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Step 3: Calculate of birth modification
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Bmod = B;

Bzero = 0 ; %zero boundary at first cell

Bmod(1) = Bzero...%B(1)*lambdaminus(1,abar(1),xcenter)*heaviSide(xcenter(1)-abar(1)) ...
          + B(1)*lambdaplus(1,abar(1),x)*heaviSide(abar(1)-x(1));
 
Bmod(n) = B(n-1)*lambdaminus(n,abar(n-1),x)*heaviSide(x(n-1)-abar(n-1)) ...
          + Bzero;%B(n-1)*lambdaplus(n,abar(n-1),xcenter)*heaviSide(abar(n-1)-xcenter(n-1))  ;


for i = 2:n-1
    Bmod(i) = B(i-1)*lambdaminus(i,abar(i-1),x)*heaviSide(abar(i-1)-x(i-1)) + ...
        B(i)*lambdaminus(i,abar(i),x)*heaviSide(x(i)-abar(i)) ...
        + B(i)*lambdaplus(i,abar(i),x)*heaviSide(abar(i)-x(i)) +...
        B(i+1)*lambdaplus(i,abar(i+1),x)*heaviSide(x(i+1)-abar(i+1)) ;
end
%Bmod(n) = Bmod(n-1) ; %zero flux at last cell
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Step 4: Calculate of death terms
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for i = 1:n
    sum = 0;
    for k = 1: n
        sum = sum + beta(i,k)*N(k) ;
    end
    D(i) = N(i)*sum;
end

%dNdt = Bmod - D; %Birth-Death
for i =1:n
    dNdt(i,1) =  Bmod(i) - D(i);
end
fprintf('t= %2.2e; maxN = %2.2e ; minN = %2.2e \n',t,max(N),min(N))
end



%% Auxillary math functions
function y = dirac_delta(i,j)
if i == j
    y = 1;
else
    y = 0;
end
end

function y = heaviSide(x)
if x > 0
    y = 1;
elseif x == 0
    y = 1/2;
else
    y = 0;
end
end

function y = lambdaplus(i,x,xcenter)
y = (x - xcenter(i+1))/(xcenter(i) - xcenter(i+1));
end

function y = lambdaminus(i,x,xcenter)
y = (x - xcenter(i-1))/(xcenter(i) - xcenter(i-1));
end