Home On modeling column crystallizers and a Hermite predictor–corrector scheme for a system of integro-differential equations
Article
Licensed
Unlicensed Requires Authentication

On modeling column crystallizers and a Hermite predictor–corrector scheme for a system of integro-differential equations

  • Khaldoun El Khaldi , Nima Rabiei and Elias G. Saleeby EMAIL logo
Published/Copyright: December 2, 2021

Abstract

Multistaged crystallization systems are used in the production of many chemicals. In this article, employing the population balance framework, we develop a model for a column crystallizer where particle agglomeration is a significant growth mechanism. The main part of the model can be reduced to a system of integrodifferential equations (IDEs) of the Volterra type. To solve this system simultaneously, we examine two numerical schemes that yield a direct method of solution and an implicit Runge–Kutta type method. Our numerical experiments show that the extension of a Hermite predictor–corrector method originally advanced in Khanh (1994) for a single IDE is effective in solving our model. The numerical method is presented for a generalization of the model which can be used to study and simulate a number of possible operating profiles of the column.

Mathematics Subject Classification (2010): 65R20; 45J05; 92E20

Corresponding author: Elias G. Saleeby, Mount Lebanon, Lebanon, E-mail:

Acknowledgments

The authors would like to thank the referee for his thoughtful suggestions to improve the paper.

  1. Author contribution: All the authors have accepted responsibility for the entire content of this submitted manuscript and approved submission.

  2. Research funding: None declared.

  3. Conflict of interest statement: The authors declare no conflicts of interest regarding this article.

Appendix A: Further details of the HPC method

In this appendix, we give the representation of Eqs. (21)(23) in Section 3.3 as linear expressions of the unknowns u t n , u t n , u t n , u t n .

We start by expanding F 0 , F 0 , F 0 , F 1 , F 1 , F 2 (see definition of F 0, F 1, F 2 in Section 3.3):

F 0 s = H t n , s = K t n , s u t n s u s F 0 s = s K t n , s u t n s u s + K t n , s u t n s u s + u t n s u s F 0 s = 2 s 2 K t n , s u t n s u s + 2 s K t n , s u t n s u s + u t n s u s + K t n , s u t n s u s 2 u t n s u s + u t n s u s F 0 s = 3 s 3 K t n , s u t n s u s + 3 2 s 2 K t n , s u t n s u s + u t n s u s + 3 s K t n , s u t n s u s 2 u t n s u s + u t n s u s + K t n , s u t n s u s + 3 u t n s u s 3 u t n s u s + u t n s u s F 1 s = t H ( t , s ) t = t n = t K ( t , s ) t = t n u t n s u s + K ( t n , s ) u t n s u s F 1 s = s t K ( t , s ) t = t n u t n s u s + t K ( t , s ) t = t n u t n s u s + u t n s u s + s K t n , s u t n s u s + K t n , s u t n s u s + u t n s u s F 1 s = 2 s 2 t K ( t , s ) t = t n u t n s u s + 2 s t K ( t , s ) t = t n u t n s u s + u t n s u s + t K ( t , s ) t = t n u t n s u s 2 u t n s u s + u t n s u s + t K ( t , s ) t = t n u t n s u s 2 u t n s u s + u t n s u s + 2 s 2 K t n , s u t n s u s + 2 s K t n , s u t n s u s + u t n s u s + K t n , s u t n s u s 2 u t n s u s + u t n s u s F 2 s = 2 t 2 H ( t , s ) t = t n = 2 t 2 K ( t , s ) t = t n u t n s u s + 2 t K ( t , s ) t = t n u t n s u s + K ( t n , s ) u t n s u s F 2 s = s 2 t 2 K ( t , s ) t = t n u t n s u s + 2 t 2 K ( t , s ) t = t n u t n s u s + u t n s u s + 2 s t K ( t , s ) t = t n u t n s u s

+ t K ( t , s ) t = t n u t n s u s + u t n s u s + s K t n , s u t n s u s + K t n , s u t n s u s + u t n s u s

With constant K(t, s) we obtain:

F 0 s = H t n , s = K u t n s u s F 0 s = K u t n s u s + u t n s u s F 0 s = K u t n s u s 2 u t n s u s + u t n s u s F 0 s = K u t n s u s + 3 u t n s u s 3 u t n s u s + u t n s u s F 1 s = K u t n s u s F 1 s = K u t n s u s + u t n s u s F 1 s = K u t n s u s 2 u t n s u s + u t n s u s F 2 s = K u t n s u s F 2 s = K u t n s u s + u t n s u s

Then we rewrite the expressions of (25)(27) in order to separate the terms F 0, F 1, F 2 and their derivatives, which depend on the unknowns u t n , u t n , u t n , u t n from the summations. Then (25) becomes

(A.1) v ( t ) = 0 t K ( t , s ) u ( t s ) u ( s ) d s 1 i n 1 2 h F 0 t i 1 + F 0 t i + 3 28 h 2 F 0 t i 1 F 0 t i + 1 84 h 3 F 0 t i 1 + F 0 t i + 1 1680 h 4 F 0 t i 1 F 0 t i 1 2 h F 0 t 0 + F 0 t 1 + 3 28 h 2 F 0 t 0 F 0 t 1 + 1 84 h 3 F 0 t 0 + F 0 t 1 + 1 1680 h 4 F 0 t 0 F 0 t 1 + 1 2 h F 0 t n 1 + F 0 t n + 3 28 h 2 F 0 t n 1 F 0 t n + 1 84 h 3 F 0 t n 1 + F 0 t n + 1 1680 h 4 F 0 t n 1 F 0 t n + C 0 ,

where

C 0 = 2 i n 1 1 2 h F 0 t i 1 + F 0 t i + 3 28 h 2 F 0 t i 1 F 0 t i + 1 84 h 3 F 0 t i 1 + F 0 t i + 1 1680 h 4 F 0 t i 1 F 0 t i .

Equation (26) becomes

(A.2) v ( t n ) = H ( t n , t n ) + 0 t n t H ( t n , s ) d s H ( t n , t n ) + 1 i n 1 2 h F 1 t i 1 + F 1 t i + 1 10 h 2 F 1 t i 1 F 1 t i + 1 120 h 3 F 1 t i 1 + F 1 t i , H ( t n , t n ) + 1 2 h F 1 t 0 + F 1 t 1 + 1 10 h 2 F 1 t 0 F 1 t 1 + 1 120 h 3 F 1 t 0 + F 1 t 1 + + 1 2 h F 1 t n 1 + F 1 t n + 1 10 h 2 F 1 t n 1 F 1 t n + 1 120 h 3 F 1 t n 1 + F 1 t n + C 1 ,

where

C 1 = 2 i n 1 1 2 h F 1 t i 1 + F 1 t i + 1 10 h 2 F 1 t i 1 F 1 t i + 1 120 h 3 F 1 t i 1 + F 1 t i .

Equation (27) becomes

(A.3) v ( t n ) = d d t n H ( t n , t n ) + t n H ( t n , s ) s = t n + 0 t n 2 t 2 H ( t n , s ) d s d d t n H ( t n , t n ) + t n H ( t n , s ) s = t n + 1 i n 1 2 h F 2 t i 1 + F 2 t i + 1 12 h 2 F 2 t i 1 F 2 t i , d d t n H ( t n , t n ) + t n H ( t n , s ) s = t n + 1 2 h F 2 t 0 + F 2 t 1 + 1 12 h 2 F 2 t 0 F 2 t 1 + 1 2 h F 2 t n 1 + F 2 t n + 1 12 h 2 F 2 t n 1 F 2 t n + C 3 ,

where C 3 = 2 i n 1 1 2 h F 2 t i 1 + F 2 t i + 1 12 h 2 F 2 t i 1 F 2 t i .

Using the general form of the quadrature rule (24), taking q = 4, let

