Loading [MathJax]/jax/output/SVG/jax.js
Research article Special Issues

High resolution finite difference schemes for a size structured coagulation-fragmentation model in the space of radon measures

  • In this paper, we develop explicit and semi-implicit second-order high-resolution finite difference schemes for a structured coagulation-fragmentation model formulated on the space of Radon measures. We prove the convergence of each of the two schemes to the unique weak solution of the model. We perform numerical simulations to demonstrate that the second order accuracy in the Bounded-Lipschitz norm is achieved by both schemes.

    Citation: Azmy S. Ackleh, Rainey Lyons, Nicolas Saintier. High resolution finite difference schemes for a size structured coagulation-fragmentation model in the space of radon measures[J]. Mathematical Biosciences and Engineering, 2023, 20(7): 11805-11820. doi: 10.3934/mbe.2023525

    Related Papers:

    [1] Azmy S. Ackleh, Rainey Lyons, Nicolas Saintier . Finite difference schemes for a structured population model in the space of measures. Mathematical Biosciences and Engineering, 2020, 17(1): 747-775. doi: 10.3934/mbe.2020039
    [2] Azmy S. Ackleh, Vinodh K. Chellamuthu, Kazufumi Ito . Finite difference approximations for measure-valued solutions of a hierarchicallysize-structured population model. Mathematical Biosciences and Engineering, 2015, 12(2): 233-258. doi: 10.3934/mbe.2015.12.233
    [3] Azmy S. Ackleh, Mark L. Delcambre, Karyn L. Sutton, Don G. Ennis . A structured model for the spread of Mycobacterium marinum: Foundations for a numerical approximation scheme. Mathematical Biosciences and Engineering, 2014, 11(4): 679-721. doi: 10.3934/mbe.2014.11.679
    [4] Horst R. Thieme . Discrete-time population dynamics on the state space of measures. Mathematical Biosciences and Engineering, 2020, 17(2): 1168-1217. doi: 10.3934/mbe.2020061
    [5] Qihua Huang, Hao Wang . A toxin-mediated size-structured population model: Finite difference approximation and well-posedness. Mathematical Biosciences and Engineering, 2016, 13(4): 697-722. doi: 10.3934/mbe.2016015
    [6] Inom Mirzaev, David M. Bortz . A numerical framework for computing steady states of structured population models and their stability. Mathematical Biosciences and Engineering, 2017, 14(4): 933-952. doi: 10.3934/mbe.2017049
    [7] Suzanne M. O'Regan, John M. Drake . Finite mixture models of superspreading in epidemics. Mathematical Biosciences and Engineering, 2025, 22(5): 1081-1108. doi: 10.3934/mbe.2025039
    [8] Ugur G. Abdulla, Vladislav Bukshtynov, Saleheh Seif . Cancer detection through Electrical Impedance Tomography and optimal control theory: theoretical and computational analysis. Mathematical Biosciences and Engineering, 2021, 18(4): 4834-4859. doi: 10.3934/mbe.2021246
    [9] József Z. Farkas, Peter Hinow . Physiologically structured populations with diffusion and dynamic boundary conditions. Mathematical Biosciences and Engineering, 2011, 8(2): 503-513. doi: 10.3934/mbe.2011.8.503
    [10] Nalin Fonseka, Jerome Goddard Ⅱ, Alketa Henderson, Dustin Nichols, Ratnasingham Shivaji . Modeling effects of matrix heterogeneity on population persistence at the patch-level. Mathematical Biosciences and Engineering, 2022, 19(12): 13675-13709. doi: 10.3934/mbe.2022638
  • In this paper, we develop explicit and semi-implicit second-order high-resolution finite difference schemes for a structured coagulation-fragmentation model formulated on the space of Radon measures. We prove the convergence of each of the two schemes to the unique weak solution of the model. We perform numerical simulations to demonstrate that the second order accuracy in the Bounded-Lipschitz norm is achieved by both schemes.



    Coagulation-fragmentation (CF) equations have been used to model many physical and biological phenomena [1,2]. In particular, when combined with transport terms, these equations can be used to model the population dynamics of oceanic phytoplankton [3,4,5]. Setting such models in the space of Radon measures allows for the unified study of both discrete and continuous structures. Not only are the classical discrete and continuous CF equations special cases of the measure valued model (as shown in [6]), but this setting allows for a mixing of the two structures, which has become of interest in particular applications [7,8].

    With the above applications in mind, numerical schemes to solve CF equations are of great importance to researchers. In particular, finite difference methods offer numerical schemes which are easy to implement and approximate the solution with a high order of accuracy. The latter benefit is especially important in the study of stability and optimal control of such equations.

    The purpose of this article is to make improvements on two of the three first-order schemes presented in [9], namely the fully explicit and semi-implicit schemes. These schemes are shown to have certain advantages and disadvantages as discussed in the aforementioned study. In particular, the fully explicit scheme has the qualitative property of conservation of mass through coagulation. On the other hand, the semi-implicit scheme has a more relaxed Courant–Friedrichs–Lewy (CFL) condition, which does not depend on the initial condition. We have decided not to attempt to improve the third scheme presented in [9] as there does not seem to be a significant advantage of the named conservation law scheme to outweigh the drastic computational cost. The main improvement here is to lift-up these two first-order schemes to second-order ones on the space of Radon measures; however, as this state space contains singular elements (including point measures), the improvement of these schemes must be handled with care. As shown in [10], discontinuities and singularities in the solution can cause drastic changes in not only the order of convergence of the scheme, but also in the behavior of the scheme. To address these issues, we turn to a high resolution scheme studied with classical structured population models (i.e. without coagulation-fragmentation) in [11,12,13]. This scheme makes use of a minmod flux limiter to control any oscillatory behavior of the scheme caused by irregularities. With this new flux, we show that it is possible for second order convergence rates to be obtained for continuous density solutions. However, as the solutions become more irregular, one should expect the convergence rate to decline. Such a phenomenon is demonstrated in [10,13], and we direct the reader to these manuscripts for more discussion.

    The layout of the paper is as follows. In Section 2, we present the notation and preliminary results about the model and state space used throughout the paper. In Section 3, we describe the model and state all assumptions imposed on the model parameters. In Section 4, we present the numerical schemes, their CFL conditions, and state the main theorem of the paper. In Section 5, we test the convergence rate of the schemes against well-known examples. In Section 6, we provide a conclusion and in the Appendix (Section 7) we provide proofs for some of our results.

    We make use of standard notations for function spaces. The most common examples of these are C1(R+) for the space of real valued continuously differentable functions and W1,(R+) for the usual Sobelov space. The space of Radon measures will be denoted with M(R+), with M+(R+) representing its positive cone. This space will be equipped with the Bounded-Lipschitz (BL) norm given by

    μBL:=supϕW1,1{R+ϕ(x)μ(dx):ϕW1,(R+)}.

    Another norm of interest to this space is the well studied Total Variation (TV) norm given by

    νTV=|ν|(R+)=supf1{R+f(x)ν(dx):fCc(R+)}.

    For more information about these particular norms and their relationship we direct the reader to [14,15]. For lucidity, we use operator notation in place of integration when we believe it necessary, namely

    (μ,f):=Af(x)μ(dx),

    where the set A is the support of the measure μ. Finally, we denote the minmod function by mm(a,b) and use the following definition

    mm(a,b):=sign(a)+sign(b)2max(|a|,|b|).

    The model of interest is the size-structured coagulation fragmentation model given by

    {tμ+x(g(t,μ)μ)+d(t,μ)μ=K[μ]+F[μ],(t,x)(0,T)×(0,),g(t,μ)(0)Ddxμ(0)=R+β(t,μ)(y)μ(dy),t[0,T],μ(0)=μ0M+(R+),, (3.1)

    where μ(t)M+(R+) represents individuals' size distribution at time t and the functions g,d,β are their growth, death, and reproduction rate, respectively. The coagulation and fragmentation processes of a population distributed according to μM+(R+) are modeled by the measures K[μ] and F[μ] given

    (K[μ],ϕ)=12R+R+κ(y,x)ϕ(x+y)μ(dx)μ(dy)R+R+κ(y,x)ϕ(x)μ(dy)μ(dx)

    and

    (F[μ],ϕ)=R+(b(y,),ϕ)a(y)μ(dy)R+a(y)ϕ(y)μ(dy)

    for any test function ϕ. Here, κ(x,y) is the rate at which individuals of size x coalesce with individuals of size y, a(y) is the global fragmentation rate of individuals of size y, and b(y,) is a measure supported on [0,y] such that b(y,A) represents the probability a particle of size y fragments to a particle with size in the Borel set A.

    Definition 3.1. Given T0, we say a function μC([0,T],M+(R+)) is a weak solution to (3.1) if for all ϕ(C1W1,)([0,T]×R+) and for all t[0,T], the following holds:

    0ϕ(t,x)μt(dx)0ϕ(0,x)μ0(dx)=t00[tϕ(s,x)+g(s,μs)(x)xϕ(s,x)d(s,μs)(x)ϕ(s,x)]μs(dx)ds+t0(K[μs]+F[μs],ϕ(s,))ds+t00ϕ(s,0)β(s,μs)(x)μs(dx)ds. (3.2)

    For the numerical scheme, we will restrict ourselves to a finite domain, [0,xmax]. Thus, we impose the following assumptions on the growth, death and birth functions:

    (A1) For any R>0, there exists LR>0 such that for all μiTVR and ti[0,) (i=1,2) the following hold for f=g,d,β,

    f(t1,μ1)f(t2,μ2)LR(|t1t2|+μ1μ2BL),

    (A2) There exists ζ>0 such that for all T>0,

    supt[0,T]supμM+(R+)g(t,μ)W1,+d(t,μ)W1,+β(t,μ)W1,<ζ,

    (A3) For all (t,μ)[0,)×M+(R+),

    g(t,μ)(0)>0andg(t,μ)(xmax)=0

    for some large xmax>0.

    We assume that the coagulation kernel κ satisfies the following assumption:

    (K1) κ is symmetric, nonnegative, bounded by a constant Cκ, and globally Lipschitz with Lipschitz constant Lκ.

    (K2) \(\kappa(x, y) = 0 \text{ whenever } x+y > x_{\max}.\)

    We assume that the fragmentation kernel satisfies the following assumptions:

    (F1) aW1,(R+) is non-negative,

    (F2) for any y0, b(y,dx) is a measure such that

    (i) b(y,dx) is non-negative and supported in [0,y], and there exist a Cb>0 such that b(y,R+)Cb for all y>0,

    (ii) there exists Lb such that for any y,ˉy0,

    b(y,)b(ˉy,)BLLb|yˉy|

    (iii) for any y0,

    (b(y,dx),x)=y0xb(y,dx)=y.

    The existence and uniqueness of mass conserving solutions of model (3.1) under these assumptions were established in [6].

    We adopt the numerical discretization presented in [6]. For some fixed mesh sizes Δx,Δt>0, we discretize the size domain [0,xmax] with the cells

    ΛΔxj:=[(j12)Δx,(j+12)Δx), for j=1,,J,

    and

    ΛΔx0:=[0,Δx2).

    We denote the midpoints of these grids by xj. The initial condition μ0M+(R+) will be approximated by a combination of Dirac measures

    μΔx0=Jj=0m0jδxj, where m0j:=μ0(ΛΔxj).

    We first approximate the model coefficients κ, a, b as follows. For the physical ingredients, we define

    aΔxi=1ΔxΛΔxia(y)dy,κΔxi,j=1Δx2ΛΔxi×ΛΔxjκ(x,y)dxdy

    for i,j1, and

    aΔx0=2ΔxΛΔx0a(y)dy,κΔx0,0=4Δx2ΛΔx0×ΛΔx0κ(x,y)dxdy

    (with the natural modifications for κΔx0,j and κΔxi,0, i1). We then let aΔxW1,(R+) and κΔxW1,(R+×R+) be the linear interpolation of the aΔxi and κΔxi,j, respectively. Finally, we define the measure bΔx(xj,)M+(ΔxN) by

    bΔx(xj,)=ijb(xj,ΛΔxi)δxj=:ijbΔxj,iδxj

    and then bΔx(x,)M+(ΔxN0) for x0 as the linear interpolate between the bΔx(xj,). When the context is clear, we omit the Δx from the notation above.

    We make use of these approximations to combine the high-resolution scheme presented in [13] with the fully explicit and semi-implicit schemes presented in [9]. Together, these schemes give us the numerical scheme

    {mk+1j=mkjΔtΔx(fkj+12fkj12)Δtdkjmkj+Δt(Cj,k+Fj,k),j=1,..,J,gk0mk0=ΔxJj=1βkjmkj:=Δx(32βk1mk1+12βkJmkJ+J1j=2βkjmkj), (4.1)

    where the flux term is given by

    fkj+12={gkjmkj+12(gkj+1gkj)mkj+12gkjmm(Δ+mkj,Δmkj)j=2,3,,J2gkjmkjj=0,1,J1,J, (4.2)

    the fragmentation term, Fj,k, is given by

    Fj,k:=Ji=jbi,jaimkiajmkj, (4.3)

    and the coagulation term, Cj, is either given by an explicit discretization as

    Cexpj,k:=12j1i=1κi,jimkimkjiJi=1κi,jmkimkj, (4.4)

    or by an implicit one as

    Cimpj,k:=12j1i=1κi,jimk+1imkjiJi=1κi,jmkimk+1j. (4.5)

    As discussed in [9], the explicit scheme which uses (4.4) to approximate the coagulation term and the semi-implicit scheme which instead uses (4.5) to approximate the coagulation term behave differently with respect to the mass conservation and have different Courant–Friedrichs–Lewy (CFL) conditions. The assumed CFL condition for the schemes are

    Explicit:Δt(Cκμ0TVexp((ζ+CbCa)T)+Camax{1,Cb}+(1+32Δx)ζ)1Semi-Implicit:ˉζ(2+32Δx)Δt1, (4.6)

    where ˉζ=max{ζ,aW1,}, Ca=a. The CFL conditions above are similar to those used in [9], but are adjusted due to the flux limiter term as in [13]. It is clear that the semi-implicit scheme has a less restrictive and simpler CFL condition than the explicit scheme. In particular, the CFL condition of the semi-implicit scheme is independent on the initial condition, unlike its counterpart. The trade off for this is a loss of qualitative behavior of the scheme in the sense of mass conservation. Indeed as shown in [9], when β=d=g=a=0, the semi-implicit coagulation term does not conserve mass, whereas the explicit term does. As shown in the appendix, this loss is controlled by the time step size, Δt.

    It is useful to define the following coefficients:

    Akj={gkjj=1,J,12(gkj+1+gkj+gkjmm(Δ+mkj,Δmkj)Δmkj)j=2,12(gkj+1+gkj+gkjmm(Δ+mkj,Δmkj)Δmkjgkj1mm(Δmkj,Δmkj1)Δmkj)j=3,,J2,12(2gkjgkj1mm(Δmkj,Δmkj1)Δmkj)j=J1,

    and

    Bkj={Δgkjj=1,J,12Δ+gkjj=2,12(Δ+gkj+Δgkj)j=3,,J2,12Δgkjj=J1..

    Notice, |Akj|3Δt2Δxζ and AkjBkj0 as

    2(AkjBkj)={2gkj1j=1,J,gkj(2+mm(Δ+mkj,Δmkj)Δmkj)j=2,gkj(1+mm(Δ+mkj,Δmkj)Δmkj)+gkj1(1mm(Δmkj,Δmkj1)Δmkj)j=3,,J2,gnj+gnj1(1mm(Δmnj,Δmnj1)Δmnj)j=J1..

    Scheme (4.1) can then be rewritten as

    {mk+1j=(1ΔtΔxAkjΔt(dkj+aj))mkj+ΔtΔx(AkjBkj)mkj1+ΔtJi=jbi,jaimki+Δt.Cj,kgk0mk0=ΔxJj=1βkjmkj.. (4.7)

    Depending on the choice of coagulation term, this formulation leads to either

    {mk+1j=(1ΔtΔxAkjΔt(dkj+aj)ΔtJi=1κi,jmki)mkj+ΔtΔx(AkjBkj)mkj1+ΔtJi=jbi,jaimki+Δt2j1i=1κi,jimkimkjigk0mk0=ΔxJj=1βkjmkj, (4.8)

    for the explicit term, Cexpj,k, or

    {(1+ΔtJi=1κi,jmki)mk+1j=(1ΔtΔxAkjΔt(dkj+aj))mkj+ΔtΔx(AkjBkj)mkj1+ΔtJi=jbi,jaimki+Δt2j1i=1κi,jimk+1imkjigk0mk0=ΔxJj=1βkjmkj,. (4.9)

    for the implicit term, Cimpj,k.

    For these, schemes, we have the following Lemmas which are proven in the appendix:

    Lemma 4.1. For each k=1,2,,ˉk,

    mkj0 for all j=1,2,J,

    μkΔxTVμ0TVexp((ζ+CbCa)T).

    Lemma 4.2. For any l,p=1,2,,ˉk,

    μlΔxμpΔxBLLT|lp|.

    Using the above two Lemmas, we can arrive at analogous results for the linear interpolation (4.10):

    μΔtΔx(t):=μ0Δxχ{0}(t)+ˉk1k=0[(1tkΔtΔt)μkΔx+tkΔtΔtμk+1Δx]χ(kΔt,(k+1)Δt](t). (4.10)

    Thus by the well know Ascoli-Arzela Theorem, we have the existence of a convergent subsequence of the net {μΔtΔx(t)} in C([0,T],M+([0,xmax]). We now need only show any convergent subsequence converges to the unique solution (3.2).

    Theorem 4.1. As Δx,Δt0 the sequence μΔtΔx converges in C([0,T],M+([0,xmax])) to the solution of (3.1).

    Proof. By multiplying (4.1) by a superfluously smooth test function ϕ(W1,C2)([0,T]×R), denoting ϕkj:=ϕ(kΔt,xj), summing over all j and k, and rearranging we arrive at

    ˉk1k=0Jj=1((mk+1jmkj)ϕkj+ΔtΔx(fkj+12fkj12)ϕkj)+Δtˉk1k=0j=1dkjmkjϕkj=Δtˉk1k=1Jj=1ϕkj(12j1i=1κi,jimkimkjiJi=1κi,jmkimkj+Ji=jbi,jaimkiajmkj). (4.11)

    The left-hand side of equation (4.11) was shown in [13] to be equivalent to

    xmax0ϕ(T,x)dμˉkΔx(x)xmax0ϕ(0,x)dμ0Δx(x)Δtˉk1k=0(xmax0tϕ(tk,x)dμkΔx(x)+xmax0xϕ(tk,x)g(tk,μkΔx)(x)dμkΔx(x)R+d(tk,μkΔx)(x)ϕ(tk,x)dμkΔx(x)+xmax0ϕ(tk,Δx)β(tk,μkΔx)(x)dμkΔx(x))+o(1),

    where o(1)0 as Δt,Δx0.

    The right-hand side of (4.11) was shown in [9] to be equal to

    Δtˉk1k=1{(K[μΔtΔx(tk)],ϕ(tk,))+(F[μΔtΔx(tk)],ϕ(tk,))}+O(Δx).

    Making use of results, it is then easy to see (4.11) is equivalent to

    xmax0ϕ(T,x)dμΔtΔx(T)(x)xmax0ϕ(0,x)dμ0Δx(x)=T0(xmax0tϕ(t,x)+xϕ(t,x)g(t,μΔtΔx(t))(x)dμΔtΔx(t)(x)xmax0d(t,μΔtΔx(t))(x)ϕ(t,x)dμΔtΔx(t)(x)+xmax0ϕ(t,Δx)β(t,μΔtΔx(t))(x)dμΔtΔx(t)(x))dt+T0(K[μΔtΔx(t)],ϕ(t,))+(F[μΔtΔx(t)],ϕ(t,))dt+o(1).

    Passing the limit as Δt,Δx0 along a converging subsequence, we then obtain that equation (3.2) holds for any ϕ(C2W1,)([0,T]×R+) with compact support. A standard density argument shows that equation (3.2) holds for any ϕ(C1W1,)([0,T]×R+). As the weak solution is unique [6], we conclude the net {μΔtΔx} converges to the solution of model (3.1).

    We point out that while these schemes are higher-order in space, they are only first order in time. To lift these schemes into a second-order in time as well, we make use of the second-order Runge-Kutta time discretization [16] for the explicit scheme and second-order Richardson extrapolation [17] for the semi-implicit scheme.

    In this section, we provide numerical simulations which test the order of the explicit and semi-implicit schemes developed in the previous sections. We test each component separately, beginning first with a pure coagulation equation in example 1 (where we set g=β=d=a=0), then a pure fragmentation equation in example 2 (where we set g=β=d=κ=0). In example 3, we consider all components of model (3.1) including the boundary term which is implemented as in scheme (4.7). For readers interested in the schemes performance in the absence of the coagulation-fragmentation processes, we direct you to [11,12,13]. For each example, we give the BL error and the order of convergence. To appreciate the gain in the order of convergence compared to those studied in [9], which are based on a first order approximation of the transport term, we add some of the numerical results from the scheme presented in [9].

    In some of the following examples, the exact solution of the model problem is given. In these cases, we approximate the order of accuracy, q, with the standard calculation:

    q=log2(ρ(μΔtΔx(T),μ(T))ρ(μ0.5Δt0.5Δx(T),μ(T)))

    where μ represents the exact solution of the examples considered. In the cases where the exact solutions are unknown, we approximate the order by

    q=log2(ρ(μΔtΔx(T),μ2Δt2Δx(T))ρ(μ0.5Δt0.5Δx(T),μΔtΔx(T)))

    and we report the numerator of the log argument as the error. The metric ρ we use here was introduced in [18] and is equivalent to the BL metric, namely

    Cρ(μ,ν)μνBLρ(μ,ν)

    for some constant C (dependent on the finite domain). As discussed in [18], this metric is more efficient to compute than the BL norm and maintains the same order of convergence. An alternative to this algorithm would be to make use of the algorithms presented in [19], where convergence in the Fortet-Mourier distance is considered.

    Example 1 In this example, we test the quality of the finite difference schemes against coagulation equations. To this end, we take κ(x,y)1 and μ0=exdx with all other ingredients set to 0. This example has an exact solution given by

    μt=(22+t)2exp(22+tx)dx

    see [20] for more details. The simulation is performed over the truncated domain (x,t)[0,20]×[0,0.5]. We present the BL error and the numerical order of convergence for both schemes in Table 1.

    Table 1.  Error, order, and computation time for example 1. Here, Nx and Nt represent the number of points in x and t, respectively. The numerical result in the last row for the 1st order variant is generated from the scheme presented in [9].
    Explicit Semi-Implicit
    Nx Nt BL Error Order Time (secs) BL Error Order Time (secs)
    100 250 0.0020733 1.0374 0.0020886 0.79633
    200 500 0.00054068 1.9391 6.8224 0.00054408 1.9407 5.1724
    400 1000 0.00013802 1.9699 98.525 0.00013883 1.9705 73.298
    800 2000 3.4842e-05 1.9860 2430.2 3.5040e-05 1.9862 1792.5
    1600 4000 8.7417e-06 1.9948 43381 8.7906e-06 1.9950 32361
    Explicit (1st order) Semi-Implicit (1st order)
    800 2000 0.015675 0.96974 523.11 0.010996 0.97418 1393.3

     | Show Table
    DownLoad: CSV

    Example 2 In this example, we test the quality of the finite difference scheme against fragmentation equations. We point out that in this case, the two schemes are identical in the spacial component. For this demonstration, we take μ0=exdx, b(y,)=2ydx and a(x)=x. As given in [21], this problem has an exact solution of

    μt=(1+t)2exp(x(1+t))dx.

    The simulation is performed over the finite domain (x,t)[0,20]×[0,0.5]. We present the BL error and the numerical order of convergence for both schemes in Table 2. Note as compared to coagulation, the fragmentation process is more affected by the truncation of the domain. This results in the numerical order of the scheme being further from 2 than example 1.

    Table 2.  Error, order, and computation time for example 2. Here, Nx and Nt represent the number of points in x and t, respectively. The numerical result in the last row for the 1st order variant is generated from the scheme presented in [9].
    Explicit Semi-Implicit
    Nx Nt BL Error Order Time (secs) BL Error Order Time (secs)
    100 250 0.0053857 1.0148 0.0053836 0.78499
    200 500 0.0014548 1.8883 6.7398 0.0014536 1.8890 5.1448
    400 1000 0.00037786 1.9449 99.38 0.00037753 1.9449 73.587
    800 2000 9.6317e-05 1.9720 2369.4 9.6322e-05 1.9707 1763.3
    1600 4000 2.4468e-05 1.9769 43512 2.4514e-05 1.9743 32585
    Explicit (1st order) Semi-Implicit (1st order)
    800 2000 0.059804 0.9128 574.91 0.096943 0.86667 1368.9

     | Show Table
    DownLoad: CSV

    Example 3 In this example, we test the schemes against the complete model (i.e., with all biological and physical processes). To this end, we take μ0=exdx, g(x)=22ex20, β(x)=2, d(x)=1, κ(x,y)=1, a(x)=x, and b(y,)=2y. The simulation is performed over the finite domain (x,t)[0,20]×[0,0.5]. To our knowledge, the solution of this problem is unknown.

    Table 3.  Error, order, and computation time for example 3. Here, Nx and Nt represent the number of points in x and t, respectively. The numerical result in the last row for the 1st order variant is generated from the scheme presented in [9].
    Explicit Semi-Implicit
    Nx Nt BL Error Order Time (secs) BL Error Order Time (secs)
    100 250 0.0023026 1.0332 0.0028799 0.74398
    200 500 0.00085562 1.4282 6.8831 0.00076654 1.9096 5.5104
    400 1000 0.0002743 1.6412 100.57 0.00076654 1.9549 75.055
    800 2000 7.5404e-05 1.8631 2371.2 5.021e-05 1.9775 1739.1
    1600 4000 1.9495e-05 1.9515 43779 1.2651e-05 1.9887 32286
    Explicit (1st order) Semi-Implicit (1st order)
    800 2000 0.0092432 0.97728 625.3 0.0014192 0.98355 1112.8

     | Show Table
    DownLoad: CSV

    Example 4 As mentioned in [9], the mixed discrete and continuous fragmentation model studied in [7,8], with adjusted assumptions, is a special case of model (3.1). Indeed, by removing the biological and coagulation terms and letting the kernel

    (b(y,),ϕ)=Ni=1bi(y)ϕ(ih)+yNhϕ(x)bc(y,x)dx

    with supp bc(y,)[Nh,y] for some h>0, we have the mixed model in question. We wish to demonstrate the finite difference scheme presented here maintains this mixed structure.

    To this end, we take the fragmentation kernel

    bc(y,x)=2y,bi(y)=2y, and a(x)=x1,

    with initial condition μ=5i=1δi+χ[5,15](x)dx, where χA represents the characteristic function over the set A. This is similar to some examples in [8], where more detail and analysis are provided. In Figure 1, we present the simulation of this example. Notice, the mixed structure is preserved in finite time. For examples of this type, the scheme could be improved upon by the inclusion of mass conservative fragmentation terms similar to those presented in [6].

    Figure 1.  Initial condition and numerical solution at time T=4 of example 4.

    In this paper, we have lifted two of the first order finite difference schemes presented in [9] to second order high resolution schemes using flux limiter methods. The difference between both schemes is only found in the coagulation term, where the semi-implicit scheme is made linear. In context of standard structured population models (i.e. without coagulation or fragmentation), these type of schemes have been shown to be well-behaved in the presences of discontinuities and singularities. This quality makes them a well fit tool for studying PDEs in spaces of measures. We prove the convergence of both schemes under the assumption of natural CFL conditions. The order of convergence of both schemes is then tested numerically with previously used examples.

    In summary, the schemes preform as expected in the presence of smooth initial conditions. In all such simulations, the numerical schemes presented demonstrate a convergence rate of order 2. For simulations with biological terms, this convergence rate is expected to drop when singularities and discontinuities occur, as demonstrated in [13]. Mass conservation of the schemes, an important property for coagulation/fragmentation processes, is discussed in detail in [6,9].

    The research of ASA is supported in part by funds from R.P. Authement Eminent Scholar and Endowed Chair in Computational Mathematics at the University of Louisiana at Lafayette. RL is grateful for the support of the Carl Tryggers Stiftelse via the grant CTS 21:1656.

    In this section, we present the proofs of Lemmas 4.1 and 4.2 for the explicit coagulation term. The semi-implicit term follows from similar arguments in the same fashion as [9].

    Proof of Lemma 4.1

    Proof. We first prove via induction that for any k=1,2,,ˉk, μkΔx satisfies the following:

    (i) μkΔxM+(R+) i.e. mkj0 for all j=1,,J,

    (ii) μkΔxTVμ0ΔxTV(1+(ζ+CbCa)Δt)k.

    Then, the TV bound in the Lemma follows from standard arguments (see e.g. Lemma 4.1 in [9]). We prove this Theorem for the choice of the explicit coagulation term, Cexpj,k, as the implicit case is similar and more straight forward.

    We begin by showing that mk+1j0 for every j=1,2,,J. Notice by way of (4.8), this reduces down to showing

    ΔtΔxAkj+Δt(dkj+aj)+ΔtJi=1κi,jmki1.

    Indeed, by the CFL condition (4.6), induction hypothesis, and

    Ji=1κi,jmkiCκJi=1mki=CκμkΔxTVCκμ0ΔxTVexp((ζ+CbCa)T),

    we arrive at the result.

    For the TV bound, we have since the mkj are non-negative, μkΔxTV=Jj=1mkj. By rearranging (4.8) and summing over j=1,2,,J we have

    μk+1ΔxTVJj=1mkj+ΔtΔxJj=1(fkj12fkj+12)+ΔtJj=1Ji=jbi,jaimki+Δt(12Jj=1j1i=1κi,jimkimkjiJj=1Ji=1κi,jmkimkj). (7.1)

    To bound the right-hand side of equation (7.1), we directly follow the arguments of Lemma 4.1 in [9] which yields

    μk+1ΔxTV(1+(ζ+CaCb)Δt)Jj=1mkj=(1+(ζ+CaCb)Δt)μkΔxTV.

    Using the induction hypothesis, we obtain μk+1ΔxTVμ0ΔxTV(1+(ζ+CbCa)Δt)k+1 as desired.

    Proof of Lemma 4.2

    Proof. For ϕW1,(R+) with ϕW1,1, and denoting ϕj:=ϕ(xj), we have for any k,

    (μk+1ΔxμkΔx,ϕ)=Jj=1(mk+1jmkj)ϕjΔtJj=1ϕj(1Δx(fkj12fkj+12)dkjmkjajmkj+12j1i=1κi,jimkimkjiJi=1κi,jmkimkj+Ji=jbi,jaimki).

    Let C be the right-hand side of the TV-bound from Lemma 4.1, we then see

    (μk+1ΔxμkΔx,ϕ)ΔtΔxJj=1ϕj(fkj12fkj+12)+Δt(ζ+Ca+CbCa+32CκC)C.

    Moreover, since gkJ=0 the sum in the right-hand side takes the form

    ϕ1gk0mk0+J1j=1(ϕj+1ϕj)fkj+12=Δxϕ1Jj=1βkjmkj+J1j=1(ϕj+1ϕj)fkj+123.5ΔxζC.

    We thus obtain

    (μk+1ΔxμkΔx,ϕ)LΔt,L:=(3.5ζ+Ca+CbCa+32CκC)C.

    Taking the supremum over ϕ gives μk+1ΔxμkΔxBLLΔt for any k. The result follows.

    In this section, we consider the semi-implicit scheme (4.9) without any biological ingredients or fragmentation (i.e., a,g,d,β=0) where mass is not conserved (observed by the numerical experiments in Section 5) and show this change is controlled by the time step. For a bound on the loss of mass via fragmentation, we direct the reader to Section 6.1 of [6]. Multiplying (4.9) by xj and summing over j, we can arrive at

    Jj=1xjmk+1j=Jj=1xjmkj+Δt2Jj=1j1i=1xjκi,jimk+1imkjiΔtJj=1Ji=1xjκi,jmkimk+1j.

    In [6, Section 6.1] it was shown that the explicit scheme conserves mass through the coagulation process, i.e.,

    Δt2Jj=1j1i=1xjκi,jimkimkjiΔtJj=1Ji=1xjκi,jmkimkj=0.

    Adding this to the previous equation we have

    Jj=1xjmk+1j=Jj=1xjmkj+Δt2Jj=1j1i=1xjκi,ji(mk+1imki)mkjiΔtJj=1Ji=1xjκi,jmki(mk+1jmkj)=Jj=1xjmkj+Δt2Ji=1Jl=1xl+iκi,l(mk+1imki)mklΔtJj=1Ji=1xjκi,jmki(mk+1jmkj),

    where in the last equality, we change the order of integration and introduce the new index l=ji. Noticing that due to the uniform mesh size xl+i=xl+xi, we can split the second term on the right-hand side and obtain the equation

    Jj=1xjmk+1j=Jj=1xjmkj+Δt2Ji=1Jl=1xiκi,l(mkl(mk+1imki)mki(mk+1lmkl)).

    Since xjxmax, we can bound the last term on the right-hand side

    |Δt2Ji=1Jl=1xiκi,l(mkl(mk+1imki)mki(mk+1lmkl))|ΔtCκxmaxμkΔxTVμk+1ΔxμkΔxBL.

    Using Lemmas 4.1 and 4.2, we have the estimate

    |Jj=1xjmk+1jxjmkjΔt|CκxmaxLexp((ζ+CaCb)T)μ0TVΔt.


    [1] A. B. Burd, G. A. Jackson, Particle aggregation, Ann. Rev. Marine Sci., 1 (2009), 65–90. https://doi.org/10.1146/annurev.marine.010908.163904 doi: 10.1146/annurev.marine.010908.163904
    [2] D. J. Aldous, Deterministic and stochastic models for coalescence (aggregation and coagulation): A review of the Mean-field theory for probabilists, Bernoulli, (1999), 3–48. https://doi.org/10.2307/3318611 doi: 10.2307/3318611
    [3] A. S. Ackleh, B. G. Fitzpatrick, Modeling aggregation and growth processes in an algal population model: Analysis and computations, J. Math. Biol., 35 (1997), 480–502. https://doi.org/10.1007/s002850050062 doi: 10.1007/s002850050062
    [4] A. S. Ackleh, Parameter estimation in a structured algal Coagulation-fragmentation model, Nonlinear Anal. Theory Methods Appl., 28 (1997), 837–854. https://doi.org/10.1016/0362-546X(95)00195-2 doi: 10.1016/0362-546X(95)00195-2
    [5] R. Rudnicki, R. Wieczorek, Fragmentation-Coagulation Models of Phytoplankton, Bulletin Polish Acad. Sci. Math. 54 (2006), 175–191. https://doi.org/10.4064/ba54-2-9 doi: 10.4064/ba54-2-9
    [6] A. S. Ackleh, R. Lyons, N. Saintier, Structured Coagulation-Fragmentation Equation in the Space of Radon Measures: Unifying Discrete and Continuous Models, ESAIM Math. Model. Numer. Anal., 55 (2021). https://doi.org/10.1051/m2an/2021061 doi: 10.1051/m2an/2021061
    [7] G. Baird, E. Süli, A mixed discrete-continuous fragmentation model, J. Math. Anal. Appl., 473 (2019), 273–296 https://doi.org/10.1016/j.jmaa.2018.12.048 doi: 10.1016/j.jmaa.2018.12.048
    [8] G. Baird, E. Süli, A finite volume scheme for the solution of a mixed discrete-continuous fragmentation model, ESAIM Math. Model. Numer. Anal., 55 (2021), 1067–1101. https://doi.org/10.1051/m2an/2020088 doi: 10.1051/m2an/2020088
    [9] A. S. Ackleh, R. Lyons, N. Saintier, Finite difference schemes for a size structured coagulation-fragmentation model in the space of Radon measures, IMA J. Numer. Anal., (2022). https://doi.org/10.1093/imanum/drac071 doi: 10.1093/imanum/drac071
    [10] R. LeVeque, Numerical Mehtods for Conservation Laws, Springer Basel AG, 1992. https://doi.org/10.1007/978-3-0348-8629-1
    [11] J. Shen, C. W. Shu, M. Zhang, High resolution schemes for a hierarchical size-structured model, SIAM J. Numer. Anal., 45 (2007), 352–370. https://doi.org/10.1137/050638126 doi: 10.1137/050638126
    [12] A. S. Ackleh, V. K. Chellamuthu, K. Ito, Finite difference approximations for measure-valued solutions of a hierarchically size-structured population model, Math. Biosci. Eng., 12 (2015), 233–258. https://doi.org/10.3934/mbe.2015.12.233 doi: 10.3934/mbe.2015.12.233
    [13] A. S. Ackleh, R. Lyons, N. Saintier, Finite Difference Schemes for a Structured Population Model in the Space of Measures, Math. Biosci. Eng., 17 (2020), 747–775. https://doi.org/10.3934/mbe.2020039v doi: 10.3934/mbe.2020039v
    [14] C. Düll, P. Gwiazda, A. Marciniak-Czochra, J. Skrzeczkowski, Spaces of measures and their applications to structured population models, Cambridge University Press, 36 (2021). https://doi.org/10.1017/9781009004770
    [15] P. Gwiazda, A. Marciniak-Czochra, H. R. Thieme, Measures Under the Flat Norm as Ordered Normed Vector Space, Positivity, 22 (2017), 105–138. https://doi.org/10.1007/s11117-017-0503-z doi: 10.1007/s11117-017-0503-z
    [16] C. W. Shu, S. Osher, Efficient implementation of essentially non-oscillatory shock capturing schemes, J. Comput. Phys., 77 (1988), 439–471. https://doi.org/10.1016/0021-9991(88)90177-5 doi: 10.1016/0021-9991(88)90177-5
    [17] L. F. Richardson, The Approximate Arithmetical Solution by Finite Differences with an Application to Stresses in Masonry Dams, Philosoph. Transact. Royal Soc. Am., 210 (1911), 307–357. https://doi.org/10.1098/rsta.1911.0009 doi: 10.1098/rsta.1911.0009
    [18] J. Jabłoński, A. Marciniak-Czochra, Efficient Algorithms Computing Distances Between Radon Measures on R, preprint, arXiv: 1304.3501, (2013).
    [19] S. C. Hille, E. S. Theewis, Explicit Expressions and Computational Methods for the Fortet-Mourier Distance to Finite Weighted Sums of Dirac Measures, preprint, arXiv: 2206.12234, (2022).
    [20] D. D. Keck, D. M. Bortz, Numerical Simulation of Solutions and Moments of the Smoluchowski Coagulation Equation, preprint, arXiv: 1312.7240, (2013).
    [21] R. Singh, J. Saha, J. Kumar, A Domain Decomposition Method for Solving Fragmentation and Aggregation Population Balance Equations, J. Appl. Math. Comput., 48 (2015), 265–292. https://doi.org/10.1007/s12190-014-0802-5 doi: 10.1007/s12190-014-0802-5
  • This article has been cited by:

    1. Carlo Bianca, Nicolas Saintier, Thermostatted kinetic theory in measure spaces: Well-posedness, 2025, 251, 0362546X, 113666, 10.1016/j.na.2024.113666
  • Reader Comments
  • © 2023 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(1805) PDF downloads(74) Cited by(1)

Figures and Tables

Figures(1)  /  Tables(3)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog