How to code a certain part of this equation?
1 vue (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Angshuman Podder
le 17 Mai 2022
Modifié(e) : Angshuman Podder
le 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:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1001315/image.png)
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 commentaires
Réponse acceptée
Torsten
le 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 commentaire
Angshuman Podder
le 18 Mai 2022
Modifié(e) : Angshuman Podder
le 18 Mai 2022
Plus de réponses (0)
Voir également
Catégories
En savoir plus sur Special Functions dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!