Q i , ζ = 1 2 h ζ x i 1 + ζ x i + 3 28 h 2 ζ x i 1 ζ x i + 1 84 h 3 ζ x i 1 + ζ x i + 1 1680 h 4 ζ 3 x i 1 ζ 3 x i ;

taking q = 3, let

Q 1 i , ζ = 1 2 h ζ x i 1 + ζ x i + 1 10 h 2 ζ x i 1 ζ x i + 1 120 h 3 ζ x i 1 + ζ x i ;

and finally taking q = 2, let

Q 2 i , ζ = 1 2 h ζ x i 1 + ζ x i + 1 12 h 2 ζ x i 1 ζ x i .

Equation (A.1) becomes

(A.4) v ( t n ) Q 1 , F 0 + Q n , F 0 + C 0 .

Equation (A.2) becomes

(A.5) v ( t n ) H ( t n , t n ) + Q 1 1 , F 1 + Q 1 n , F 1 + C 1 .

Equation (A.3) becomes

(A.6) v ( t n ) d d t H ( t , t ) t = t n + t H ( t , s ) s = t = t n + Q 2 1 , F 2 + Q 2 n , F 2 + C 2 .

We then expand the terms Q, Q 1, Q 2 in (A.4)(A.6). We obtain from (A.4)

Q 1 , F 0 = 1 2 h F 0 t 0 + F 0 t 1 + 3 28 h 2 F 0 t 0 F 0 t 1 + 1 84 h 3 F 0 t 0 + F 0 t 1 + 1 1680 h 4 F 0 t 0 F 0 t 1 = 1 2 h K t n , t 0 u t n t 0 u t 0 + K t n , t 1 u t n t 1 u t 1 + 3 28 h 2 s K t n , t 0 u t n t 0 u t 0 K t n , t 0 u t n t 0 u t 0 + K t n , t 0 u t n t 0 u t 0 s K t n , t 1 u t n t 1 u t 1 + K t n , t 1 u t n t 1 u t 1 K t n , t 1 u t n t 1 u t 1 + 1 84 h 3 2 s 2 K t n , t 0 u t n t 0 u t 0 2 s K t n , t 0 u t n t 0 u t 0 + 2 s K t n , t 0 u t n t 0 u t 0 + K t n , t 0 u t n t 0 u t 0 2 K t n , t 0 u t n t 0 u t 0 + K t n , t 0 u t n t 0 u t 0 + 2 s 2 K t n , t 1 u t n t 1 u t 1 2 s K t n , t 1 u t n t 1 u t 1 + 2 s K t n , t 1 u t n t 1 u t 1 + K t n , t 1 u t n t 1 u t 1 2 K t n , t 1 u t n t 1 u t 1 + K t n , t 1 u t n t 1 u t 1 + 1 1680 h 4 3 s 3 K t n , t 0 u t n t 0 u t 0 3 2 s 2 K t n , t 0 u t n t 0 u t 0 + 3 2 s 2 K t n , t 0 u t n t 0 u t 0 + 3 s K t n , t 0 u t n t 0 u t 0 6 s K t n , t 0 u t n t 0 u t 0 + 3 s K t n , t 0 u t n t 0 u t 0 K t n , t 0 u ( 3 ) t n t 0 u t 0 + 3 K t n , t 0 u t n t 0 u t 0 3 K t n , t 0 u t n t 0 u t 0 + K t n , t 0 u t n t 0 u ( 3 ) t 0 3 s 3 K t n , t 1 u t n t 1 u t 1 + 3 2 s 2 K t n , t 1 u t n t 1 u t 1 3 2 s 2 K t n , t 1 u t n t 1 u t 1 3 s K t n , t 1 u t n t 1 u t 1 + 6 s K t n , t 1 u t n t 1 u t 1 3 s K t n , t 1 u t n t 1 u t 1 + K t n , t 1 u ( 3 ) t n t 1 u t 1 3 K t n , t 1 u t n t 1 u t 1 + 3 K t n , t 1 u t n t 1 u t 1 K t n , t 1 u t n t 1 u ( 3 ) t 1

with constant K t n , s , we obtain:

1 2 h K u t n t 0 u t 0 + K u t n t 1 u t 1 + 3 28 h 2 K u t n t 0 u t 0 + K u t n t 0 u t 0 + K u t n t 1 u t 1 K u t n t 1 u t 1 + 1 84 h 3 K u t n t 0 u t 0 2 K u t n t 0 u t 0 + K u t n t 0 u t 0 + K u t n t 1 u t 1 2 K u t n t 1 u t 1 + K u t n t 1 u t 1 + 1 1680 h 4 K u ( 3 ) t n t 0 u t 0 + 3 K u t n t 0 u t 0 3 K u t n t 0 u t 0 + K u t n t 0 u ( 3 ) t 0 + K u ( 3 ) t n t 1 u t 1 3 K u t n t 1 u t 1 + 3 K u t n t 1 u t 1 K u t n t 1 u ( 3 ) t 1 , Q n , F 0 = 1 2 h F 0 t n 1 + F 0 t n + 3 28 h 2 F 0 t n 1 F 0 t n + 1 84 h 3 F 0 t n 1 + F 0 t n + 1 1680 h 4 F 0 t n 1 F 0 t n = 1 2 h K t n , t n 1 u t n t n 1 u t n 1 + K t n , t n u 0 u t n + 3 28 h 2 s K t n , t n 1 u t n t n 1 u t n 1 K t n , t n 1 u t n t n 1 u t n 1 + K t n , t n 1 u t n t n 1 u t n 1 s K t n , t n u 0 u t n + K t n , t n u 0 u t n K t n , t n u 0 u t n + 1 84 h 3 2 s 2 K t n , t n 1 u t n t n 1 u t n 1 2 s K t n , t n 1 u t n t n 1 u t n 1 + 2 s K t n , t n 1 u t n t n 1 u t n 1 + K t n , t n 1 u t n t n 1 u t n 1 2 K t n , t n 1 u t n t n 1 u t n 1 + K t n , t n 1 u t n t n 1 u t n 1 + 2 s 2 K t n , t n u 0 u t n 2 s K t n , t n u 0 u t n + 2 s K t n , t n u 0 u t n + K t n , t n u 0 u t n 2 K t n , t n u 0 u t n + K t n , t n u 0 u t n + 1 1680 h 4 3 s 3 K t n , t n 1 u t n t n 1 u t n 1 3 2 s 2 K t n , t n 1 u t n t n 1 u t n 1 + 3 2 s 2 K t n , t n 1 u t n t n 1 u t n 1 + 3 s K t n , t n 1 u t n t n 1 u t n 1 6 s K t n , t n 1 u t n t n 1 u t n 1 + 3 s K t n , t n 1 u t n t n 1 u t n 1 K t n , t n 1 u ( 3 ) t n t n 1 u t n 1 + 3 K t n , t n 1 u t n t n 1 u t n 1 3 K t n , t n 1 u t n t n 1 u t n 1 + K t n , t n 1 u t n t n 1 u ( 3 ) t n 1 3 s 3 K t n , t n u 0 u t n + 3 2 s 2 K t n , t n u 0 u t n 3 2 s 2 K t n , t n u 0 u t n 3 s K t n , t n u 0 u t n + 6 s K t n , t n u 0 u t n 3 s K t n , t n u 0 u t n + K t n , t n u ( 3 ) 0 u t n 3 K t n , t n u 0 u t n + 3 K t n , t n u 0 u t n K t n , t n u 0 u ( 3 ) t n

with constant K t n , s , we obtain:

1 2 h K u t n t n 1 u t n 1 + K u 0 u t n + 3 28 h 2 K u t n t n 1 u t n 1 + K u t n t n 1 u t n 1 + K u 0 u t n K u 0 u t n + 1 84 h 3 K u t n t n 1 u t n 1 2 K u t n t n 1 u t n 1 + K u t n t n 1 u t n 1 + K u 0 u t n 2 K u 0 u t n + K u 0 u t n + 1 1680 h 4 K u ( 3 ) t n t n 1 u t n 1 + 3 K u t n t n 1 u t n 1 3 K u t n t n 1 u t n 1 + K u t n t n 1 u ( 3 ) t n 1 + K u ( 3 ) 0 u t n 3 K u 0 u t n + 3 K u 0 u t n K u 0 u ( 3 ) t n

