µli,j=
bk¯m,l+1c
X
k=b˘ k¯m,l+1c k−1˘
X
b=b˘ b¯m,l+1c
k˘i·b˘j·n(t,˘ k,˘ b).˘ (4.28)
Furthermore, one needs to substitute the Dirac deltas in the approximation of the NDF by Kronecker deltas which requires pivot at integer values. The pivot is almost always taken to be in the middle of the element which would not be at an integer value. This problem is solved by simply rounding. For the monovariate case
k˘l=round(k¯m,l,k¯m,l+1) (4.29)
˘ n(t,k)˘ ≈
N
X
i=1
¯ wi·δ˘˘k
l,k˘ (4.30)
Not having pivots on the boundary linek−1=bcould introduce problems and, accordingly, pivots are forced to be on this lines. This leads to
k˘l=round(k¯m,l,k¯m,l+1) (4.31) b˘l=round(b¯m,l,b¯m,l+1) (4.32)
˘ n(t,k)˘ ≈
N
X
i=1
¯
wi·δ˘˘kl,k˘·δ˘b˘l,b˘. (4.33)
The moments can again be expressed by a linear product of the zeroth order moment of the death term and a time-independent multiplicator.
4.3 Implementation
In this section the details of the numerical implementation are described.
4.3.1 Meshing
The division of the computational domain into elements is required for applying any sectional technique.
A mesh is good if a high accuracy is obtained at a moderate computational cost.
Linear polymer
Here, the meshing strategy proposed by Ho et al. [127] is used. The mesh is visualized in Figure 4.1a.
Ho et al. [127] proposed a mesh divided in three parts: The first part up tok˘c has exactly represented polymers. It is important to know the products (monomer and small oligomers) precisely. Therefore, the equations describing those oligomers are solved directly without using any class method. The second
part starts with a bin size of 1 and each bin increases in size by 1. The third part has a geometric mesh. This means that the width of element l+1is equal to the width of element ltimes the growth factor. In order to not have any jumps in the mesh size increase it is demanded that the first increase of the geometric mesh is the equal to the increase of the linear mesh. The meshing is continued until a specified amount of monomer units is within the meshed area.
Branched polymer
Several authors have studied how to best set up a multi-dimensional grid [224, 237]. However, nobody has studied breakage or depolymerization in great detail. Here a simple meshing strategy is used.
This strategy is visualized in Figure 4.1b The meshing strategy for the linear mesh is used for both coordinates with the same parameters, except that the exactly represented branching bonds start from 0 and go up tok˘c−1and that the maximal amount of branching bonds is set such that almost no polymer has a greater amount of branching bonds than the maximal amount. This results in one completely exact region, one region in which the class method is only used in the amount of monomer unit, and the largest region where the class method is used for both coordinates. This domain has rectangular and triangular elements. After this meshing is done, the elements above a specified amount of branching bonds, which is chosen to be sufficiently large, are cut away.
4.3.2 Computations prior to the solution of the Ordinary Differential Equations
The matrix describing the birth term for the Fixed Pivot technique and the initial condition can be eval-uated prior to solving the ODEs.
Distributing the birth term
For all pivots with less than or equal to k˘c monomer units no distribution is necessary. And for all pivots with less than or equal to k˘c−1branching bonds the equations have already been provided in Equations (4.15) and (4.16). For the remaining pivots a space filling Delaunay triangulation using the MATLAB functiondelaunayis used. Then the MATLAB algorithmtriangulationis used to identify the element in which the birth occurs and also to provide already the distribution to each corner. This algo-rithm, however, returns small negative values if the birth occurs close to the boundary of the elements.
In this case the birth is shifted onto the boundary and the distribution of the birth is performed using the 1D algorithm. For CAT and FP the same algorithm is used. However, using CAT the distribution of the birth term has to computed with every evaluation of the hydrolysis rate.
Initial condition
As the first step in initialization the amount of polymers, monomer units, and branching bonds in each bin is obtained. Then using the algorithm to distribute the birth term, which has been described just
1 5 10 15 20 25 30 35 40 45 50 55 60 Amount of monomer units
Exactly represented
Linear part Geometric part
(a) The mesh for a linear polymer withk˘c=15. In green the linear part of the mesh and in orange the geometric part.
10 20 30 40 50
0 5 10 15
Amount of monomer units
Amount of branching bonds
(b) The mesh for a branched polymer with k˘c=5 and maximal 10 branching bonds. In green the exact representation in the amount of branching bonds and in orange the approximation for both co-ordinates.
Figure 4.1: The growth factor is 1.3. The maximal amount of monomer units is 55. The blue dots represent the pivots and the black bars are separators between the bins.
before this, this is distributed to the bins. This guarantees that the amount of polymers, monomer units, and branching bonds is computed to values close to the machine precision. By just assigning the amount of polymer units to each pivot, the amount of monomer units would not be obtained with a high accuracy with a coarse mesh.
For a small amount of monomer units and linear polymers the required moments can be obtained by just evaluating the initial distribution in each point in a reasonable time (<10 minutes). However, for branched polymer the computational time and memory requirement scales quadratically with the amount of monomer units. This would lead to requiring days with large branched polymer. Therefore, for such large systems the initial distribution in monomer directions is approximated as a Gamma distri-bution for more than1×105 monomer units and in branching bond direction with a normal distribution for more than1×105 branching bonds. Then instead of evaluating at every point and summation, this can be replaced by a numerical integration with the MATLAB integratorintegral2.
Normalization
The polymer concentration is low which would cause numerical problems. Furthermore, large polymers have more monomer units than small polymers and should, therefore, be weighted more strongly in decisions about the step size. For FP this can be implemented easily offline. However, for CAT parts of the computations would have to be performed online, because of the non-linear nature of the distribution of the birth term. Therefore, only normalization with the polymer concentration is performed
Normalization for the Fixed Pivot technique In order to satisfy this requirements, the initial amount at each pivot is normalized by multiplication with the amount of monomer units of the pivot and division
by the total initial concentration
˜
wl= wl
PN
i=1wi(t=0)·kl. (4.34)
Accordingly, the inverse of the Michaelis constants needs to be normalized as well I(k˜ l,bl) =I˜l= Il
Furthermore, the matrices relating the zeroth order moment of the death term to the other moments need to be normalized
This leads to a normalized distribution matrixP
P˜•=
Normalization for the Cell Average technique The initial amount of polymer at each pivot is nor-malized by division with the total initial number concentration
˜
wl= wl PN
i=1wi(t=0). (4.38)
Only the inverse of the Michealis constants needs to be normalized as well I(k˜ l,bl) =I˜l = Il
4.3.3 Solving the system of Ordinary Differential Equations
The resulting system of ODEs is solved using the MATLAB solver ode45with a relative and absolute tolerance of1×10−6demanding that the polymer concentration is positive. If due to rounding errors a negative value occurs, it is set to zero.