How to code a certain part of this equation?

2 Ansichten (letzte 30 Tage)
Angshuman Podder am 17 Mai 2022
Bearbeitet: Angshuman Podder am 18 Mai 2022
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
0 Kommentare-2 ältere Kommentare anzeigen-2 ältere Kommentare ausblenden

Melden Sie sich an, um zu kommentieren.

Akzeptierte Antwort

Torsten am 17 Mai 2022
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 !
1 Kommentar-1 ältere Kommentare anzeigen-1 ältere Kommentare ausblenden
Angshuman Podder am 18 Mai 2022
Bearbeitet: Angshuman Podder am 18 Mai 2022
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

Melden Sie sich an, um zu kommentieren.

Kategorien

Mehr zu Special Functions finden Sie in Help Center und File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by