From (A.5)

Q 1 1 , F 1 = 1 2 h F 1 t 0 + F 1 t 1 + 1 10 h 2 F 1 t 0 F 1 t 1 + 1 120 h 3 F 1 t 0 + F 1 t 1 = 1 2 h u t 0 t K t , t 0 t = t n u t n t 0 + K t n , t 0 u t n t 0 + u t 1 t K t , t 1 t = t n u t n t 1 + K t n , t 1 u t n t 1 + 1 10 h 2 u t 0 t K t , t 0 t = t n u t n t 0 + u t 0 K t n , t 0 u t n t 0 + u t 0 s t K t , t 0 t = t n u t n t 0 u t 0 t K t , t 0 t = t n u t n t 0 + u t 0 s K t n , t 0 u t n t 0 u t 0 K t n , t 0 u t n t 0 u t 1 t K t , t 1 t = t n u t n t 1 u t 1 K t n , t 1 u t n t 1 u t 1 s t K t , t 1 t = t n u t n t 1 + u t 1 t K t , t 1 t = t n u t n t 1 u t 1 s K t n , t 1 u t n t 1 + u t 1 K t n , t 1 u t n t 1 + 1 120 h 3 u t 0 t K t , t 0 u t n t 0 + u t 0 K t n , t 0 u t n t 0 + 2 u t 0 s t K t , t 0 u t n t 0 2 u t 0 t K t , t 0 u t n t 0 + 2 u t 0 s K t n , t 0 u t n t 0 2 u t 0 K t n , t 0 u t n t 0 + u t 0 2 s 2 t K t , t 0 u t n t 0 2 u t 0 s t K t , t 0 u t n t 0 + u t 0 t K t , t 0 u t n t 0 + u t 0 2 s 2 K t n , t 0 u t n t 0 2 u t 0 s K t n , t 0 u t n t 0 + u t 0 K t n , t 0 u ( 3 ) t n t 0 + u t 1 t K t , t 1 u t n t 1 + u t 1 K t n , t 1 u t n t 1 + 2 u t 1 s t K t , t 1 u t n t 1 2 u t 1 t K t , t 1 u t n t 1 + 2 u t 1 s K t n , t 1 u t n t 1 2 u t 1 K t n , t 1 u t n t 1 + u t 1 2 s 2 t K t , t 1 u t n t 1 2 u t 1 s t K t , t 1 u t n t 1 + u t 1 t K t , t 1 u t n t 1 + u t 1 2 s 2 K t n , t 1 u t n t 1 2 u t 1 s K t n , t 1 u t n t 1 + u t 1 K t n , t 1 u ( 3 ) t n t 1

with constant K t n , s , we obtain:

1 2 h u t 0 K u t n t 0 + u t 1 K u t n t 1 + 1 10 h 2 u t 0 K u t n t 0 u t 0 K u t n t 0 u t 1 K u t n t 1 + u t 1 K u t n t 1 + 1 120 h 3 u t 0 K u t n t 0 2 u t 0 K u t n t 0 + u t 0 K u ( 3 ) t n t 0 + u t 1 K u t n t 1 2 u t 1 K u t n t 1 + u t 1 K u ( 3 ) t n t 1 ,

Q 1 n , F 1 = 1 2 h F 1 t n 1 + F 1 t n + 1 10 h 2 F 1 t n 1 F 1 t n + 1 120 h 3 F 1 t n 1 + F 1 t n = 1 2 h u t n 1 t K t , t n 1 t = t n u t n t n 1 + K t n , t n 1 u t n t n 1 + u t n t K t , t n t = t n u 0 + K t n , t n u 0 + 1 10 h 2 u t n 1 t K t , t n 1 t = t n u t n t n 1 + u t n 1 K t n , t n 1 u t n t n 1 + u t n 1 s t K t , t n 1 t = t n u t n t n 1 u t n 1 t K t , t n 1 t = t n u t n t n 1 + u t n 1 s K t n , t n 1 u t n t n 1 u t n 1 K t n , t n 1 u t n t n 1 u t n t K t , t n t = t n u 0 u t n K t n , t n u 0 u t n s t K t , t n t = t n u 0 + u t n t K t , t n t = t n u 0 u t n s K t n , t n u 0 + u t n K t n , t n u 0 + 1 120 h 3 u t n 1 t K t , t n 1 t = t n u t n t n 1 + u t n 1 K t n , t n 1 u t n t n 1 + 2 u t n 1 s t K t , t n 1 t = t n u t n t n 1 2 u t n 1 t K t , t n 1 t = t n u t n t n 1 + 2 u t n 1 s K t n , t n 1 u t n t n 1 2 u t n 1 K t n , t n 1 u t n t n 1 + u t n 1 2 s 2 t K t , t n 1 t = t n u t n t n 1 2 u t n 1 s t K t , t n 1 t = t n u t n t n 1 + u t n 1 t K t , t n 1 t = t n u t n t n 1 + u t n 1 2 s 2 K t n , t n 1 u t n t n 1 2 u t n 1 s K t n , t n 1 u t n t n 1 + u t n 1 K t n , t n 1 u ( 3 ) t n t n 1 + u t n t K t , t n t = t n u 0 + u t n K t n , t n u 0 + 2 u t n s t K t , t n t = t n u 0 2 u t n t K t , t n t = t n u 0 + 2 u t n s K t n , t n u 0 2 u t n K t n , t n u 0 + u t n 2 s 2 t K t , t n t = t n u 0 2 u t n s t K t , t n t = t n u 0 + u t n t K t , t n t = t n u 0 + u t n 2 s 2 K t n , t n u 0 2 u t n s K t n , t n u 0 + u t n K t n , t n u ( 3 ) 0

with constant K t n , s , we obtain:

1 2 h u t n 1 K u t n t n 1 + u t n K u 0 + 1 10 h 2 u t n 1 K u t n t n 1 u t n 1 K u t n t n 1 u t n K u 0 + u t n K u 0 + 1 120 h 3 u t n 1 K u t n t n 1 2 u t n 1 K u t n t n 1 + u t n 1 K u ( 3 ) t n t n 1 + u t n K u 0 2 u t n K u 0 + u t n K u ( 3 ) 0

From (A.6)

Q 2 1 , F 2 = 1 2 h F 2 t 0 + F 2 t 1 + 1 12 h 2 F 2 t 0 F 2 t 1 = 1 2 h 2 t 2 K t , t 0 t = t n u t n t 0 u t 0 + 2 t K t , t 0 t = t n u t n t 0 u t 0 + K t n , t 0 u t n t 0 u t 0 + 2 t 2 K t , t 1 t = t n u t n t 1 u t 1 + 2 t K t , t 1 t = t n u t n t 1 u t 1 + K t n , t 1 u t n t 1 u t 1 + 1 12 h 2 s 2 t 2 K t , t 0 t = t n u t n t 0 u t 0 2 t 2 K t , t 0 t = t n u t n t 0 u t 0 + 2 t 2 K t , t 0 t = t n u t n t 0 u t 0 + 2 s t K t , t 0 t = t n u t n t 0 u t 0 2 t K t , t 0 t = t n u t n t 0 u t 0 + 2 t K t , t 0 t = t n u t n t 0 u t 0 + s K t n , t 0 u t n t 0 u t 0 K t n , t 0 u ( 3 ) t n t 0 u t 0 + K t n , t 0 u t n t 0 u t 0 s 2 t 2 K t , t 1 t = t n u t n t 1 u t 1 + 2 t 2 K t , t 1 t = t n u t n t 1 u t 1 2 t 2 K t , t 1 t = t n u t n t 1 u t 1 2 s t K ( t , s ) t = t n t , t 1 u t n t 1 u t 1 + 2 t K t , t 1 t = t n u t n t 1 u t 1 2 t K t , t 1 t = t n u t n t 1 u t 1 s K t n , t 1 u t n t 1 u t 1 + K t n , t 1 u ( 3 ) t n t 1 u t 1 K t n , t 1 u t n t 1 u t 1

