Acta Polytechnica https://doi.org/10.14311/AP.2025.65.0566 Acta Polytechnica 65(5):566–577, 2025 © 2025 The Author(s). Licensed under a CC-BY 4.0 licence Published by the Czech Technical University in Prague METHOD OF LINES FOR REACTION-DIFFUSION SYSTEMS ADMITTING INVARIANT REGIONS Niels van der Meer, Michal Beneš∗ Czech Technical University in Prague, Faculty of Nuclear Sciences and Physical Engineering, Department of Mathematics, Trojanova 13, 120 00 Prague, Czech Republic ∗ corresponding author: michal.benes@fjfi.cvut.cz Abstract. Systems of nonlinear reaction-diffusion equations arise in various fields, including chemistry, population dynamics, pattern formation, phase transitions, and image processing. With the exception of few analytically solvable cases, they are treated by numerical methods carefully adjusted to capture the nonlinear phenomena exhibited by the solution. This article shows how to extend the notion of invariant regions generalizing the maximum principle for diffusion equations to the finite-difference method of lines, and how to consequently prove convergence of the underlying numerical scheme. We also provide two particular examples of reaction-diffusion systems in one-dimensional space with a diagonal diffusion operator, which are solved by the presented numerical method. Keywords: Reaction-diffusion dynamics, finite-difference method, method of lines, invariant regions. 1. Introduction Reaction-diffusion systems arise in a variety of natural science disciplines. Originally, the reaction of several chemical components in a volume of space is described in terms of the mass balance by a system of reaction equations with spatial transport by diffusion. This is summarized, for example in [1–3]. Similarly, such sys- tems describe the competitive coexistence of several biological species within a spatial volume (see, e.g. [4– 6]). Phase transitions such as those during metal solid- ification at microscale, or the conversion of crystalline lattice can also lead to systems of reaction-diffusion equations known as phase-field [7–10]. Geometric methods of image processing based on the motion of level sets employ a reaction-diffusion equation of the Allen-Cahn type [11]. Numerical solution of such systems is being per- formed using a variety of methods. The full finite- difference discretization is summarized, for exam- ple, in [2, 12], the finite-element method is used e.g. in [13, 14], the non-linear spectral Galerkin method is discussed, for example, in [15, 16]. This article is devoted to the method of lines based on space discretization by the finite-difference method and subsequent treatment of the ODE system in time. We show that the property of invariant regions prop- agates to the same property for such a system of ODE’s generated by the finite-difference method and allows to show sufficient a priori estimates needed for convergence of the numerical scheme. The text is organized into a mathematical introduc- tion with the state of the art, design of the method of lines, a demonstration of the invariant-region property, convergence analysis and applications to two differ- ent systems of reaction-diffusion equations, including computational results. We also indicate possible gen- eralizations of the presented approach. 2. Reaction-diffusion systems In the following, we consider an initial-boundary-value problem for a system of d reaction-diffusion equations on a space interval (a, b) and a time interval (0, T ) in the form: ∂tU = D∂2 xxU + F(U) in (0, T ) × (a, b), U |x=a = 0, U |x=b = 0, U |t=0 = Uini, (1) where D ∈ Rd×d, d ∈ N, is the positive definite, for simplicity diagonal, matrix, F : Rd → Rd is the vector- valued C(1)(Rd) map, U : [0, T ] × [a, b] → Rd is the solution. As discussed in literature (see e.g. [1, 3, 13, 15, 17]), the system (1) is studied using the weak formula- tion. For this purpose, we introduce the Lebesgue space H = L2((a, b);Rd) of square-integrable vector functions U, V , with the scalar product: (U, V ) = d∑ i=1 ∫ b a ui(x)vi(x)dx, and the Sobolev space of vector-valued functions V = W̊ (1) 2 ((a, b);Rd) with the scalar product: (U, V )V = d∑ i=1 ∫ b a dui(x) dx dvi(x) dx dx, for U = [u1, . . . , ud]T, V = [v1, . . . , vd]T. As in literature (see e.g. [3, 15]), problem (1) has the weak solution U : (0, T ) → V, provided: d dt (U, V ) + (DU, V )V = (F(U), V ), for all V ∈ V, in the sense of D′((0, T )), (2) U |t=0 = Uini. 566 https://doi.org/10.14311/AP.2025.65.0566 https://creativecommons.org/licenses/by/4.0/ https://www.cvut.cz/en vol. 65 no. 5/2025 Method of lines for reaction-diffusion systems . . . Remark In the given framework, the weak solution, provided it exists, is in fact a continuous map, due to the corresponding embedding of the Sobolev space (see [3]). The existence as well as convergence of the numer- ical solution are studied using suitable bounds for norms of the solution (see [10, 15]). One way to ob- tain them is based on a contractive property of the reaction field (F(U), V ), known as the invariant region (see e.g. [1]). Definition 1. The system (1) is said to possess an invariant region O ⊂ Rd, provided O is a bounded closed convex set, and for each initial condition Uini ∈ C((a, b);Rd), with values in O, the solution values U(t, x) also are in O for all x ∈ (a, b) and all t > 0 for which the solution exists. We aim to show that such a property valid for the numerical solution guarantees the convergence of the numerical solution to the weak solution. 3. Method of lines 3.1. Finite-difference discretization Here, we summarize the use of the Finite-Difference Method for space discretization of (1), while the time variable remains continuous. This widely used ap- proach is known as the method of lines (see [2, 5, 15, 18– 20]). For this purpose we introduce corresponding notations (as in [11, 21, 22]). Using m ∈ N as the number of meshes, and h = b−a m as the mesh size, we denote ωh = {a + jh | j = 0, . . . ,m} the finite- difference grid, ωh := {a + jh | j = 1, . . . ,m − 1} the internal nodes. We abbreviate the notation of values of functions V : ωh → Rd defined on ωh – grid functions – as Vj = V (a+ jh). The space of all grid functions with zero boundary values is denoted as Hh = {V : ωh → Rd | V0 = 0, Vm = 0}. Discretization of the second space derivatives is obtained by the central difference with corresponding error of approximation O(h2). Vx,j = Vj+1 − Vj h , Vx̄,j = Vj − Vj−1 h , Vx̄x,j = Vj+1 − 2Vj + Vj−1 h2 . Next, for the grid functions V,W ∈ Hh, we denote the products and norms on the grid: (V,W )h = m−1∑ j=1 hvjwj , ∥V ∥h = √ (V, V )h, (V,W ] = m∑ j=1 hvjwj , ||V ]| = √ (V, V ]. We notice that the integration by parts can be used: (V,Wx)h = vmwm − v0w1 − (Vx̄,W ] , which obviously simplifies on Hh. We also recall the maximum norm: ∥V ∥0,h = max{|vj | | j = 0, . . . ,m}, and the embedding inequalities discussed in [11, 22]: ∥V ∥0,h ≤ C(a, b)||V ]|. (3) We then discretize (1) and design the semi-discrete scheme: dZ dt = DZx̄x + F(Z), in (0, T ) × ωh, Z|t=0= Zini, in ωh, (4) Z0 = 0, Zm = 0, in (0, T ), whose solution is a time dependent vector-valued func- tion Z = Z(t) with values in Hh. 3.2. Invariant regions for the semi-discrete scheme We now show that the idea of an invariant region for the original system (1) can be extended to semi- discrete scheme (4). For simplicity, we consider a pris- matic shape of the invariant region, even though a general shape as in [1] can be treated correspond- ingly. To this end, we denote, for given constants Ak < Bk, k = 1, . . . , d, the set: Σ = {W ∈ Rd|Ak ≤ wk ≤ Bk, k = 1, . . . , d}. We say that Σ is the invariant region for system (4), provided for each initial condition Zini : ωh → Rd having values in Σ, the solution of (4) available for t ∈ [0, T ) (including T = +∞) also has values in Σ. Consequently, we formulate the statement related to it. Lemma 1. Consider the semi-discrete scheme (4) and the region Σ. If F points strictly into Σ along ∂Σ, then Σ is the invariant region. Proof. We denote F(V ) = [F 1(V ), . . . , F d(V )], D = diag(D1, . . . , Dd). By contradiction, we assume that Zini has values in Σ, however, there is t1 ∈ (0, T ) at which the solution Z(t) of (4) leaves Σ. Assume that this happens on the boundary of Σ given by the value Zk0 = Bk0 for a specific k0. As F k0(Z) < 0 along this edge of Σ, there is a band of values of thickness δ > 0 such that Fk0(Z) < 0 for Zk0 ∈ [Bk0 , Bk0 + δ]. This means that there is a time t0 ≤ t1 a node j0 for which: Zk0 j0 (t0) > Bk0 , and such a value Zk0 j0 (t0) is maximal on ωh, i.e. Zk0 j0 (t0) ≥ Zk0 j (t0), j = 0, . . . ,m. The central second difference Zk0 x̄x(t0) is, therefore, non-positive: Zk0 x̄x,j0 (t0) = Zk0 j0+1(t0) − 2Zk0 j0 (t0) + Zk0 j0−1(t0) h2 ≤ 0. As Zk0 j0 ∈ [Bk0 , Bk0 + δ], the inward pointing field F satisfies: F k0(Zj0(t0)) < 0. 567 Niels van der Meer, Michal Beneš Acta Polytechnica We then obtain: dZk0 j0 dt (t0) = Dk0Zk0 x̄x,j0 (t0) + F k0(Zj0(t0)) < 0, which means that Zk0 j0 (t) cannot grow to any higher value, and is still maximal over all other nodes. Such values, therefore, cannot reach Zk0(t1) > Bk0 . This lemma will be used in the following section to derive estimates of the grid function found by semi- discrete scheme (4), and then facilitate the proof of convergence for scheme (4). 4. Numerical analysis of method of lines Scheme (4), as a system of first order ODE’s with a convenient right-hand side, possesses the unique solution Z = Z(t) defined on (0, Tm) for a Tm > 0, following the Picard theorem. Due to Lemma 1 with the invariant region Σ independent of m, we see that all solutions of (4) are uniformly bounded in their values. The theory of ODE’s [23] then yields Tm = +∞. Solutions of (4) are then available on a selected time interval (0, T ) for all h. 4.1. A priori estimates We multiply (4) by its solution in terms of the scalar product (., .)h, use the integration by parts over ωh, and obtain: 1 2 d∥Z∥2 h dt + (DZx̄, Zx̄] = (F(Z), Z)h. Due to Lemma 1, the C1 map F is bounded with bounded derivative on Σ, i.e. |F(Z)| ≤ L|Z| for an L > 0. As D is a positive diagonal matrix with values bounded from below by d0 > 0, we have: 1 2 d∥Z∥2 h dt + d0||Zx̄]|2 ≤ L∥Z∥2 h. The Grönwall argument then yields: ∥Z∥2 h(t)+2d0 ∫ t 0 ||Zx̄]|2(t)dt ≤ ∥Z∥2 h(0)(1+e2Lt). (5) Multiplying (4) by dZ dt = Ż in terms of the scalar product (., .)h, and using the integration by parts over ωh yield: ∥Ż∥2 h(t) + 1 2 d(DZx̄, Zx̄] dt = ( F(Z), Ż ) h . Again, due to Lemma 1, |F(Z)| ≤ L|Z| for an L > 0. Using the Cauchy-Schwarz and Young inequalities, we have: ∥Ż∥2 h(t) + 1 2 d(DZx̄, Zx̄] dt ≤ L2 2 ∥Z∥2 h + 1 2∥Ż∥2 h, and: ∥Ż∥2 h(t) + d(DZx̄, Zx̄] dt ≤ L2∥Z∥2 h. The integration and the Grönwall argument then yield:∫ t 0 ∥Ż∥2 h(t)dt+ (DZx̄, Zx̄] (t) ≤ (DZx̄, Zx̄] (0) + L2∥Z∥2 h(0)e2Lt. (6) 4.2. Interpolation results As in [21, 22], we introduce the interpolation opera- tors mapping the grid functions defined on ωh to con- venient Lebesgue-integrable functions on (a, b). We define: • Qh : Hh → C([a, b]) such that for each u ∈ Hh (Qhu)(x) = uj−1 + ux̄,j−1(x− a− (j − 1)h), for x ∈ [a+ (j − 1)h, a+ jh]; • Sh : Hh → Lp((a, b)) for a given p ≥ 1 such that for each u ∈ Hh (Shu)(x) = uj , for x ∈ (a+ (j − 1 2 )h, a+ (j + 1 2 )h); • Ph : C([a, b]) → Hh such that for each u ∈ H (Phu)j = u(a+ jh), for j = 0, . . . ,m. Remark The operator Ph is linear and continu- ous from C([a, b]) to Hh, and can be extended to W (1) 2 ((a, b)) via density argument. Qhu is a continu- ous piecewise linear function, ∂x(Qhu) exists a.e. in (a, b). We proceed by determining basic properties of the above defined maps as proven in [21], taking V,W ∈ Hh: (1.) The scalar products in L2 and Hh are related as:∫ b a ShV ShWdx = (V,W )h. (7) (2.) The products of gradients are related as:∫ b a ∂x(QhV )∂x(QhW )dx = (Vx̄,Wx̄]. (8) (3.) The relation of norms is: ∥QhV ∥L2((a,b)) ≤ ∥ShV ∥L2((a,b)). (9) (4.) The difference of extrapolation operators is:∫ b a |QhV − ShV |2dx ≤ h2 6 ∥Vx̄]|2. (10) (5.) Let V ∈ C0,ν((a, b);Rd), ν ∈ (0, 1). Then: Sh(PhV ) → V in Ls((a, b);Rd), whenever h → 0, (11) for s > 1. 568 vol. 65 no. 5/2025 Method of lines for reaction-diffusion systems . . . (6.) Let V ∈ V ∩W (2) 2 ((a, b);Rd). Then: Qh(PhV ) → V (12) in V, if h → 0. Extending the results of Section 4.1 to the contin- uum of (a, b), we see that ∂xQh(PhUini) are bounded in L2((a, b)) (by Equation (12)), and Sh(PhUini) is bounded in L2((a, b)) (by Equation (7)) independently of h. Therefore, solution of Equation (4) obtained for Zini = PhUini satisfies: ∂xQhZ ∈ L∞(0, T ; L2((a, b))), ShZ ∈ L∞(0, T ; L2((a, b))), from which: ∂xQhZ ∈ L2(0, T ; L2((a, b))), ShZ ∈ L2(0, T ; L2((a, b))), are bounded independently of h. Moreover, we obtain that: ShŻ ∈ L2(0, T ; L2((a, b))), are bounded independently of h as follows from (6). We conclude that: QhZ ∈ L∞(0, T ; H1 0((a, b))), QhZ ∈ L2(0, T ; H1 0((a, b))), are bounded independently on h. According to (9), QhŻ ∈ L2(0, T ; L2((a, b))). 4.3. Passage to the limit Passing to a subsequence, we have: • QhnZ ⇀ U in L2(0, T ; V); • ShnŻ ⇀ ∂tU in L2(0, T ; V−1); • Shn Z ⇀ U in L2(0, T ; H). The non-linear terms in Equation (1) require a stronger convergence result. Using the lemma on the compact embedding (see e.g. [3, 24]), we conclude that QhnZ converges strongly in L2(0, T ; H). Denote their common limit as U and the weak limit of ShnŻ in L2(0, T ; H) as U (1). These limits exist as a conse- quence of the above-mentioned facts. The fact that U (1) = ∂tU is implied by the uniquenes of the limit in D′(0, T ), as:∫ T 0 (Shn Ż − Qhn Ż, q)ψ(t)dt = − ∫ T 0 (Shn Z − Qhn Z, q)ψ̇(t)dt, where q ∈ D((a, b)), ψ ∈ D(0, T ). Lemma 2. If U denotes the weak limit of Shn Z in L2(0, T ; H), and the strong limit of Qhn Z in L2(0, T ; H), then: F(Shn Z) → F(U) weakly in L2(0, T ; H). Proof. The argument is provided by the Lipschitz continuity of F within the invariant region. To pass to the limit, we proceed as in [3, 25, 26], and multiply (4) by the test functions PhnV , where V ∈ D((a, b)). We integrate it over ωh. Then, we have, in terms of H: (Shn Ż,Shn Phn V ) + (∂xQhn Z, ∂xQhn Phn V ) = (F(Shn Z),Shn Phn V ). (13) Knowing that (1.) Shn Ż converges weakly in L2(0, T ; H) to ∂tU ; (2.) ∂xQhn Z converges strongly in L2(0, T ; H) to ∂xU ; (3.) Shn Phn Uini converges strongly to Uini in H, we multiply (13) by a scalar function ψ(t) ∈ C1([0, T ]), for which ψ(T ) = 0, and integrate by parts. Taking into account all previous results, the fact that: Shn Z(0) = Shn Phn Uini, and the Lebesgue theorem, we are able to pass to the limit: (Uini, V )ψ(0) − ∫ T 0 (U, V )ψ̇dt = ∫ T 0 ψ(t) (−(∂xu, ∂xV ) + (F(U), V )) dt. (14) If ψ ∈ D(0, T ), we have: d(U, V ) dt + (∂xU, ∂xV ) = (F(U), V ), U |t=0 = Uini. (15) It remains to show that the weak solution satisfies the initial condition. Multiplying (15) by a scalar function ψ(t) ∈ C1([0, T ]), for which ψ(T ) = 0, and integrating by parts, we obtain: (U, V )|t=0ψ(0) − ∫ T 0 (U, V )ψ̇dt = ∫ T 0 ψ(t) (−(∂xu, ∂xV ) + (F(Z), V )) dt. (16) Subtracting (16) from (14), we get: (Uini − U |t=0, V )ψ(0) = 0, for all V ∈ D((a, b)). From this we see that U |t=0 = Uini in L2((a, b)). To prove uniqueness, consider two solutions of the prob- lem (2), denoted as Z and Z̄. Subtracting two systems of equations and denoting R = Z − Z̄, multiplying the system by R via the semi-discrete scheme (4), we have: 1 2 d∥R∥2 dt + (∂xR, ∂xR) = (F(Z) − F(Z̄), R), R(0) = 0. (17) As: |(F(Z) − F(Z̄), R)| ≤ L∥R∥2, Equation (17) can be integrated by the Grönwall ar- gument to see that R = 0. 569 Niels van der Meer, Michal Beneš Acta Polytechnica 4.4. Regularity Lemma 3. The solution U of (2) is in L2(0, T ;W (2) 2 ((a, b);Rd), its time derivative in L2(0, T ; H), and due to the embedding property, it is in C(0, T ;C1([a, b])). Proof. From the above facts, we see that ∂tU ∈ L2(0, T ; H). Then from (15), we obtain: (∂xU, ∂xV ) = (F(U), V ) − (dU dt , V ), which in terms of D((a, b)) yields: ∂xxU = F(U) − dU dt , and then ∂xxU ∈ L2(0, T ; H). The embedding then provides continuity. 4.5. Error estimates Consider the weak solution U of Equation (2), and the solution Z of the semi-discrete scheme Equation (4). We project U onto ωh within the weak equality Equa- tion (15), producing the approximation error Ψ: d(PhnU,PhnV ) dt + (Phn U x̄,Phn V x̄) = (F(Phn U),Phn V ) + Ψ, (18) Phn U |t=0 = Phn Uini. Then we subtract it from (4) to obtain:( d(Phn U − Z) dt ,Phn V ) + ((Phn U − Z)x̄,Phn V x̄) = (F(Phn U) − F(Z),Phn V ) + Ψ, PhnU − Z|t=0 = 0. Selecting Phn V = Phn U − Z, we get:( d(Phn U − Z) dt ,Phn U − Z ) + ((Phn U − Z)x̄,Phn (Phn U − Z)x̄) = (F(Phn U) − F(Z),Phn U − Z) + Ψ, PhnU − Z|t=0 = 0. from which the Grönwall argument yields: ∥PhnU − Z∥2(t) ≤ ∫ T 0 ∥Ψ(t)∥2eL(T −t)dt. As the difference and interpolation approximation errors imply ∥Ψ(t)∥2 ≈ O(h4), we obtain that ∥Phn U − Z∥(t) ≈ O(h2). We then conclude by stating the convergence prop- erty. Theorem 1. Let the above assumption hold and Uini ∈ V. Then the solution of scheme (4) converges to the weak solution of (2). Remark The analysis of the method using invari- ant regions can be extended, namely to the following cases. The spatial domain can have a higher dimen- sion than 1 and the diffusion terms will contain the Laplace operator. In that case, we would assume the domain to be bounded with the Lipschitz boundary (e.g. see [5, 27]). The positive definite matrix D in Equation (1) does not have to be diagonal, as de- scribed in [1]. However, the invariant region can then have a more general shape than prismatic. This would cover the processes including cross-diffusion (e.g. in the phase-field models [8], or the FitzHugh-Nagumo systems [28]). The Dirichlet boundary conditions are used in a variety of reaction-diffusion models to capture fixed solution values at the domain bound- ary (e.g. concentration, voltage, phase state – see e.g. [1, 3, 5, 10, 29]). As the Neumann boundary conditions also serve in such models, e.g. to express quantity conservation or reflexion, they can be pro- cessed by the described method as well, influencing the choice of space V = W (1) 2 ((a, b);Rd) and Hh. For details, see [11, 30, 31]. However, the convergence rate in Section 4.5 will be just O(h), provided the usual one-sided differences have been used to approximate the Neumann boundary conditions. More recently, the reaction-diffusion systems are considered on curves or surfaces (see e.g. [14, 32]). 5. Examples We accompany the numerical analysis of the method of lines by selecting two distinct examples of systems – the Brusselator reaction diffusion model from chem- istry introduced in [33] and discussed in [2, 34], and the FitzHugh-Nagumo model of excitable medium from biophysical context introduced in [35, 36], and recently discussed, for example, in [5, 28]. The method is well applicable to any other reaction-diffusion model satisfying underlying assumptions, e.g. to the Grey- Scott model (see [13]), the phase-field model [10], or the competition-diffusion systems such as [37]. 5.1. Numerical error measurement To study convergence of the numerical solution ob- tained by scheme (4), the numerical solution is com- puted using several grids with decreasing mesh size h, (increasing number of meshes m = b−a h ), and com- pared to the numerical solution computed on a very fine mesh with h̄, m̄ = b−a h̄ while projecting solutions on sparser meshes to the finest mesh using linear inter- polation. The solution on each mesh is stored at fixed time levels using the output time step τ̄ , NT = T τ̄ . Denoting Zk = Zk(t) the vector-valued grid func- tion representing the finest grid projection of the nu- merical solution computed on the grid with parameters hk, mk, and Zd = Zd(t) the numerical solution on the very fine mesh, we can express their distance measured 570 vol. 65 no. 5/2025 Method of lines for reaction-diffusion systems . . . in corresponding norm as: E2(hk) = max 0≤l≤NT  m̄∑ j=1 |Zk j (lτ̄) − Zd j (lτ̄)|2h̄ 1 2 . Convergence of the numerical solution is assessed us- ing the experimental order of convergence (EoC), as in [38], calculated from two numerical solutions ob- tained on grids with meshes h1, h2 as: EoC(h1, h2) = log E2(h1) E2(h2) . 5.2. Brusselator model As an example of a reaction diffusion system admit- ting invariant region, we recall the Brusselator system, which describes the following fictitious system of chem- ical reactions of two chemicals A and B, transforming to the chemicals D and E with the side products X and Y in an inert medium (see [2]): A k1−→ X, B + X k2−→ Y + D, 2X + Y k3−→ 3X, X k4−→ E. The reaction takes place in a reactor, characterized by the length L. The mentioned system of RDE has been proposed in [39], and widely studied by many authors. We refer, for example, to [34, 40, 41]. The Brusselator equations have the form: ∂u ∂t = Du L2 ∂2u ∂x2 +A− (B + 1)u+ u2v, ∂v ∂t = Dv L2 ∂2v ∂x2 +Bu− u2v, (19) where A, B, Du, Dv, and L are positive constants, t ∈ [0,+∞), x ∈ [0, 1]. The equations are endowed by the boundary conditions: u(t, 0) = A, u(t, 1) = A, v(t, 0) = B A , v(t, 1) = B A , (20) and the initial conditions: u(0, x) = u0(x), v(0, x) = v0(x). (21) The parameters A, B express the rescaled constant concentrations of the reactants in the reactor, Du, Dv are the diffusion coefficients of X and Y. The functions u(t, x), v(t, x) are the concentrations of X and Y rescaled with respect to A, B, and the reaction rates ki. We convert problem (19)–(21) to have homogeneous boundary conditions. Defining the transformation: X(t, x) = u(t, x) −A, Y (t, x) = v(t, x) − B A , (22) we obtain: ∂X ∂t = Du L2 ∂2X ∂x2 + (B − 1)X +A2Y + 2AXY + B A X2 +X2Y, ∂Y ∂t = Dv L2 ∂2Y ∂x2 −BX −A2Y − 2AXY − B A X2 −X2Y, (23) with homogeneous boundary conditions and with ini- tial conditions in the form: X|t=0 = u0 −A, Y |t=0 = v0 − B A . (24) Denoting: U(t) = ( X(t, ·) Y (t, ·) ) , F(U) = CU + B(U) + T(U), D = ( Du L2 , 0 0, Dv L2 ) , C = ( B−1, A2 −B, −A2 ) , B(U) = ( 2AXY + B AX 2 −2AXY − B AX 2 ) , T(U) = ( X2Y −X2Y ) , problem (23)–(24) can be written as: ∂U ∂t = D∆U + CU + B(U) + T(U), U |∂Ω = 0, (25) U |t=0 = Uini. The invariant regions for system (25) have been found by several authors, e.g. by Eslerová in [42] or by [43, 44]. Under such circumstances, convergence of numerical methods have been analyzed, for example, the nonlinear Galerkin method has been studied in [13, 15, 17]. For an example of numerical solution, we set the space and time intervals as (a, b) = (0, 1), (0, T ) = (0, 100), the reaction parameters A = 2.0, B = 5.45, the diffusion parameters Du = 0.008, Dv = 0.004, and L = 1.42. This choice is motivated by computational studies, for example, in [34, 40, 42]. More specifi- cally, the selected value of L generates an oscillatory behavior of the solution in time. The initial condition is set to: u0(x) = A+ sin(2πx), v0(x) = B A + sin(2πx). It means that it is a sinusoidal perturbation of the fixed point [A, B A ]. As known from [27, 34, 40, 42, 44], the solution approaches a periodic trajectory for t → +∞. We remark that the amplitude, the fre- quency of this perturbation, the spatial symmetry of the initial condition, or localization of its profile can influence the dynamics of the solution. Convergence with respect to a very fine solution obtained for m = 2 000 is summarized in terms of 571 Niels van der Meer, Michal Beneš Acta Polytechnica Mesh level i m E2(hi) EoC(hi−1, hi) 0 50 0.0688107 – 1 100 0.0165081 2.059 2 200 0.0040801 2.016 3 400 0.0009881 2.046 4 800 0.0002161 2.193 Table 1. Brusselator model – variable u: Table of numerical parameters and convergence errors. Mesh level i m E2(hi) EoC(hi−1, hi) 0 50 0.0789331 – 1 100 0.0188665 2.065 2 200 0.0046606 2.017 3 400 0.0011289 2.046 4 800 0.0002469 2.193 Table 2. Brusselator model – variable v: Table of numerical parameters and convergence errors. Figure 1. Brusselator dynamics. Time evolution of the u(t, 0.5)-component of the solution. errors and the experimental order of convergence in Tables 1 and 2. As expected, the convergence rate is close to 2 given by the order of approximation of the second derivative. The profile of the solution components u, and v is shown in Figures 1, and 2. Convergence in space profiles is depicted in Figures 3 and 4. 5.3. FitzHugh Nagumo model Some of mathematical models used in electrocardi- ology and electrophysiology, are systems of reaction- diffusion equations. One of them designed for the conduction of nerve impulses along an axon is the FitzHugh-Nagumo (FHN) model. First proposed in 1961 by FitzHugh [35], it is a simplification of the pioneering Hodgkin-Huxley model from 1956 (see the original work [45]). In 1962, Nagumo et al. [36] de- rived the equations from an active pulse transmission line simulating an animal nerve axon. Since then, many results regarding the qualitative behaviour of the FHN model have been published, for example, by Keener [46], and in [47, 48]. In this work, we consider the FitzHugh-Nagumo system in the form ∂tv = D∂2 xxv + f(v) − w + Iext, ∂tw = δD∂2 xxw + ϵ(βv − γw), (26) v|x=a = g1, v|x=b = g1, w|x=a = g2, w|x=b = g2, v|t=0 = vini, w|t=0 = wini, where D, δ, ϵ, β, γ > 0, and f is the reaction term, f(v) = v(1 − v)(v − α), α ∈ (0, 1). 572 vol. 65 no. 5/2025 Method of lines for reaction-diffusion systems . . . Figure 2. Brusselator dynamics. Time evolution of the v(t, 0.5)-component of the solution. Figure 3. Brusselator dynamics. Convergence of the space profile of the u-component of the solution for t = 100. The external excitation current Iext ∈ R can provide signal from neighbouring neurons. The model reflects important qualitative properties of excitable medium: • Excitability: small stimulus of a suitable form gen- erates a larger response of medium leading, after a while, to the initial state. • Refractoriness: if such a reaction occurs, the medium cannot be excited again for some time. Variable v represents a (normalized) membrane po- tential and is called the excitable. The variable w serves as a gating variable (see [46]). We also remark that system (26) can be formulated in several different forms in literature. The small parameter ϵ quantifies relative rates of excitation and recovery. The invariant regions for (26) are available and have been studied by several authors, for example, in [1] and recently in [5]. We, therefore, apply the approach described in this text and have the convergence of the numerical scheme guaranteed. A computational example has been performed on the space interval (0, 50) for the initial conditions: vini(x) = 1 2 exp ( − (x− 9)2 3 ) , wini(x) = 1 5 exp ( − (x− 7)2 3 ) , 573 Niels van der Meer, Michal Beneš Acta Polytechnica Figure 4. Brusselator dynamics. Convergence of the space profile of the v-component of the solution for t = 100. Mesh level i m E2(hi) EoC(hi−1, hi) 0 50 37.66 – 1 100 7.89 2.25 2 200 1.39 2.50 3 400 0.246 2.50 4 800 0.0433 2.50 5 1 600 0.0107 2.01 6 3 200 0.0035 1.61 7 6 400 0.00106 1.73 Table 3. FHN model – variable v: Table of numerical parameters and convergence errors. Mesh level i m E2(hi) EoC(hi−1, hi) 0 50 2.29 – 1 100 0.409 2.48 2 200 0.0728 2.49 3 400 0.0139 2.39 4 800 0.00368 1.92 5 1 600 0.00121 1.60 6 3 200 0.000393 1.62 7 6 400 0.000119 1.73 Table 4. FHN model – variable w: Table of numerical parameters and convergence errors. and for the boundary conditions: v(0, t) = 0, v(50, t) = 0, w(0, t) = 0, w(50, t) = 0. The remaining settings for this computation are: D = 0.1, α = 0.1, β = 0.3, γ = 1.0, δ = 0.001, ϵ = 0.01, Iext = 0.0. The time extent is (0, 200). The choice of the initial condition is motivated by modeling excitation in myocardium [5, 46]. The convergence with respect to a very fine solution obtained for m = 25 600 is summarized in terms of errors in Table 3, and in terms of the experimental order of convergence in Table 4. Convergence in space profiles is depicted in Figures 5 and 6. Again, as expected, the convergence rate is close to 2 given by the order of approximation of the second derivative. 6. Conclusion The finite-difference method of lines turns out to be ef- ficient, easy to implement, and accurate in approximat- ing nonlinear dynamics of reaction diffusion systems. Its convergence is shown by means of the generalized 574 vol. 65 no. 5/2025 Method of lines for reaction-diffusion systems . . . Figure 5. FHN dynamics. Convergence of the space profile of the v-component of the solution for t = 160. Figure 6. FHN dynamics. Convergence of the space profile of the w-component of the solution for t = 160. maximum principle known as the invariant-region concept. We have shown that the invariant-region property carries over to the method of lines, thereby facilitating the proofs. In two examples, we show how the method can be applied. 7. Dedication The authors devote the text to the memory of Prof. Ing. Miloslav Havlíček, DrSc. Acknowledgements The authors were partly supported by the project 25- 18265S of the Czech Science Foundation, and by the CTU under the grant No. SGS23/188/OHK4/3T/14. The second author is grateful to Prof. Ing. Miloslav Havlíček, DrSc. for his encouragement, interest in results of the research, kind and friendly attitude, and mediation of international contacts which led to many results at the scientific as well as educational level. References [1] J. Smoller. Shock waves and reaction – Diffusion equations. Springer, New York, 2nd edn., 1994. https://doi.org/10.1007/978-1-4612-0873-0 [2] M. Holodniok, M. Kubíček, M. Marek. Metody analýzy nelineárních dynamických modelů [In Czech; Methods of the analysis of nonlinear dynamical models]. Academia, Prague, 1986. 575 https://doi.org/10.1007/978-1-4612-0873-0 Niels van der Meer, Michal Beneš Acta Polytechnica [3] R. Temam. Infinite-dimensional dynamical systems in mechanics and physics. Springer, New York, 1988. ISBN 0-387-96638-2. [4] J. D. Murray. Mathematical biology. Springer, Berlin, Heidelberg, 1993. https://doi.org/10.1007/978-3-662-08542-4 [5] J. Kantner, M. Beneš. Mathematical model of signal propagation in excitable media. Discrete and Continuous Dynamical Systems – Series S 14(3):935– 951, 2021. https://doi.org/10.3934/dcdss.2020382 [6] F. Verhulst. Nonlinear differential equations and dynamical systems. Springer, Berlin, Heidelberg, 1996. https://doi.org/10.1007/978-3-642-61453-8 [7] T. Ohta, M. Mimura, R. Kobayashi. Higher- dimensional localized patterns in excitable media. Physica D: Nonlinear Phenomena 34(1):115–144, 1989. https://doi.org/10.1016/0167-2789(89)90230-3 [8] O. Penrose, P. C. Fife. On the relation between the standard phase-field model and a “thermodynamically consistent” phase-field model. Physica D: Nonlinear Phenomena 69(1):107–113, 1993. https://doi.org/10.1016/0167-2789(93)90183-2 [9] M. Beneš. Numerical solution of phase-field equations with a gradient coupling term. In W. Jäger, J. Nečas, O. John, et al. (eds.), Partial Differential Equations – Theory and Numerical Solution, pp. 25–33. New York, 2000. [10] M. Beneš. Mathematical analysis of phase-field equations with numerically efficient coupling terms. Interfaces and Free Boundaries 3(2):201–212, 2001. https://doi.org/10.4171/IFB/38 [11] M. Beneš, V. Chalupecký, K. Mikula. Geometrical image segmentation by the Allen-Cahn equation. Applied Numerical Mathematics 51(2–3):187–205, 2004. https://doi.org/10.1016/j.apnum.2004.05.001 [12] D. Hoff. Stability and convergence of finite difference methods for systems of nonlinear reaction-diffusion equations. SIAM Journal on Numerical Analysis 15(6):1161–1177, 1978. https://doi.org/10.1137/0715077 [13] J. Mach, M. Beneš, P. Strachota. Nonlinear Galerkin finite element method applied to the system of reaction- diffusion equations in one space dimension. Computers & Mathematics with Applications 73(9):2053–2065, 2017. https://doi.org/10.1016/j.camwa.2017.02.032 [14] M. Frittelli, A. Madzvamuse, I. Sgura, C. Venkataraman. Preserving invariance properties of reaction–diffusion systems on stationary surfaces. IMA Journal of Numerical Analysis 39(1):235–270, 2017. https://doi.org/10.1093/imanum/drx058 [15] J. Šembera, M. Beneš. Nonlinear Galerkin method for reaction–diffusion systems admitting invariant regions. Journal of Computational and Applied Mathematics 136(1–2):163–176, 2001. https://doi.org/10.1016/S0377-0427(00)00582-3 [16] M. Kolář. Computational studies of reaction-diffusion systems by nonlinear Galerkin method. American Journal of Computational Mathematics 3(2):137–146, 2013. https://doi.org/10.4236/ajcm.2013.32022 [17] J. Mach. Application of the nonlinear Galerkin FEM method to the numerical solution of a reaction-diffusion system in two dimensions. RIMS Kokyuroku Bessatsu B35:95–113, 2012. [18] M. Beneš, K. Mikula, T. Oberhuber, D. Ševčovič. Comparison study for level set and direct Lagrangian methods for computing Willmore flow of closed planar curves. Computing and Visualization in Science 12(6):307–317, 2009. https://doi.org/10.1007/s00791-008-0112-2 [19] P. Pauš. Numerical simulation of dislocation dynamics. In M. Vajsáblová, P. Struk (eds.), Proceedings of Slovak–Austrian Congress, Magia, pp. 45–52. Bratislava, 2007. ISBN 978-80-227-2796-9. [20] M. Kolář, D. Ševčovič. Evolution of multiple closed knotted curves in space. In P. Frolkovič, K. Mikula, D. Ševčovič (eds.), Proceedings of the Algoritmy 2024 Conference, pp. 129–138. Bratislava, 2024. [21] M. Beneš. Phase-field model of microstructure growth in solidification of pure substances. Ph.D. thesis, Czech Technical University in Prague, Faculty of Nuclear Sciences and Physical Engineering, 1997. [22] M. Beneš. Analysis of equations in the phase field model. In Z. Došlá, J. Kuben, J. Vosmanský (eds.), Proceedings of Equadiff 9, pp. 17–35. Brno, 1998. [2025-06-30]. https://dml.cz/handle/10338.dmlcz/700304 [23] L. S. Pontryagin. Ordinary differential equations. Addisson Wesley Publishing, Palo Alto, 1962. [24] T. Roubíček. A generalization of the Lions-Temam compact imbedding theorem. Časopis pro pěstování matematiky 115(4):338–342, 1990. [2025-06-30]. http://eudml.org/doc/21776 [25] J.-L. Lions. Quelques méthodes de résolution des problèmes aux limites non linéaires [In French; Several solution methods of nonlinear boundary-value problems]. Dunod, Gauthiers-Villars, Paris, 1969. [26] R. Temam. Navier-Stokes equations, theory and numerical analysis. North-Holland, Amsterdam, 1979. [27] Y. Ji, J. Shen, X. Mao. Pattern formation of Brusselator in the reaction-diffusion system. Discrete and Continuous Dynamical Systems – Series S 16(3–4):434– 459, 2023. https://doi.org/10.3934/dcdss.2022103 [28] G. Gambino, V. Giunta, M. C. Lombardo, G. Rubino. Cross-diffusion effects on stationary pattern formation in the FitzHugh-Nagumo model. Discrete and Continuous Dynamical Systems – Series B 27(12):7783– 7816, 2022. https://doi.org/10.3934/dcdsb.2022063 [29] V. Thomée, L. B. Wahlbin. On the existence of maximum principles in parabolic finite element equations. Mathematics of Computation 77(261):11–19, 2008. https://doi.org/10.1090/S0025-5718-07-02021-2 [30] A. A. Samarskii, V. B. Andreev. Raznostnyje metody dlja ellipticeskich uravnenij [In Russian; Difference methods for elliptic equations]. Nauka, Moscow, 1976. [31] A. A. Samarskii. Teorija raznostnych schem [In Russian; Theory of difference schemes]. Nauka, Moscow, 1977. 576 https://doi.org/10.1007/978-3-662-08542-4 https://doi.org/10.3934/dcdss.2020382 https://doi.org/10.1007/978-3-642-61453-8 https://doi.org/10.1016/0167-2789(89)90230-3 https://doi.org/10.1016/0167-2789(93)90183-2 https://doi.org/10.4171/IFB/38 https://doi.org/10.1016/j.apnum.2004.05.001 https://doi.org/10.1137/0715077 https://doi.org/10.1016/j.camwa.2017.02.032 https://doi.org/10.1093/imanum/drx058 https://doi.org/10.1016/S0377-0427(00)00582-3 https://doi.org/10.4236/ajcm.2013.32022 https://doi.org/10.1007/s00791-008-0112-2 https://dml.cz/handle/10338.dmlcz/700304 http://eudml.org/doc/21776 https://doi.org/10.3934/dcdss.2022103 https://doi.org/10.3934/dcdsb.2022063 https://doi.org/10.1090/S0025-5718-07-02021-2 vol. 65 no. 5/2025 Method of lines for reaction-diffusion systems . . . [32] D. Žurek. Electromechanical model of excitable medium. Master’s thesis, Czech Technical University in Prague, Faculty of Nuclear Sciences and Physical Engineering, Prague, 2024. [33] I. Prigogine, R. Lefever. Symmetry breaking instabilities in dissipative systems. II. The Journal of Chemical Physics 48(4):1695–1700, 1968. https://doi.org/10.1063/1.1668896 [34] M. Marek, I. Schreiber. Chaotic behaviour of deterministic dissipative systems, vol. 1. Cambridge University Press & Academia Prague, 1991. https://doi.org/10.1017/CBO9780511608162 [35] R. FitzHugh. Impulses and physiological states in theoretical models of nerve membrane. Biophysical Journal 1(6):445–466, 1961. https://doi.org/10.1016/S0006-3495(61)86902-6 [36] J. Nagumo, S. Arimoto, S. Yoshizawa. An active pulse transmission line simulating nerve axon. Proceedings of the IRE 50(10):2061–2070, 1962. https://doi.org/10.1109/JRPROC.1962.288235 [37] S.-I. Ei, R. Ikota, M. Mimura. Segregating partition problem in competition-diffusion systems. Interfaces and Free Boundaries 1(1):57–80, 1999. https://doi.org/10.4171/IFB/4 [38] G. Dziuk. Convergence of a semi-discrete scheme for the curve shortening flow. Mathematical Models and Methods in Applied Sciences 4(4):589–606, 1994. https://doi.org/10.1142/S0218202594000339 [39] G. Nicolis, I. Prigogine. Self-organisation in nonequilibrium systems. Wiley and Sons, New York, 1977. [40] P. Raschman, M. Kubíček, M. Marek. Concentration waves in reaction-diffusion systems. Scientific Papers of the Prague Institute of Chemical Technology K17:151–175, 1982. [41] M. Holodniok, P. Knedlík, M. Kubíček. Continuation of periodic solutions in parabolic partial differential equations. In T. Küpper, R. Seydel, H. Troger (eds.), Bifurcation: Analysis, Algorithms, Applications, pp. 122–130. Birkhäuser, Basel, 1987. https://doi.org/10.1007/978-3-0348-7241-6_13 [42] L. Čížková-Eslerová. Numerická analýza dynamiky reakčně-difúzních rovnic Brusselátor [In Czech; Numerical analysis of dynamics of Brusselator reaction-diffusion equations]. Master’s thesis, Czech Technical University in Prague, Faculty of Nuclear Sciences and Physical Engineering, Prague, 1996. https://doi.org/10.13140/RG.2.2.16452.37765 [43] S. Kouachi. Global existence for some strongly coupled reaction-diffusion systems non-dissipative via invariant regions techniques. Tech. Rep. hal-01429082, HAL, Open Science, 2017. [2025-06-30]. https://hal.archives-ouvertes.fr/hal-01429082 [44] S. M. Stoltz. Pattern formation in the Brusselator model of chemical reactions. Master’s thesis, University of Pretoria, Pretoria, 2016. [45] A. L. Hodgkin, A. F. Huxley, B. Katz. Measurement of current-voltage relations in the membrane of the giant axon of Loligo. The Journal of Physiology 116(4):424–448, 1952. https://doi.org/10.1113/jphysiol.1952.sp004716 [46] J. Keener, J. Sneyd. Mathematical physiology. Springer, New York, 1998. https://doi.org/10.1007/b98841 [47] B. Deng. The existence of infinitely many traveling front and back waves in the FitzHugh–Nagumo equations. SIAM Journal on Mathematical Analysis 22(6):1631– 1650, 1991. https://doi.org/10.1137/0522102 [48] J. Rauch, J. Smoller. Qualitative theory of the FitzHugh-Nagumo equations. Advances in Mathematics 27(1):12–24, 1978. https://doi.org/10.1016/0001-8708(78)90075-0 577 https://doi.org/10.1063/1.1668896 https://doi.org/10.1017/CBO9780511608162 https://doi.org/10.1016/S0006-3495(61)86902-6 https://doi.org/10.1109/JRPROC.1962.288235 https://doi.org/10.4171/IFB/4 https://doi.org/10.1142/S0218202594000339 https://doi.org/10.1007/978-3-0348-7241-6_13 https://doi.org/10.13140/RG.2.2.16452.37765 https://hal.archives-ouvertes.fr/hal-01429082 https://doi.org/10.1113/jphysiol.1952.sp004716 https://doi.org/10.1007/b98841 https://doi.org/10.1137/0522102 https://doi.org/10.1016/0001-8708(78)90075-0 Acta Polytechnica 65(5):566–577, 2025 1 Introduction 2 Reaction-diffusion systems 3 Method of lines 3.1 Finite-difference discretization 3.2 Invariant regions for the semi-discrete scheme 4 Numerical analysis of method of lines 4.1 A priori estimates 4.2 Interpolation results 4.3 Passage to the limit 4.4 Regularity 4.5 Error estimates 5 Examples 5.1 Numerical error measurement 5.2 Brusselator model 5.3 FitzHugh Nagumo model 6 Conclusion 7 Dedication Acknowledgements References