with K t n , s constant, we obtain:

1 2 h K u t n t 0 u t 0 + K u t n t 1 u t 1 + 1 12 h 2 K u ( 3 ) t n t 0 u t 0 + K u t n t 0 u t 0 + K u ( 3 ) t n t 1 u t 1 K u t n t 1 u t 1

Q 2 n , F 2 = 1 2 h F 2 t n 1 + F 2 t n + 1 12 h 2 F 2 t n 1 F 2 t n = 1 2 h 2 t 2 K t , t n 1 t = t n u t n t n 1 u t n 1 + 2 t K t , t n 1 t = t n u t n t n 1 u t n 1 + K t n , t n 1 u t n t n 1 u t n 1 + 2 t 2 K t , t n t = t n u 0 u t n + 2 t K t , t n t = t n u 0 u t n + K t n , t n u 0 u t n + 1 12 h 2 s 2 t 2 t , t n 1 t = t n u t n t n 1 u t n 1 2 t 2 K t , t n 1 t = t n u t n t n 1 u t n 1 + 2 t 2 K t , t n 1 t = t n u t n t n 1 u t n 1

+ 2 s t K t , t n 1 t = t n u t n t n 1 u t n 1 2 t K t , t n 1 t = t n u t n t n 1 u t n 1 + 2 t K t , t n 1 t = t n u t n t n 1 u t n 1 + s K t n , t n 1 u t n t n 1 u t n 1 K t n , t n 1 u ( 3 ) t n t n 1 u t n 1 + K t n , t n 1 u t n t n 1 u t n 1 s 2 t 2 K t , t n t = t n u 0 u t n + 2 t 2 K t , t n t = t n u 0 u t n 2 t 2 K t , t n t = t n u 0 u t n 2 s t K t , t n t = t n u 0 u t n + 2 t K t , t n t = t n u 0 u t n 2 t K t , t n t = t n u 0 u t n s K t n , t n u 0 u t n + K t n , t n u ( 3 ) 0 u t n K t n , t n u 0 u t n

with constant K t n , s , we obtain:

1 2 h K u t n t n 1 u t n 1 + K u 0 u t n + 1 12 h 2 K u ( 3 ) t n t n 1 u t n 1 + K u t n t n 1 u t n 1 + K u ( 3 ) 0 u t n K u 0 u t n

From (A.4) we obtain:

v ( t n ) Q 1 , F 0 + Q n , F 0 + C 0 , w h e r e

Q 1 , F 0 = K 1 2 h u 0 + 3 28 h 2 u 0 + 1 84 h 3 u 0 + 1 1680 h 4 u ( 3 ) 0 u t n 3 28 h 2 u 0 + 2 84 h 3 u 0 + 3 1680 h 4 u 0 u t n + 1 84 h 3 u 0 + 3 1680 h 4 u 0 u t n 1 1680 h 4 u 0 u ( 3 ) t n + K 1 2 h u t n t 1 + 3 28 h 2 u t n t 1 + 1 84 h 3 u t n t 1 + 1 1680 h 4 u ( 3 ) t n t 1 u t 1 3 28 h 2 u t n t 1 + 2 84 h 3 u t n t 1 + 3 1680 h 4 u t n t 1 u t 1 + 1 84 h 3 u t n t 1 + 3 1680 h 4 u t n t 1 u t 1 1 1680 h 4 u t n t 1 u ( 3 ) t 1 ,

Q n , F 0 = K 1 2 h u 0 + 3 28 h 2 u 0 + 1 84 h 3 u 0 + 1 1680 h 4 u ( 3 ) 0 u t n 3 28 h 2 u 0 + 2 84 h 3 u 0 + 3 1680 h 4 u 0 u t n + 1 84 h 3 u 0 + 3 1680 h 4 u 0 u t n 1 1680 h 4 u 0 u ( 3 ) t n + K 1 2 h u t n t n 1 3 28 h 2 u t n t n 1 + 1 84 h 3 u t n t n 1 1 1680 h 4 u ( 3 ) t n t n 1 u t n 1 + 3 28 h 2 u t n t n 1 2 84 h 3 u t n t n 1 + 3 1680 h 4 u t n t n 1 u t n 1 + 1 84 h 3 u t n t n 1 3 1680 h 4 u t n t n 1 u t n 1 + 1 1680 h 4 u t n t n 1 u ( 3 ) t n 1 .

The v(t) expression represents a linear combination of four unknowns

u t n , u t n , u t n , u t n :

v ( t n ) A v 1 A v 2 A v 3 A v 4 u t n u t n u t n u t n + B v , w h e r e

A v 1 = 2 K 1 2 h u 0 + 3 28 h 2 u 0 + 1 84 h 3 u 0 + 1 1680 h 4 u ( 3 ) 0 , A v 2 = 2 K 3 28 h 2 u 0 + 2 84 h 3 u 0 + 3 1680 h 4 u 0 , A v 3 = 2 K 1 84 h 3 u 0 + 3 1680 h 4 u 0 , A v 4 = 2 K 1 1680 h 4 u 0 u ( 3 ) t n , B v = K 1 2 h u t n t 1 + 3 28 h 2 u t n t 1 + 1 84 h 3 u t n t 1 + 1 1680 h 4 u ( 3 ) t n t 1 u t 1 3 28 h 2 u t n t 1 + 2 84 h 3 u t n t 1 + 3 1680 h 4 u t n t 1 u t 1 + 1 84 h 3 u t n t 1 + 3 1680 h 4 u t n t 1 u t 1 1 1680 h 4 u t n t 1 u ( 3 ) t 1 + 1 2 h u t n t n 1 3 28 h 2 u t n t n 1 + 1 84 h 3 u t n t n 1 1 1680 h 4 u ( 3 ) t n t n 1 u t n 1 + 3 28 h 2 u t n t n 1 2 84 h 3 u t n t n 1 + 3 1680 h 4 u t n t n 1 u t n 1 + 1 84 h 3 u t n t n 1 3 1680 h 4 u t n t n 1 u t n 1 + 1 1680 h 4 u t n t n 1 u ( 3 ) t n 1 + C 0 ,

n ≥ 2.

From (A.5) we obtain:

v ( t n ) H ( t n , t n ) + Q 1 1 , F 1 + Q 1 n , F 1 + C 1 , w h e r e

H t , t = K u ( 0 ) u ( t n ) ;

Q 1 1 , F 1 = K 1 2 h u 0 + 1 10 h 2 u 0 + 1 120 h 3 u 0 u t n 1 10 h 2 u 0 + 2 120 h 3 u 0 u t n + 1 120 h 3 u 0 u ( 3 ) t n + K 1 2 h u t n t 1 + 1 10 h 2 u t n t 1 + 1 120 h 3 u ( 3 ) t n t 1 u t 1 1 10 h 2 u t n t 1 + 2 120 h 3 u t n t 1 u t 1 + 1 120 h 3 u t n t 1 u t 1 ; Q 1 n , F 1 = K 1 2 h u 0 + 1 10 h 2 u 0 + 1 120 h 3 u ( 3 ) 0 u t n 1 10 h 2 u 0 + 2 120 h 3 u 0 u t n + 1 120 h 3 u 0 u t n + K 1 2 h u t n t n 1 1 10 u t n t n 1 + 1 120 h 3 u ( 3 ) t n t n 1 u t n 1 + 1 10 h 2 u t n t n 1 2 120 h 3 u t n t n 1 u t n 1 + 1 120 h u t n t n 1 3 u t n 1 .

The v′(t) expression represents a linear combination of four unknowns

u t n , u t n , u t n , u t n :

v ( t n ) A v p 1 A v p 2 A v p 3 A v p 4 u t n u t n u t n u t n + B v p , w h e r e

A v p 1 = K u ( 0 ) + 1 2 h u 0 + 1 10 h 2 u 0 + 1 120 h 3 u ( 3 ) 0 ; A v p 2 = K 1 2 h u 0 + 1 10 h 2 u 0 + 1 120 h 3 u 0 1 10 h 2 u 0 + 2 120 h 3 u 0 ; A v p 3 = K 1 10 h 2 u 0 + 2 120 h 3 u 0 + 1 120 h 3 u 0 ; A v p 4 = K + 1 120 h 3 u 0 ; and B v p = K 1 2 h u t n t 1 + 1 10 h 2 u t n t 1 + 1 120 h 3 u ( 3 ) t n t 1 u t 1 1 10 h 2 u t n t 1 + 2 120 h 3 u t n t 1 u t 1 + 1 120 h 3 u t n t 1 u t 1 + 1 2 h u t n t n 1 1 10 u t n t n 1 + 1 120 h 3 u ( 3 ) t n t n 1 u t n 1 + 1 10 h 2 u t n t n 1 2 120 h 3 u t n t n 1 u t n 1 + 1 120 h u t n t n 1 3 u t n 1 + C 1 ,

n ≥ 2.

From (A.6) we obtain:

v ( t n ) d d t H ( t , t ) t = t n + t H ( t , s ) s = t = t n + Q 2 1 , F 2 + Q 2 n , F 2 + C 2 , w h e r e

d d t H ( t n , t n ) = K u ( 0 ) u ( t n )

t H t n , s s = t n = K u 0 u t n

Q 2 1 , F 2 = K 1 2 h u 0 + 1 12 h 2 u 0 u t n 1 12 h 2 u 0 u ( 3 ) t n + K 1 2 h u t n t 1 + 1 12 h 2 u ( 3 ) t n t 1 u t 1 1 12 h 2 u t n t 1 u t 1 ; Q 2 n , F 2 = K 1 2 h u 0 + 1 12 h 2 u ( 3 ) 0 u t n 1 12 h 2 u 0 u t n + K 1 2 h u t n t n 1 1 12 h 2 u ( 3 ) t n t n 1 u t n 1 + 1 12 h 2 u t n t n 1 u t n 1 .

The v″(t) expression represents a linear combination of four unknowns:

u t n , u t n , u t n , u t n :

v 2 p ( t n ) A v 2 p 1 A v 2 p 2 A v 2 p 3 A v 2 p 4 u t n u t n u t n u t n + B v 2 p , w h e r e

A v 2 p 1 = K u 0 + 1 2 h u 0 + 1 12 h 2 u ( 3 ) 0 , A v 2 p 2 = K u ( 0 ) 1 12 h 2 u 0 , A v 2 p 3 = K 1 2 h u 0 + 1 12 h 2 u 0 , A v 2 p 4 = K 1 12 h 2 u 0 , B v 2 p = K 1 2 h u t n t 1 + 1 12 h 2 u ( 3 ) t n t 1 u t 1 1 12 h 2 u t n t 1 u t 1 + 1 2 h u t n t n 1 1 12 h 2 u ( 3 ) t n t n 1 u t n 1 + 1 12 h 2 u t n t n 1 u t n 1 + C 2 , n 2 .

Appendix B: Extending the single-stage to a multi-stage

In this appendix, we give some further details that could be helpful with developing and implementing a code for the multistage system.

Adjust Eq. (28) to become

g m ( t ) = u m ( t ) + a m ( t ) u m ( t ) + b m ( t ) + b m , m 1 ( t ) u m 1 t + b m , m + 1 ( t ) u m + 1 t + 0 t k ( t , s ) u m ( t s ) u m ( s ) d s ,

where m denotes the stage number and the term b m , j ( t ) u j t represents the coupling between the stages m and j. For the first and last stages, we have

g 1 ( t ) = u 1 ( t ) + a 1 ( t ) u 1 ( t ) + b 1 ( t ) + b 1,2 ( t ) u 2 t + 0 t k ( t , s ) u 1 ( t s ) u 1 ( s ) d s , g M ( t ) = u M ( t ) + a M ( t ) u M ( t ) + b M ( t ) + b M , M 1 ( t ) u M 1 t + 0 t k ( t , s ) u M ( t s ) u M ( s ) d s .

Let c m , j t = b m , j ( t ) u j t , the first two derivatives are:

c m , j t = b m , j ( t ) u j t + b m , j ( t ) u j t , c m , j t = b m , j ( t ) u j t + 2 b m , j ( t ) u j t + b m , j ( t ) u j t .

Let c m , j t = b m , j ( t ) u j t and differentiate twice, we obtain the following system of equations for a stage m:

g m ( t ) = u m ( t ) + a m ( t ) u m ( t ) + b m ( t ) + c m , m 1 ( t ) + c m , m + 1 ( t ) + 0 t k ( t , s ) u m ( t s ) u m ( s ) d s g m ( t ) = f m ( t ) + a m ( t ) f m ( t ) + a m ( t ) f m ( t ) + b m ( t ) + c m , m 1 ( t ) + c m , m + 1 ( t ) + H m ( t , t ) + 0 t t H m ( t , s ) d s g m ( t ) = f m ( t ) + a m ( t ) f m ( t ) + 2 a m ( t ) f m ( t ) + a m ( t ) f m ( t ) + b m ( t ) + c m , m 1 ( t ) + c m , m + 1 ( t ) + d d t H m ( t , t ) + t H m ( t , s ) s = t + 0 t 2 t 2 H m ( t , s ) d s .

By extending the development of the linear system representation in Section 3.3 to the multi-stage system represented by Eq. (11) we obtain a general matrix with the following pattern:

u m t n u m t n u m t n u m t n
g m−1 b m−1,m (t n ) 0 0 0
g m 1 b m 1 , m ( t n ) b m−1,m (t n ) 0 0
g m 1 b m 1 , m ( t n ) 2 b m 1 , m ( t n ) b m−1,m (t n ) 0
hc m−1 0 0 0 0
g m Ag 1 Ag 2 Ag 3 Ag 4
g m Agp 1 Agp 2 Agp 3 Agp 4
g m Ag2p 1 Ag2p 2 Ag2p 3 Ag2p 4
hc m Ahcf 1 Ahcf 2 Ahcf 3 Ahcf 4
g m+1 b m+1,m (t n ) 0 0 0
g m + 1 b m + 1 , m ( t n ) b m+1,m (t n ) 0 0
g m + 1 b m + 1 , m ( t n ) 2 b m + 1 , m ( t n ) b m+1,m (t n ) 0
hc m+1 0 0 0 0
hc denotes the Hermite correction formula.

The linear system AU = B extends for the M-stage system as follows:

Denoting the stage number using m, we rename A to A m , B to B m , and U to U m . Then, using the development done for the multi-stage system, we define the following matrices:

C m = b m 1 , m ( t n ) 0 0 0 b m 1 , m ( t n ) b m 1 , k ( t n ) 0 0 b m 1 , m ( t n ) 2 b m 1 , m ( t n ) b m 1 , m ( t n ) 0 0 0 0 0 , m = 2 , , M D m = b m + 1 , m ( t n ) 0 0 0 b m + 1 , m ( t n ) b m + 1 , m ( t n ) 0 0 b m + 1 , m ( t n ) 2 b m + 1 , m ( t n ) b m + 1 , m ( t n ) 0 0 0 0 0 , m = 1 , , M 1 .

We obtain the multi-stage linear system AU = B, where A is a block tridiagonal matrix of size 4M × 4M, B and U are column vectors of size 4M × 1.

A = A 1 C 2 0 D 1 A 2 C 3 D m 1 A m C m + 1 D M 2 A M 1 C M 0 D M 1 A M , B = B 1 B M , U = U 1 U M .

Appendix C: Convergence

In this appendix, we investigate the convergence of the method. For simplicity and without loss of generality, we consider a 3-stage system with constant coefficients. We consider first the method applied to a 3-stage system of the linear ode part associated with our system, that is, we set k i ≡ 0 in (11) We have

g 1 ( t ) = u 1 ( t ) + a 1 u 1 ( t ) + b 1,2 u 2 t g 2 ( t ) = u 2 ( t ) + a 2 u 2 ( t ) + b 2,3 u 3 t + b 2,1 u 1 t g 3 ( t ) = u k ( t ) + a 3 u 3 ( t ) + b 3,2 u 2 t .

We assume the uniqueness of solutions for the initial value problems. Let u n , j , u n , j , u n , j , u n , j be the approximation for the theoretical solutions u j t n , u j t n , u j t n , u j t n , respectively, j = 1, 2, 3. Now define the errors

ε n , 1 = u 1 t n u n , 1 , ε n , 1 = u 1 t n u n , 1 , ε n , 1 = u 1 t n u n , 1 , ε n , 1 = u 1 t n u n , 1 , ε n , 2 = u 2 t n u n , 2 , ε n , 2 = u 2 t n u n , 2 , ε n , 2 = u 2 t n u n , 2 , ε n , 2 = u 2 t n u n , 2 , ε n , 3 = u 3 t n u n , 3 , ε n , 3 = u 3 t n u n , 3 , ε n , 3 = u 3 t n u n , 3 , ε n , 3 = u 3 t n u n , 3 .

This Milne–Obrechkoff predictor–corrector type method for odes leads to the following systems of equations and the corresponding error equations. We ignore at first the error in the Hermite corrector.

g 1 ( t n ) = u 1 ( t n ) + a 1 u 1 ( t n ) + b 1,2 u 2 t n g 1 ( t n ) = u 1 ( t n ) + a 1 u 1 ( t n ) + b 1,2 u 2 t n g 1 ( t n ) = u 1 ( t n ) + a 1 u 1 ( t n ) + b 1,2 u 2 t n u 1 t n = u 1 t n 1 + h 2 u 1 t n 1 + u 1 t n + h 2 10 u 1 t n 1 u 1 t n + h 3 120 u 1 t n 1 + u 1 t n , ε n , 1 a 1 ε n , 1 + b 1,2 ε n , 2 ε n , 1 a 1 ε n , 1 + b 1,2 ε n , 2 ε n , 1 a 1 ε n , 1 + b 1,2 ε n , 2 ε n , 1 ε n 1,1 + h 2 ε n 1,1 + ε n , 1 + h 2 10 ε n 1,1 + ε n , 1 + h 3 120 ε n 1,1 + ε n , 1 . g 2 ( t n ) = u 2 ( t n ) + a 2 u 2 ( t n ) + b 2,3 u 3 t n + b 2,1 u 1 t n g 2 ( t n ) = u 2 ( t n ) + a 2 u 2 ( t n ) + b 2,3 u 3 t n + b 2,1 u 1 t n g 2 ( t n ) = u 2 ( t n ) + a 2 u 2 ( t n ) + b 2,3 u 3 t n + b 2,1 u 1 t n u 2 t n = u 2 t n 1 + h 2 u 2 t n 1 + u 2 t n + h 2 10 u 2 t n 1 u 2 t n + h 3 120 u 2 t n 1 + u 2 t n , ε n , 2 a 2 ε n , 2 + b 2,3 ε n , 3 + b 2,1 ε n , 1 ε n , 2 a 2 ε n , 2 + b 2,3 ε n , 3 + b 2,1 ε n , 1 ε n , 2 a 2 ε n , 2 + b 2,3 ε n , 3 + b 2,1 ε n , 1 ε n , 2 ε n 1,2 + h 2 ε n 1,2 + ε n , 2 + h 2 10 ε n 1,2 + ε n , 2 + h 3 120 ε n 1,2 + ε n , 2 . g 3 ( t n ) = u 3 ( t n ) + a 3 u 3 ( t n ) + b 3,2 u 2 t n g 3 ( t n ) = u 3 ( t n ) + a 3 u 3 ( t n ) + b 3,2 u 2 t n g 3 ( t n ) = u 3 ( t n ) + a 3 u 3 ( t n ) + b 3,2 u 2 t n u 3 t n = u 3 t n 1 + h 2 u 3 t n 1 + u 3 t n + h 2 10 u 3 t n 1 u 3 t n + h 3 120 u 3 t n 1 + u 3 t n , ε n , 3 a 3 ε n , 3 + b 3,2 ε n , 2 , ε n , 3 a 3 ε n , 3 + b 3,2 ε n , 2 , ε n , 3 a 3 ε n , 3 + b 3,2 ε n , 2 , ε n , 3 ε n 1,3 + h 2 ε n 1,3 + ε n , 3 + h 2 10 ε n 1,3 + ε n , 3 + h 3 120 ε n 1,3 + ε n , 3 .

Adding the errors from the first 3 predictor equations in each system, we obtain

L H S = ε n , 1 + ε n , 1 + ε n , 1 + ε n , 2 + ε n , 2 + ε n , 2 + ε n , 3 + ε n , 3 + ε n , 3 ,

R H S = a 1 + b 2,1 ε n , 1 + ε n , 1 + ε n , 1 + b 1,2 + a 2 + b 3,2 ε n , 2 + ε n , 2 + ε n , 2 + a 3 + b 2,3 ε n , 3 + ε n , 3 + ε n , 3 .

Since we solve the approximating equations simultaneously (predictor and corrector—each is applied once), we add the errors from the corrector, and use bounds on the previous step estimates

ε n , j ε n 1 , j + h 2 ε n 1 , j + ε n , j + h 2 10 ε n 1 , j + ε n , j + h 3 120 ε n 1 , j + ε n , j E j ε n , j + ε n , j + ε n , j + F j ε n 1 , j + ε n 1 , j + ε n 1 , j + ε n 1 , j

where E j , F j are constants. Hence,

L H S = j = 1 3 ε n , j + ε n , j + ε n , j + ε n , j ,

and

R H S = a 1 + b 2,1 ε n , 1 + ε n , 1 + ε n , 1 + b 1,2 + a 2 + b 3,2 ε n , 2 + ε n , 2 + ε n , 2 + a 3 + b 2,3 ε n , 3 + ε n , 3 + ε n , 3 + j = 1 3 E j ε n , j + ε n , j + ε n , j + F j ε n 1 , j + ε n 1 , j + ε n 1 , j + ε n 1 , j A j = 1 3 ε n , j + ε n , j + ε n , j + ε n , j + O h 6 ,

where A is a large enough constant, and ε 0 , i = ε 0 , i = ε 0 , i = ε 0 , i = 0 .

Now consider the generic stage equation

g j ( t ) = u j ( t ) + a j ( t ) u j ( t ) + b j , j + 1 ( t ) u j + 1 t + b j , j 1 t u j 1 t + 0 t k j ( t , s ) u j ( t s ) u j ( s ) d s .

Under our general smoothness assumptions on 0 , T , which gives us uniform bounds on the coefficients, we can for simplicity consider the convergence analysis of a generic stage with constant coefficients, which we write as

g j ( t ) = u j t + a j u t + b j , j + 1 u j + 1 t + b j , j 1 u j 1 t + k j 0 t u j ( t s ) u j ( s ) d s ,

where a i , b j,j+1, b j,j−1, k j are constants. For 0 ≤ in, let ε i , j = u j t i u i , j , ε i , i = u j t i u i , j , ε i , j = u j t i u i , j , ε i , j = u j t i u i , j . Obviously, ε 0 , j = 0 , ε 0 , j = 0 , ε 0 , j = 0 , ε 0 , j = 0 . We have

(C.1) g j ( t n ) = u j ( t n ) + a j u j t n + b j , j + 1 u j + 1 t n + b j , j 1 u j 1 t n + k j 0 t n u j ( t n s ) u j ( s ) d s u j ( t n ) + a j u j t n + b j , j + 1 u j + 1 t n + b j , j 1 u j 1 t n + 1 i n 1 2 h F 0 t i 1 + F 0 t i + 3 28 h 2 F 0 t i 1 F 0 t i + 1 84 h 3 F 0 t i 1 + F 0 t i + 1 1680 h 4 F 0 t i 1 F 0 t i = u j ( t n ) + a j u j t n + b j , j + 1 u j + 1 t n + b j , j 1 u j 1 t n + h 0 i n α i u ( t n t i ) u ( t i ) + h 2 0 i n β i s u ( t n t i ) u ( t i ) + h 3 0 i n γ i 2 s 2 u ( t n t i ) u ( t i ) + h 4 0 i n δ i 3 s 3 u ( t n t i ) u ( t i ) ,

where α i , β i , γ i , δ i , 1 ≤ in are constant coefficients. In a similar manner, we get

(C.2) g j ( t n ) = u n , j + a j u n , j + b j , j + 1 u n , j + 1 + b j , j 1 u n , j 1 + h 0 i n α i u n i , j u i , j + h 2 0 i n β i s u n i , j u i , j + h 3 0 i n γ i 2 s 2 u n i , j u i , j + h 4 0 i n δ i 3 s 3 u n i , j u i , j .

Then (C.1) and (C.2) yield

(C.3) ε j , n = u j ( t n ) u n , j a j u j t n u n , j + b j , j + 1 u j + 1 t n u n , j + 1 + b j , j 1 u j 1 t n u n , j 1 + h 0 i n α i u j ( t n t i ) u j ( t i ) u n i , j u i , j + h 2 0 i n β i s u j ( t n t i ) u j ( t i ) s u n i , j u i , j + h 3 0 i n γ i 2 s 2 u j ( t n t i ) u j ( t i ) 2 s 2 u n i , j u i , j + h 4 0 i n δ i 3 s 3 u j ( t n t i ) u j ( t i ) 3 s 3 u n i , j u i , j + h 2 F 0 2 q T h 2 q 2 q + 1 q ! 2 q ! 2 , q = 4 .

Then we obtain 0 i n α i u ( t n t i ) u ( t i ) u n i u i C 0 i n ε i + ε n i , where C is a positive constant, as

α i u j t n t i u j t i u n i , j u i , j = α i u j t n t i u j t i u j t n t i u i , j + u j t n t i u i , j u n i , j u i , j α i u j t n t i u j t i u i , j + u i , j u j t n t i u n i , j C 1 u j t i u i , j + C 2 u j t n t i u n i , j C 3 ε i , j + ε n i , j ,

where C 1, C 2, C 3 are positive constants. Similarly we obtain for the fourth term in (C.3)

β i u j t n t i u j t i u n i , j u i , j + u j t n t i u j t i u n i , j u i , j C 4 ε i , j + ε n i , j + C 5 ε i , j + ε n i , j ,

where C 4, C 5 are positive constants.

Next, in a similar manner, we obtain for the fifth term in (C.3)

γ i 2 s 2 u j ( t n t i ) u j ( t i ) 2 s 2 u n i , j u i , j L 1 ε n i , j + L 2 ε i , j + L 3 ε n i , j + L 4 ε i , j + L 5 ε n i , j + L 6 ε i , j ,

where L m , m = 1, …, 6 are positive constants; and for the sixth term in (C.3)

δ i 3 s 3 u j ( t n t i ) u j ( t i ) 3 s 3 u n i , j u i , j M 1 ε i , j + M 2 ε n i , j + M 3 ε i , j + M 4 ε n i , j + M 5 ε i , j + M 6 ε n i , j + M 7 ε i , j + M 8 ε n i , j

where M m , m = 1, …, 8, are positive constants. Noting that i = 0 n ε i , j = 0 i n ε n i , j , and similarly for ɛ replaced by ɛ′, ɛ″, and ɛ‴, then (C.3) yields

ε n , j N 1 h 8 + a j ε n , j + b 1 , j ε n . j + 1 + b 2 , j ε n , j 1 + Q 1 h 0 i n ε i , j + ε i , j + ε i , j + ε i , j ,

where N 1 and Q 1 are positive constants.

Applying a similar analysis to the equations of g j , g j , and u j , we obtain,

ε n , j N 2 h 6 + a j ε n , j + b 1 , j ε n , j + 1 + b 2 , j ε n , j 1 + Q 2 h 0 i n ε i , j + ε i , j + ε i , j + ε i , j , ε n , j N 3 h 4 + a j ε n , j + b 1 , j ε n , j + 1 + b 2 , j ε n , j 1 + Q 3 h 0 i n ε i , j + ε i , j + ε i , j + ε i , j , ε n , j N 0 h 6 + Q 0 h 0 i n ε i , j + ε i , j + ε i , j + ε i , j ,

where N 0, N 2, N 3, Q 0, Q 2, and Q 3 are positive constants.

Put j = 2, and then obtain analogously the estimates for the first stage j = 1 , and those for the last stage j = 3 . In view of the errors we have derived for the ode systems above, we see that after summing up the errors from the 3 stages, we obtain

j = 1 3 ε n , j + ε n , j + ε n , j + ε n , j N h 4 + R h j = 1 3 0 i n ε i , j + ε i , j + ε i , j + ε i , j ,

where N, Q, and R are positive constants.

By the discrete version of Gronwall’s inequality, we obtain

j = 1 3 ε n , j + ε n , j + ε n , j + ε n , j N h 4 1 R h exp R T 1 R h ,

which clearly indicates that the errors go to zero as h → 0 (with nh = T).

References

[1] A. D. Randolph, C. Deepak, and M. Iskander, “On the narrowing of particle-size distribution in staged vessels with classified product removal,” AIChE J., vol. 14, pp. 827–830, 1968. https://doi.org/10.1002/aic.690140532.Search in Google Scholar

[2] A. D. Randolph and C. S. Tan, “Numerical design techniques for staged classification recycle crystallizers,” Ind. Eng. Chem. Process Des. Dev., vol. 17, pp. 189–200, 1978. https://doi.org/10.1021/i260066a013.Search in Google Scholar

[3] K. Tamura and H. Tsuge, “Characteristics of multistage column crystallizer for gas–liquid reactive crystallization of calcium carbonate,” Chem. Eng. Sci., vol. 61, pp. 5818–5826, 2006. https://doi.org/10.1016/j.ces.2006.05.002.Search in Google Scholar

[4] R. W. Farmer and J. R. Beckman, “Particle size improvement by countercurrent tower crystallizer,” AIChE J., vol. 32, pp. 1099–1106, 1986. https://doi.org/10.1002/aic.690320706.Search in Google Scholar

[5] R. W. Farmer and J. R. Beckman, “Evaluation of crystal size distribution enlargement in a multistage column crystallizer,” Ind. Eng. Chem. Res., vol. 30, pp. 196–201, 1991. https://doi.org/10.1021/ie00049a029.Search in Google Scholar

[6] H. Brunner and P. J. van der Houwen, The Numerical Solution of Volterra Equations, CWI Monographs 3, Amesterdam, North-Holland, 1986.Search in Google Scholar

[7] L. M. Delves and J. L. Mohamed, Computational Methods for Integral Equations, Cambridge, Cambridge University Press, 1988.Search in Google Scholar

[8] S. H. Chung, H. W. Lee, and E. G. Saleeby, “Modeling a cascade of continuous MSMPR crystallizers with agglomeration,” Chem. Eng. Commun., vol. 63, pp. 177–200, 1998. https://doi.org/10.1080/00986449808912350.Search in Google Scholar

[9] K. El Khaldi, M. Khasawneh and E. G. Saleeby, “A Rothe-modified shooting method for solving a nonlinear boundary value problem arising in particulate processes,” Comput. Appl. Math., vol. 38, pp. 12, 2019. https://doi.org/10.1007/s40314-019-0780-1.Search in Google Scholar

[10] K. El Khaldi, N. Mourany and E. G. Saleeby, “On the identification of parameters for a nonlinear integrodifferential population balance equation,” Int. J. Comput. Math., vol. 90, no. 9, pp. 2019–2035, 2013. https://doi.org/10.1080/00207160.2013.772593.Search in Google Scholar

[11] H. M. Hulbert and S. Katz, “Some problems in particle technology,” Chem. Eng. Sci., vol. 19, pp. 555–574, 1964. https://doi.org/10.1016/0009-2509(64)85047-8.Search in Google Scholar

[12] B. D. Khanh, “Hermite predictor–corrector scheme for regular Volterra integral equations and for some integro-differential equations for turbulent diffusion,” J. Comput. Appl. Math., vol. 51, pp. 305–316, 1994. https://doi.org/10.1016/0377-0427(92)00012-x.Search in Google Scholar

[13] B. A. Velkison, “Solution of a nonlinear integro-differential equation,” USSR Comput. Math. Math. Phys., vol. 65, pp. 256–259, 1975.10.1016/0041-5553(75)90157-3Search in Google Scholar

[14] Y. Wei and T. Tao, “The numerical analysis of implicit Runge–Kutta methods for a certain nonlinear integro-differential equation,” Math. Comput., vol. 54, pp. 155–168, 1990. https://doi.org/10.2307/2008687.Search in Google Scholar

[15] S. H. Chang and J. T. Day, “On the numerical solution of a certain non-linear integro-differential equation,” J. Comput. Phys., vol. 26, pp. 162–168, 1978. https://doi.org/10.1016/0021-9991(78)90088-8.Search in Google Scholar

[16] H. Brunner, “Implicit Runge–Kutta methods of optimal order for Volterra integro-differential equations,” Math. Comput., vol. 42, pp. 95–109, 1984. https://doi.org/10.1090/s0025-5718-1984-0725986-6.Search in Google Scholar

[17] T. Tao and Y. Wei, “The further study of a certain nonlinear integro-differential equation,” J. Comput. Phys., vol. 72, pp. 486–497, 1987. https://doi.org/10.1016/0021-9991(87)90095-7.Search in Google Scholar

[18] A. D. Randolph and M. A. Larson, Theory of Particulate Processes, 2nd ed., New York, Academic Press, 1988.10.1016/B978-0-12-579652-1.50007-7Search in Google Scholar

[19] C. T. H. Baker, “A perspective on the numerical treatment of Volterra equations,” J. Comput. Appl. Math., vol. 125, pp. 217–249, 2000. https://doi.org/10.1016/s0377-0427(00)00470-2.Search in Google Scholar

[20] P. Linz, Analytical and Numerical Methods for Volterra Equations, Studies in Applied Numerical Mathematics, SIAM, Philadelphia, 1985.10.1137/1.9781611970852Search in Google Scholar

[21] C. Lubich, “Runge–Kutta theory for Volterra integrodifferential equations,” Numer. Math., vol. 40, pp. 119–135, 1982. https://doi.org/10.1007/bf01459081.Search in Google Scholar

[22] C. A. Dorao and H. A. Jakobsen, “hp-adaptive least squares spectral element method for population balance equations,” Appl. Numer. Math., vol. 58, pp. 563–576, 2008. https://doi.org/10.1016/j.apnum.2006.12.005.Search in Google Scholar

[23] F. Gelbard and J. H. Seinfeld, “Numerical solution of the dynamic equation for particulate systems,” J. Comput. Phys., vol. 28, pp. 357–376, 1978. https://doi.org/10.1016/0021-9991(78)90058-x.Search in Google Scholar

[24] D. Ramkrishna, Population Balances, Theory and Application to Particulate Systems in Engineering, SanDiego, Academic Press, 2000.Search in Google Scholar

[25] S. K. Bhatia and D. Chakraborty, “Modified MWR approach: application to agglomerative precipitation,” AIChE J., vol. 38, pp. 868–878, 1992. https://doi.org/10.1002/aic.690380608.Search in Google Scholar

[26] Y. Liu and T. Cameron, “A new wavelet-based method for the solution of the population balance,” Chem. Eng. Sci., vol. 56, pp. 5283–5294, 2001. https://doi.org/10.1016/s0009-2509(01)00196-8.Search in Google Scholar

[27] L. F. L. R. Silva, R. C. Rodrigues, J. F. Mitre, and P. L. C. Lage, “Comparison of the accuracy and performance of quadrature-based methods for population balance problems with simultaneous breakage and aggregation,” Comput. Chem. Eng., vol. 34, pp. 286–297, 2010. https://doi.org/10.1016/j.compchemeng.2009.11.005.Search in Google Scholar

[28] N. S. Tavare, J. Garside, and M. A. Larson, “Crystal size distribution from a cascade of MSMPR crystallizers with magma recycle,” Chem. Eng. Commun., vol. 47, pp. 185–199, 1986. https://doi.org/10.1080/00986448608911763.Search in Google Scholar

Received: 2021-10-21
Revised: 2021-08-22
Accepted: 2021-11-04
Published Online: 2021-12-02

© 2021 Walter de Gruyter GmbH, Berlin/Boston

Articles in the same Issue

  1. Frontmatter
  2. Original Research Articles
  3. Numerical solutions for variable-order fractional Gross–Pitaevskii equation with two spectral collocation approaches
  4. Threshold dynamics of an HIV-1 model with both virus-to-cell and cell-to-cell transmissions, immune responses, and three delays
  5. Numerical scheme and analytical solutions to the stochastic nonlinear advection diffusion dynamical model
  6. On modeling column crystallizers and a Hermite predictor–corrector scheme for a system of integro-differential equations
  7. On a discrete fractional stochastic Grönwall inequality and its application in the numerical analysis of stochastic FDEs involving a martingale
  8. A new self adaptive Tseng’s extragradient method with double-projection for solving pseudomonotone variational inequality problems in Hilbert spaces
  9. Solitary waves of the RLW equation via least squares method
  10. Optical solitons and stability regions of the higher order nonlinear Schrödinger’s equation in an inhomogeneous fiber
  11. Weighted pseudo asymptotically Bloch periodic solutions to nonlocal Cauchy problems of integrodifferential equations in Banach spaces
  12. Modified inertial subgradient extragradient method for equilibrium problems
  13. Propagation of dark-bright soliton and kink wave solutions of fluidized granular matter model arising in industrial applications
  14. Existence and uniqueness of positive solutions for fractional relaxation equation in terms of ψ-Caputo fractional derivative
  15. Investigation of nonlinear fractional delay differential equation via singular fractional operator
  16. Travelling peakon and solitary wave solutions of modified Fornberg–Whitham equations with nonhomogeneous boundary conditions
  17. The (3 + 1)-dimensional Wazwaz–KdV equations: the conservation laws and exact solutions
  18. Stability and ψ-algebraic decay of the solution to ψ-fractional differential system
  19. Some new characterizations of a space curve due to a modified frame N , C , W in Euclidean 3-space
  20. Numerical simulation of particulate matter propagation in an indoor environment with various types of heating
  21. Inertial accelerated algorithms for solving split feasibility with multiple output sets in Hilbert spaces
  22. Stability analysis and abundant closed-form wave solutions of the Date–Jimbo–Kashiwara–Miwa and combined sinh–cosh-Gordon equations arising in fluid mechanics
  23. Contrasting effects of prey refuge on biodiversity of species
Downloaded on 18.9.2025 from https://www.degruyterbrill.com/document/doi/10.1515/ijnsns-2020-0239/html
Scroll to top button