Keywords: Semilinear parabolic problem, Nonsmooth initial data, Galerkin finite element, Exponential Runge-Kutta

Mathematics Subject Classification (2010): 35K58, 65J08, 65M15, 65M60

1 Introduction↩︎

We consider the following semilinear parabolic equation \[\label{eq1} \left\{ \begin{align} &\frac{du(t)}{dt} = Au(t) + f(u(t)), \qquad 0<t\leq T, \\ &u(0) = u_0, \end{align}\right.\tag{1}\] in the Hilbert space \(X=L^2(\Omega)\), where \(\Omega\) is a convex polygonal domain in \(\mathbb{R}^d\) for \(d \in \{1,2,3\}\). The linear operator \(A:D(A)\subset X \rightarrow X\) is self-adjoint and negative definite, and it generates an analytic semigroup \(S(t)=e^{At}\) for \(t\geq 0\). The nonlinear term \(f:D(A^{\eta/2})\rightarrow X\) is a Nemytskii operator induced by a continuous function on \(\mathbb{R}\), where \(0 <\eta< 2\). The initial value satisfies \(u_0\in D(A^{\gamma/2})\) for \(0<\gamma<2\). Further details on problem 1 are given in Section 2.

This work aims to analyze the fully discrete error of this problem. Linear finite elements and high-order exponential Runge-Kutta (EERK) methods are employed for spatial and temporal discretizations, respectively, with an emphasis on the latter. Exponential integrators have proven to be highly effective for the time integration of parabolic equations. These integrators effectively mitigate the impact of stiffness by exactly handling the linear stiff terms within the scheme. With advancements in computational efficiency, significant efforts have been dedicated to the development of exponential integrators for semilinear parabolic problems. Various types of exponential integrators have been explored, including exponential Runge-Kutta methods [1]–[3], exponential multistep methods [4], [5], and exponential Rosenbrock methods [6], [7]. A comprehensive overview of exponential integrators is available in [8], [9].

Conventional error estimation techniques for numerical methods require boundedness assumptions on the derivatives of the solution and the nonlinear term in semilinear parabolic equations. These assumptions hold when the initial data and nonlinear terms are sufficiently smooth, and appropriate compatibility conditions are satisfied. However, if the initial data lacks smoothness, these boundedness conditions fail, leading to severe order reduction. To clearly illustrate this phenomenon, Table ¿tbl:tab:example? presents the temporal discretization errors and convergence orders for the Field-Noyes model, which serves as a representative semilinear parabolic system, under initial conditions of varying regularities (see [10]). Here, the parameter \(\gamma\) characterizes the exact regularity of the initial data, satisfying \(u_0 \in D(A^{\gamma/2})\).

max width=

The upper part of Table ¿tbl:tab:example? shows the numerical results using a second-order EERK method. The theoretical error analysis for this second-order scheme has already been rigorously established in our prior work [10], which proves that the convergence order strictly adapts to the initial regularity as \(\min(1+\gamma, 2)\). Conversely, the lower part of the table displays the numerical results when applying a third-order EERK method. Compared to the second-order case, the theoretical analysis for this higher-order scheme is significantly more complex, as it requires delicately bounding the higher-order Fréchet derivatives of the nonlinear terms under nonsmooth data. Developing a comprehensive analytical framework to resolve this complexity and establishing sharp error estimates for general high-order EERK methods is the primary objective of the present study.

In recent decades, extensive research has been conducted on the nonsmooth error analysis of abstract semilinear parabolic equations, covering both spatial and temporal discretization approaches. In the context of spatial discretization, Galerkin approximation techniques have been thoroughly analyzed, with comprehensive error estimates consolidated in the monograph [11]. For temporal discretization, various numerical schemes—including fully implicit, semi-implicit, exponential Rosenbrock, and implicit-explicit methods [12]–[16]—have been investigated, revealing the phenomenon of order reduction under nonsmooth initial conditions. However, these abstract analyses often yield suboptimal order estimates, primarily due to restrictive assumptions imposed on the nonlinear term. This gap motivates the necessity of rigorous nonsmooth error analysis tailored to specific PDEs.

For instance, the Navier-Stokes equations have been extensively studied, as summarized in [17]. Notable contributions include first-order convergence results for the Euler implicit/explicit scheme in [18], and suboptimal \(1.5\)-order convergence for a Crank-Nicolson/Adams-Bashforth scheme in [19], both requiring \(H^1\) initial regularity. Under weaker \(L^2\) initial conditions, [20] proved first-order convergence for a variable-stepsize semi-implicit method. Additionally, [21] extended the analysis to the Burgers equation, demonstrating \(1.5\)-order convergence under \(H^1\) initial conditions. These error estimates are derived via energy method techniques.

In our prior work [10], by applying alternative methodologies, we obtained sharp error estimates for a reaction-diffusion equation with initial values in \(D(A^{\gamma/2})\) (\(0<\gamma<2\)). As shown in the upper part of Table ¿tbl:tab:example?, when \(\gamma > 1\), the convergence order reaches the upper bound of \(2\). A natural question arises: whether employing a third-order EERK method can further improve the convergence order and whether it exhibits behavior similar to the second-order method when \(0<\gamma<1\). Experiments using a third-order method (detailed in 16 ) affirmatively answer both points (see the bottom of Table ¿tbl:tab:example?). Moreover, the analytical framework developed in [10] is sufficiently general to be extended to a broader class of nonlinear equations by representing the nonlinear terms via Nemytskii operators.

The present work develops a comprehensive framework for the error analysis of the general EERK method 15 , applied to a class of semilinear parabolic problems with nonlinearities, such as polynomial and rational functions, satisfying Assumptions 1, 2 and 4. The primary challenges, stemming from nonlinearities, are outlined in Remarks 2 and 3. The analysis ultimately establishes a temporal convergence order of \(\min(1+\gamma/2+\rho_1(\gamma)/2,\: p)\), where \(\rho_1(\gamma)\) characterizes the boundedness of \(f'(u(t))\).

The paper is organized as follows. Section 2 establishes the abstract framework for the class of semilinear parabolic problems and analyzes their well-posedness. Building on this foundation, Section 3 develops the spatial discretization scheme and establishes spatial error estimates. Section 4 formulates the EERK scheme and proves its stability, followed by a detailed Taylor expansion analysis to examine the structure of higher-order local error terms. Section 5 derives the key estimates concerning the Fréchet derivatives of the Nemytskii operator, ultimately leading to the sharp convergence rates. Section 6 presents numerical experiments that validate the theoretical findings. Finally, Section 7 provides concluding remarks and discusses future research directions.

2 Analysis framework↩︎

Let us begin by introducing standard notation. For \(s\geq 0\), we denote by \(\|\cdot\|_s\) the norm of the Sobolev spaces \(H^s=H^s(\Omega)\) over the domain \(\Omega\) (see, e.g., [22], [23]). When \(s=0\), \(H^0\) is equivalent to \(L^2=L^2(\Omega)\) with norm \(\|\cdot\|\) and inner product \((\cdot,\cdot)\). The \(L^\infty\) space consists of all bounded measurable functions on \(\Omega\). For \(s\geq 1\), the space \(H_N^s=H_N^s(\Omega)\) denotes the Sobolev space \(H^s\) subject to homogeneous Neumann boundary conditions. Given any real-valued function space \(Y\) and an interval \(I\subset \mathbb{R}\), we introduce the restricted space \(Y|_I=\{v\in Y \mid v(x)\in I \text{ a.e. on }\Omega \}\). Throughout this paper, we denote by \(C\) a generic positive constant and by \(\varepsilon\) a sufficiently small positive number, both of which may vary across instances.

With these notations established, we consider the semilinear parabolic problem 1 . The underlying space \(X\) of 1 is \(L^2\). Let \(\mathcal{L}(X)\) denote the Banach space of bounded linear operators from \(X\) to \(X\) with operator norm \(\Vert\cdot\Vert_{\mathcal{L}(X)}\). The linear operator \(-A\) is generated by the sesquilinear form \[a(u,v)=\sum_{i,j=1}^d\int_{\Omega}a_{ij}(x)D_iuD_j{v}\:\text{d}x+\int_{\Omega}c(x)u{v}\:\text{d}x,\quad u,v\in V\subset H^1,\] where \(a_{ij}(x)\in L^{\infty}\), \(c(x)\in L^{\infty}\). Moreover, \(a_{ij}(x)\) are symmetric and the following conditions hold: \[\begin{align} \sum_{i,j=1}^d a_{ij}(x)&\xi_i\xi_j\geq C|\xi|^2,&& \xi=(\xi_1,\ldots,\xi_n)\in\mathbb{R}^d,\text{ a.e. }x\in\Omega,&\\ c(x)&\geq c_0>0,&&\mathrm{~a.e.~}x\in\Omega.& \end{align}\] Then the operator \(-A\) is a positive definite self-adjoint sectorial operator on \(X\). According to [22], when \(V=H^1\), the domains of fractional powers of \(A\) are specified as \[\left\{\; \begin{align} &D(A^\theta)=H^{2\theta},\qquad \text{if}\:\: 0\leq\theta<\frac{3}{4},\\ &D(A^\theta)=H_N^{2\theta},\qquad \text{if}\:\:\frac{3}{4}<\theta\leq 1, \end{align}\right.\] with norm equivalence \[\begin{align} \label{NormEq1} C^{-1}\Vert u\Vert_{2\theta}\leq\Vert A^\theta u\Vert \leq C\Vert u\Vert _{2\theta}, \:\:u\in D(A^\theta). \end{align}\tag{2}\] While alternative choices of \(V\) would lead to operators with different boundary conditions (Dirichlet, Robin, etc.), we restrict our analysis to the \(H^1\)-case for simplicity of presentation.

Turning to the nonlinear term \(f(u)\), a key step identified in [10] for deriving fully discrete error estimates for a semilinear parabolic problem with quadratic nonlinearities was to establish the boundedness of its Fréchet derivatives in \(D(A^{s/2})\). By leveraging the specific structure of \(f(u)\), this was achieved using pointwise multiplication estimates in fractional Sobolev spaces. Building upon this methodology, the present study extends the analysis to general semilinear parabolic problems under the following smoothness assumption:

This assumption ensures that the Nemytskii operator \(f:L^\infty|_I\subset L^\infty \rightarrow L^\infty\) is Fréchet differentiable to an arbitrary order, and satisfies \[\label{Df95def} D_u^{(m)}f(u)(v_1,v_2,...,v_m) = f^{(m)}(u)\prod_{i=1}^m v_i, \quad m>0\text{, } u\in L^{\infty}|_I \text{, }v_i\in L^{\infty},\tag{3}\] where \(D_u^{(m)}f(u)\) is the \(m\)-th Fréchet derivative of \(f(u)\), \(f^{(m)}(\tau)\) is the \(m\)-th derivative of the scalar function \(f(\tau)\), and the right-hand side is defined by classical pointwise multiplication. Unless otherwise specified, the composite operator \(f(u)\) and the scalar function \(f(\tau)\) use the same symbol \(f\).

According to 3 , analyzing the boundedness of the Fréchet derivative of \(f(u)\) is equivalent to estimating the pointwise multiplication on its right-hand side. To this end, we first define the relation set \(\mathcal{B}\): \[\label{Set95B} \begin{align} \mathcal{B}=&\Big \{(s,s_1,s_2)\mid -\frac{3}{2}< s\leq s_i< \frac{3}{2}\:\: \text{for}\:\: i=1,2, \\ &\quad s_1+s_2> \frac{d}{2}+s ,\:\:\:s_1+s_2>0, \:\: s_1>0 \:\Big\}, \end{align}\tag{4}\] where the parameter bounds \((-3/2, 3/2)\) are imposed to facilitate the extension of Lemma 1 to its discrete version 26 . Then we introduce the binary pointwise multiplication estimate from [10].

The \(m\)-ary case in 3 can be generalized from the binary case. Additionally, the boundedness of the Nemytskii operator \(f^{(m)}(u)\) is unknown, hence the following assumption is made.

The reasonableness of this assumption will be illustrated with examples at the end of this section. Under the above lemma and assumptions, we can obtain the local Lipschitz condition for \(f(u)\). Selecting \(s\), \(s_1\) and \(s_2\) such that \((s,\rho_1(s_1),s_2)\in \mathcal{B}\) and expanding \(f(u)\) in a Taylor series at \(v\) gives \[\label{Lips1} \begin{align} \Vert A^{s/2} (f(u)-f(v))\Vert &= \left\Vert A^{s/2} \int_0^1 f'(v+\theta (u-v))(u-v) \text{d}\theta \right\Vert \\ &\leq C \int_{0}^{1} \left\Vert A^{\rho_1(s_1)/2} f'(v+\theta (u-v))\right\Vert \text{d}\theta \: \Vert A^{s_2/2} (u-v)\Vert \\ &\leq C \phi_1(\Vert A^{s_1/2}u\Vert+\Vert A^{s_1/2}v\Vert) \Vert A^{s_2/2} (u-v)\Vert, \end{align}\tag{5}\] where \(u,v \in L^\infty|_I \cap D(A^{s_1/2}) \cap D(A^{s_2/2})\). In particular, 5 implies the bound \[\label{Nf95bound} \Vert A^{s/2} f(u)\Vert\leq C \phi_1(\Vert A^{s_1/2}u\Vert) \Vert A^{s_2/2} u\Vert+C.\tag{6}\]

For a semilinear parabolic system with the operator \(A\) and a nonlinear term \(f(u)\) which satisfies the local Lipschitz condition 5 with \((0,\rho_1(\beta_1),\eta)\in \mathcal{B}\) (\(\eta>\beta_1\)), the local existence of solutions follows from [22]. Here, the parameter \(\beta_1\) determines the range of initial value, and we hope it is as small as possible. Therefore, we make the following assumptions about the initial value.

Thus, the semilinear parabolic problem 1 is rigorously established. By [22], the solution of 1 has the following regularity \[u\in \mathcal{C}((0,T];D(A))\cap \mathcal{C}([0,T];D(A^{\gamma/2}))\cap \mathcal{C}^1((0,T];X).\] It also holds that \[\begin{align} \label{bound95u} \|A^{s/2}u(t)\|\leq C&t^{\gamma/2-s/2}+C, \quad 0\leq s\leq 2, \end{align}\tag{7}\] where \(C\) is independent of \(t\). Using the variation-of-constants formula, we obtain \[\begin{align} u(t)=S(t) u_0+\int_{0}^{t}S(t-\tau)f(u(\tau))\text{d}\tau,\qquad t\in[0,T].\label{VoCF1} \end{align}\tag{8}\]

Next, we provide two examples of nonlinear terms to illustrate the reasonableness of the assumptions.

3 Error analysis for spatial discretization scheme↩︎

While the primary focus of this paper is on the temporal discretization error, a practical convergence analysis must also incorporate spatial discretization. In this section, we discretize problem 1 in space using a linear Galerkin finite element method and analyze the resulting error. The analysis of the temporal discretization for the resulting semi-discrete system 9 is presented in the subsequent section.

3.1 Spatial discretization scheme↩︎

Consider a regular triangulation \(\pi_h\) of the domain \(\Omega\) with maximum element diameter \(h\). We define \(V_h \subset V\) as the finite-dimensional subspace consisting of continuous, piecewise linear functions over \(\pi_h\), which serves as the approximation space \(X_h = V_h\) for our spatial discretization scheme. The projection operator \(P_h\) from \(X\) to \(X_h\) is defined as \[(P_h u,v_h)=(u,v_h),\quad \forall\: v_h\in X_h,\text{ for } u\in X.\] The discrete operator \(A_{h}:X_h\longrightarrow X_h\) is defined by \[(-A_{h}u_h,v_h)=a(u_h,v_h),\quad\forall\: v_h\in X_h,\text{ for }u_h\in X_h.\] Subsequently, we define the Ritz operator \(R_{h}:V\longrightarrow X_h\), \[{a}(R_{h}u,v_h)={a}(u,v_h),\qquad \forall\: v_h\in X_h,\text{ for }u\in V.\] Combining these components, we obtain the Galerkin finite element scheme for 9 : find \(u^h(t)\in X_h\), such that \[\label{eq2} \left\{\; \begin{align} &\frac{du^h(t)}{dt}=A_h u^h(t)+f_h(u^h(t)), \quad t\in (0,T],\\ &u^h(0)=P_h u_0, \end{align} \right.\tag{9}\] where \(f_h(u^h(t))=P_h (f(u^h(t)))\).

Our error analysis is based on the semigroup technique, exploiting the fact that the operators \(A\) and \(A_h\) generate the analytic semigroups \(S(t)=e^{tA}\) and \(S_h(t)=e^{tA_h}\) on \(X\) and \(X_h\), respectively. The necessary estimates for these semigroups are given in the following lemma.

Then, we establish the relationship between the original problem 1 and the spatial discretization problem 9 . The resulting estimates are crucial for analyzing the spatial discretization error, and will also be used later for studying the temporal error.

3.2 Spatial error analysis↩︎

In our previous work [10], we first investigate the error associated with the linear part of problem 9 , then analyze the spatial discretization error estimates for local solutions. Finally, we extend the solution to the global domain using this local error estimate. This section requires the same argument to establish our spatial discretization results. We present only the key differences and main results. First, the Lipschitz condition of \(f_h\) on the discrete space \(X_h\) follows directly from 5 , ?? and ?? : \[\begin{align} \label{Lips2} \Vert A_h^{s/2}P_h (f(u_h)-f(v_h))\Vert\leq C \phi_1(\Vert A_h^{\beta_1/2}u_h\Vert+\Vert A_h^{\beta_1/2}v_h\Vert) \Vert A_h^{s_2/2} (u_h-v_h)\Vert, \end{align}\tag{14}\] where \((s,\rho_1(\beta_1),s_2)\in \mathcal{B}\). In three dimensions, the parameter \(\beta_1\) can exceed \(1\), necessitating an extension of the parameter \(s\) in the following lemma.

The proof follows immediately from estimate ?? . Subsequently, based on the lemma above and the Lipschitz condition 14 , we can establish the well-posedness of the discrete problem and derive the spatial discretization error estimates in the following theorem.

4 Fully discrete scheme and Error representations↩︎

Building on the spatial discretization error analysis in Section 3, this section presents the fully discrete scheme and proceeds to examine the temporal error. Following the convergence analysis approach in [3], [27], we decompose the error into two components: stability and local error. While the stability of the scheme can be straightforwardly derived from the Lipschitz condition 14 , a rigorous analysis of the local error demands considerable effort, relying on a detailed Taylor expansion.

4.1 Fully discrete scheme↩︎

For the temporal discretization of semidiscrete problem 9 , we consider a class of \(\kappa\)-stage explicit exponential Runge-Kutta methods (EERK) with constant stepsize \(\delta=T/N\). \[\label{NumSche1} \begin{align} U^h_{ni}&=S_h(c_i \delta)u^h_n+\delta\sum_{j=1}^{i-1}a_{ij}(\delta A_h)f_h(U^h_{nj}),\quad 1\leq i\leq \kappa,\\ u^h_{n+1}&=S_h(\delta)u^h_n+\delta\sum^\kappa_{i=1} b_i(\delta A_h)f_h(U^h_{ni}). \end{align}\tag{15}\] The inner stages \(U^h_{ni}\) serve as approximations to \(u^h(t_n + c_i \delta)\), while the numerical solution \(u^h_{n+1}\) approximates the solution \(u^h(t)\) at time \(t_{n+1}\). The coefficients \(a_{ij}(\delta A_h)\) and \(b_i(\delta A_h)\) are chosen as linear combinations of the exponential functions \(\varphi_k(c_i\delta A_h)\) and \(\varphi_k(\delta A_h)\), respectively. These functions are given by \[\varphi_{0}(z)=e^{z},\quad\varphi_{k}(z)=\int_{0}^{1}e^{(1-\theta)z}\frac{\theta^{k-1}}{(k-1)!}\text{d}\theta,\quad k\geq 1,\] and thus satisfy the recurrence relation \[\varphi_{k+1}(z)=\frac{\varphi_k(z)-\varphi_k(0)}{z},\quad k\geq0.\] Similar to Lemma 2, for \(0\leq\alpha<1\), we get \[\Vert A_h^\alpha b_i(tA_h)\Vert_{\mathcal{L}(X_h)}+\Vert A_h^\alpha a_{ij}(tA_h)\Vert_{\mathcal{L}(X_h)} \leq Ct^{-\alpha},\quad t>0.\]

To keep the notation concise, we introduce the abbreviations \(a_{ij}=a_{ij}(\delta A_h)\), \(b_{i}=b_{i}(\delta A_h)\), \(\varphi_i=\varphi_{i}(\delta A_h)\), and \(\varphi_{ji}=\varphi_{j}(c_i\delta A_h)\). For a \(p\)th-order scheme with \(s\) stages, the coefficients can be compactly represented by a Butcher tableau. For higher accuracy, we consider the classical three-stage third-order exponential Runge-Kutta method (EERK3), originally proposed in [1]. The corresponding coefficients are detailed in 16 . While this scheme provides superior accuracy for smooth problems, our primary interest lies in its convergence behavior and order reduction phenomena when applied to equations with nonsmooth initial data. \[\label{EERK3} \renewcommand{\arraystretch}{1.2} \begin{array}{c|ccc} 0 & & & \\ 1/2 & \frac{1}{2}\varphi_{1,2} & & \\ 2/3 & \frac{2}{3}\varphi_{1,3} - \frac{8}{9}\varphi_{2,3} & \quad\frac{8}{9}\varphi_{2,3} & \\ \hline & \varphi_1 - \frac{3}{2}\varphi_2 & 0 & \quad\frac{3}{2}\varphi_2 \end{array}\tag{16}\]

Furthermore, we define \(g(t)=f_h(u^h(t))\), \(F(t)=A_h u^h(t)+g(t)\), and \(\tilde{u}^h_{n}=u^h(t_{n})\). In the single-step numerical scheme 15 , the local approximations \(\hat{u}^h_{n+1}\) and \(\hat{U}^h_{ni}\) are obtained after one step starting from \(\tilde{u}^h_n\). The temporal error is then split as follows: \[\begin{align} e_{n+1}=\tilde{u}^h_{n+1}-u^h_{n+1}=\tilde{u}^h_{n+1}-\hat{u}^h_{n+1}+\hat{u}^h_{n+1}-u^h_{n+1}=:\tilde{e}_{n+1}+\hat{e}_{n+1}, \end{align}\] where \(\tilde{e}_{n+1}\) and \(\hat{e}_{n+1}\) represent the local error and stability error of the numerical scheme, respectively.

4.2 Stability↩︎

By repeatedly applying ?? , there exist \(r_1\in (0,\frac{d}{2})\) such that \[\label{err95L2Spli} \begin{align} \Vert e_{n+1}\Vert &=\Vert S_h(\delta)e_n+\delta \mathcal{N}(e_n)+\tilde{e}_{n+1}\Vert \\ &=\left\Vert \sum_{k=1}^{n+1}S_h(t_{n+1-k})\tilde{e}_{k}+ \sum_{k=1}^{n}\delta S_h(t_{n-k})\mathcal{N}(e_k)\right\Vert \\ &\leq \left\Vert \sum_{k=1}^{n+1} S_h(t_{n+1-k})\tilde{e}_{k}\right\Vert + \delta\sum_{k=1}^{n-1} t_{n-k}^{-r_1/2} \Vert e_k \Vert + \delta^{1-r_1/2} \Vert e_n\Vert. \end{align}\tag{17}\] To estimate the first term on the right-hand side, it is necessary to describe the local error.

4.3 Taylor expansion of the local error↩︎

In this subsection, we apply the approach from [27] to expand the local error. Upon expansion, if the numerical method satisfies the corresponding order conditions, the local error will contain some higher-order remainder terms. The primary challenge in nonsmooth error analysis lies in accurately estimating these terms. We therefore analyze the expansion process to examine their composition. Like 3 , the nonlinear term \(f_h:X_h|_I\subset X_h\rightarrow X_h\) is Fréchet differentiable of arbitrary order, and satisfies \[\label{Dfh95def} D_u^{(m)}f_h(u_h)(v^h_1,v^h_2,...,v^h_m) = P_h(f^{(m)}(u_h)\prod_{i=1}^m v^h_i),\quad m>0\text{, }u_h, v_i^h\in X_h.\tag{18}\] Expressing the exact solution \(u^h(t)\) at time \(t_{n+1}\) by the variation-of-constants formula, \[\tilde{u}^h_{n+1}=S_h(\delta)\tilde{u}^h_n+\delta\int_{0}^{1}S_h((1-\theta)\delta)g(t_n+\theta \delta)\text{d}\theta,\] then expanding \(g(t_n +\theta \delta)\) in a Taylor series at \(t_n\) gives \[\label{Expan95exa} \begin{align} \tilde{u}^h_{n+1}=&\:\tilde{u}^h_n+\delta \varphi_1F(t_n)+\sum_{i=2}^{p}\delta^i \varphi_i g^{(i-1)}(t_n)\\ &+\delta^{p+1} \int_{0}^{1}S_h((1-\theta)\delta) \frac{\theta^p}{(p-1)!} \int_{0}^{1}(1-s)^{p-1} g^{(p)}(t_n+s \theta \delta)\text{d}s \text{d}\theta. \end{align}\tag{19}\]

After the expansion of the exact solution \(\tilde{u}^h_{n+1}\), we turn to expand the local numerical solution \(\hat{u}^h_{n+1}\). It is desired that the numerical scheme preserves the equilibrium. This preservation property can be guaranteed by requiring \[\label{orderCond951} \sum_{j=1}^{i-1}a_{ij}=c_i\varphi_1,\:\: i=1,\ldots,\kappa,\quad\sum_{i=1}^\kappa b_i=\varphi_1.\tag{20}\] Taking into account the above conditions, we reformulate the EERK method 15 and consider one step with initial value \(\widetilde{u}^h_n\). \[\label{NumSche2} \begin{align} \widehat{U}^h_{ni}&=\tilde{u}^h_{n}+c_{i}\delta\varphi_{1i}F(t_n)+\delta\sum_{j=2}^{i-1}a_{ij}\widehat{D}_{nj},\\ \hat{u}^h_{n+1}&=\tilde{u}^h_{n}+\delta\varphi_{1}F(t_n)+\delta\sum_{i=2}^{\kappa}b_{i}\widehat{D}_{ni}, \end{align}\tag{21}\] with \(\widehat{D}_{ni}=f_h(\widehat{U}^h_{ni})-f_h(\tilde{u}^h_n).\) Expanding \(\widehat{D}_{ni}\) in Taylor series at \(\widetilde{u}^h_n\), we obtain \[\label{recur195mid} \begin{align} \widehat{D}_{ni} &=\delta \int_0^1 P_h( f'(\tilde{u}^h_n+\theta \delta V_i)V_i) \text{d}\theta \\ &= \sum_{k=1}^{m}\delta^k\frac{1}{k!}P_h(f^{(k)}(\tilde{u}^h_n)(V_i)^k)+ \delta^{m+1}\int_0^1 \frac{(1-\theta)^m}{m!}P_h(f^{(m+1)}(\tilde{u}^h_n+\theta \delta V_i)(V_i)^{m+1}) \text{d}\theta, \end{align}\tag{22}\] with \[\label{recur195mid2} V_{i}=\frac{1}{\delta}\left(\widehat{U}^h_{ni}-\tilde{u}^h_{n}\right)=c_{i}\varphi_{1}(c_{i}\delta A_h)F(t_n)+\sum_{j=2}^{i-1}a_{ij}\widehat{D}_{nj}.\tag{23}\] The first term of right hand of above equality can be expanded as follow. \[\label{recur195end} \begin{align} \varphi_1(c_i \delta A_h)F(t_n) &= \varphi_1(c_i \delta A_h) (\tilde{u}^h_n)'\\ &= (\tilde{u}^h_n)' + c_i \delta \varphi_2(c_i \delta A_h) ( (\tilde{u}^h_n)''-g'(t_n)) \\ &= (\tilde{u}^h_n)'+\frac{c_i \delta}{2!}\mathbf{X}_i + c_i^2 \delta^2 \varphi_3(c_i \delta A_h)((\tilde{u}^h_n)^{(3)}-g''(t_n)), \end{align}\tag{24}\] where \(\mathbf{X}_i=(\tilde{u}^h_n)^{\prime\prime}-2!{\varphi}_2(c_i\delta A_h)g'(t_n)\) and the equality can be expanded infinitely. A finite number of nested iterations are performed among 22 , 23 , and 24 (the number of iterations depends on the stage \(\kappa\) of scheme 21 ), and then substituted into 21 to obtain the expansion of the numerical solution. Subtracting 19 yields the local error, from which the order conditions of the method can be established (see, e.g., [3], [27]). Below is an example of the expansion with a fifth-order remainder term, where “order” denotes the power of the step size \(\delta\). \[\begin{align} \tilde{e}_{n+1} &=\delta^2 \psi_2 g'(t_n)+\delta^3\psi_3 g''(t_n) +\delta^4\psi_4 g^{(3)}(t_n)+ \boldsymbol{R}_4(t_n) + \boldsymbol{O}_{5}(t_n) \end{align}\] with \[\begin{align} \boldsymbol{R}_4(t_n) &=\delta^{3}\sum_{i=2}^{\kappa}b_{i}P_h(f^{\prime}(\tilde{u}^h_{n})\psi_{2,i}g'(t_n))+\delta^{4}\sum_{i=2}^{\kappa}b_{i}P_h(f^{\prime}(\tilde{u}^h_{n})\psi_{3,i}g''(t_n))\\ &\qquad +\delta^{4}\sum_{i=2}^{\kappa}b_{i}P_h(f^{\prime}(\tilde{u}^h_{n})\sum_{j=2}^{i-1}a_{ij}P_h(f^{\prime}(\tilde{u}^h_{n})\psi_{2,j}g'(t_n)))+\delta^{4}\sum_{i=2}^{\kappa}b_{i}c_{i}P_h(f^{\prime\prime}(\tilde{u}^h_{n})(\tilde{u}^h_{n})^{\prime}\psi_{2,i}g'(t_n)), \end{align}\] where \(\psi_j=\sum_{i=2}^\kappa b_i\frac{c_i^{j-1}}{(j-1)!}-\varphi_j\), \(\psi_{ji}=\sum_{k=2}^{i-1}a_{ik}\frac{c_k^{j-1}}{(j-1)!}-c_i^j\varphi_{ji}\).

For a \(p\)th order EERK method, its local error expansion is given by: \[\label{Err95localReprens} \begin{align} \tilde{e}_{k+1}&=\delta^p\big(\psi_p(\delta A_h)-\psi_p(0)\big)g^{(p)}(t_k)+\delta^p\sum_{i=2}^\kappa\big(b_i(\delta A_h)-b_i(0)\big)\mathbf{Q}_p(t_{k})+\mathbf{O}_{p+1}(t_{k})\\ &=:\tilde{e}_{k+1}^{(p)}+\mathbf{O}_{p+1}(t_k). \end{align}\tag{25}\] Here \(\boldsymbol{Q}_p(t_k)\) denotes the terms multiplying \(\delta^p\sum_{i=2}^{\kappa}\) in expansion and \(\mathbf{O}_{p+1}(t_k)\) is the higher-order remainder.

The expansion 25 leaves some higher-order terms, which satisfy the following pattern:

1. Based on the expected order of the method, we can control the degree of expansion in 19 , 22 , and 24 ; for instance, in the error expansion of a \(p\)th-order method, the highest derivatives involved are \(f^{(p)}(u_h)\) and \(D_t^p u^h(t)\).

2. Since the estimates for \(S_h(\delta)\) and \(\varphi_i(\delta A_h)\) are identical under the \(\Vert A_h^{s/2}\cdot\Vert\) norm \((-2<s<2)\), we can equivalently substitute the latter with the former. The same applies to the exponential functions generated by \(\varphi_i\).

3. The boundedness of \(f^{(m+1)}(\tilde{u}^h_n+\theta \delta V_i)\) is identical to that of \(f^{(m+1)}(\tilde{u}^h_n)\), and the integral of the remaining terms in 22 is bounded; thus, the former can be replaced by the latter.

Below is an example of an equivalent substitution. \[\begin{align} &\delta^{4}\sum_{i=2}^{\kappa}b_{i}c_{i}\int_0^1 (1-\theta)P_h(f^{(2)}(\tilde{u}^h_n+\theta \delta V_i)(\tilde{u}^h_{n})^{\prime}\psi_{2,i}g'(t_n))\text{d}\theta \\ \sim\; & \delta^{4}S_h(\delta)P_h(f^{(2)}(\tilde{u}^h_{n})(\tilde{u}^h_{n})^{\prime}S_h(\delta)g'(t_n)) \end{align}\] It remains to estimate the derivative term on the right-hand side. Clearly, the orders of \(f^{(i)}(\tilde{u}^h)\) and \((\tilde{u}^h_{n})^{(j)}\) exhibit a specific pattern determined by the powers of \(\delta\). To characterize this pattern, we introduce a family of function sets, \(\mathcal{K}_m^p\), where the index \(m\) indicates the power of \(\delta^m\) that pairs with the corresponding derivative. Denote \(\mathcal{K}^{p}_1=\{D_t u^h(t)\}\) and \(\mathcal{K}'=\{ D_t^i u^h(t) \:|\: i\geq 1\}\). The set \(\mathcal{K}_{m+1}^p\) are constructed in a recursive way. For \(1\leq m\leq p-1\), \[\begin{align} \mathcal{K}^{{p}}_{m+1}&=\{D_t^{m+1} u^h(t)\}\cup \{ S_h(\delta) P_h (f^{(j)}(u^h(t))\Pi_{i=1}^j v^h_i(t)) \big|\\ &\qquad 1\leq j\leq {p},\: v^h_{i}(t)\in \mathcal{K}_{k_i}^{p},\:\sum_{i=1}^{j}k_i=m\}. \end{align}\] For \(m\geq {p}\), \[\begin{align} \mathcal{K}^{p}_{m+1}= \{ S_h(\delta) P_h(f^{(j)}(u^h(t))\Pi_{i=1}^j v^h_i(t)) \big| \:1\leq j\leq p,\: v^h_{i}(t)\in \mathcal{K}_{k_i}^{p},\:\sum_{i=1}^{j}k_i=m \}. \end{align}\] Define \(\text{span}(\mathcal{K}_m^p)\) as the vector space generated by the function set \(\mathcal{K}_m^p\). Then there exist \(V_i(t)\in \text{span}(\mathcal{K}_i^p)\) for \(i\geq p\), such that \(\boldsymbol{Q}_p(t_k)\) and \(\mathbf{O}_{p+1}(t_k)\) can be equivalently represented by \(V_p(t_k)\) and \(\delta^{p+1}V_{p+1}(t_k)+\delta^{p+2}V_{p+2}(t_k)+...\), respectively.

5 Error estimate↩︎

In the previous section, the time discretization error was separated into two parts, resulting in 17 . It remains to estimate the local error 25 , which will be referred to as a \(p\)th-order expansion. For higher-order expansions, the analysis is complicated by the involvement of higher-order Fréchet derivatives of \(f_h\). Therefore, we restrict our presentation to the second- and third-order cases and set \(\hat{p}=2\) or \(3\). The procedure for higher orders is analogous. Note that higher-order methods inherently satisfy lower-order conditions, thus the corresponding expansions remain valid, and the final error estimates are still applicable. For a \(\hat{p}\)th-order expansion, the terms \(\boldsymbol{Q}_{\hat{p}}(t_k)\) and \(\boldsymbol{O}_{\hat{p}+1}(t_k)\) can be equivalently represented by a linear combination of the elements in \(\mathcal{K}^{\hat{p}}_{m}\). To estimate the equivalent terms in \(\mathcal{K}_m^{\hat{p}}\), it is necessary to analyze the boundedness of the Fréchet derivative of \(f_h\). Following this derivation, the result is simplified to Assumption 4. Then the estimates concerning the equivalent terms are provided in Proposition 3. Substituting these into 17 yields the temporal error in Theorem 2.

5.1 Boundedness of the Fréchet derivatives of \(f_h(u^h(t))\)↩︎

Combining Lemma 1 with ?? and ?? yields its discrete counterpart. \[\begin{align} \label{mult95Xh} \Vert A_h^{s/2}P_h(v_1 v^h_2) \Vert\leq C\Vert A^{s_1/2} v_1\Vert \Vert A_h^{s_2/2}v^h_2\Vert, \end{align}\tag{26}\] where \((s,s_1,s_2)\in \mathcal{B}\). Note that if \(u\) and \(v\) satisfy homogeneous Neumann conditions, so does their product \(uv\). Moreover, when \(s>\frac{3}{2}\), the fractional Sobolev space is an algebra. Consequently, the parameter range \(\mathcal{B}\) in Lemma 1 can be extended. For the discrete version \(P_h(f'(u^h(t))v^h)\), the numerical solution \(u^h(t)\) lacks sufficient regularity for a direct estimation of \(f'(u^h(t))\) in the norm \(\Vert A^{s_1/2}\cdot \Vert\). To overcome this, one may substitute \(f'(u^h(t))\) with \(f'(u(t))\), and the resulting error can be controlled by applying the inverse inequality ?? .

Similarly, we can extend the above binary pointwise multiplication estimates to the m-ary case, as shown below. \[\begin{align} \Vert A_h^{(r-\varepsilon)/2}P_h(v_1 v^h_2 ... v_m^h) \Vert & \leq C(\Vert A^{r/2} w_1\Vert+1) \Vert A_h^{r/2}v^h_2\Vert ... \Vert A_h^{r/2}v^h_m\Vert. \tag{27} \\ \Vert A_h^{-r/2}P_h(v_1 v^h_2 ... v_m^h) \Vert \leq C(\Vert &A^{r/2} w_1\Vert+1) \Vert A_h^{r/2}v^h_2\Vert ... \Vert A_h^{r/2}v^h_{m-1}\Vert \Vert A_h^{(-r+\varepsilon)/2}v^h_m\Vert.\tag{28} \end{align}\] where \(v_i^h\in X_h\) for \(i=2,3,..,m.\)

We now proceed to estimate the Fréchet derivatives of \(f_h(u^h)\). Define a shift for \(s\): \[\label{eq:sigma} \sigma(s) = \max(s+d/2-\rho_1(\gamma),-\rho_1(\gamma),s)+\varepsilon, \quad -2\leq s<\rho_1(\gamma).\tag{29}\] Here, the inequality \(\sigma(s) - s < 2\) holds, which will be frequently employed in subsequent arguments. We want to show that \[\begin{align} \label{frech95first} \Vert A_h^{s/2}P_h (f'(u^h(t))v^h) \Vert\leq C \Vert A_h^{\sigma(s)/2}v^h \Vert, \qquad -2\leq s<\rho_1(\gamma). \end{align}\tag{30}\] For the case \(|s|<{3}/{2}\) and \(\gamma<3/2\) where \((s,\rho_1(\gamma),\sigma(s))\in \mathcal{B}\), estimate 30 follows directly from ?? and 26 . Note that through norm embedding, estimate 26 remains valid for \((s-c_1, \rho_1(\gamma)+c_2, \sigma(s)+c_3)\) even outside the set \(\mathcal{B}\) for any \(c_1,c_2,c_3>0\). We will omit this special case in what follows. When \(|s|>3/2\) and \(\gamma>3/2\), we have \(\rho_1(\gamma) = \rho_2(\gamma) = \gamma>3/2\). Let \(w_1=f'(u(t))\) and apply Lemma 1 with \((\varepsilon,d/2-\varepsilon,\varepsilon)\in \mathcal{B}\), combined with ?? , to derive \[\Vert f'(u^h(t))-f'(u(t))\Vert_\varepsilon\leq C \Vert A^{(d/2-\varepsilon)/2}f''(u^h(t)+u(t))\Vert \Vert u(t)-u^h(t)\Vert_{\varepsilon} \leq C h^{\gamma-\varepsilon}.\] This, together with the boundedness of \(\Vert A^{\rho_1(\gamma)/2}f'(u(t))\Vert\) and Lemma 5, yields 30 . It implies that \[\label{bound95g} \Vert A_h^{s/2} f_h(u^h(t))\Vert\leq \Vert A_h^{s/2} P_h(f'(u^h(t))u^h(t))\Vert+C\leq C(t^{\gamma/2-\sigma(s)/2}+1).\tag{31}\]

For the estimation of the second- and third-order Fréchet derivatives, involving the parameters \(\rho_2(\gamma)\) and \(\rho_3(\gamma)\), the derivation is quite tedious. To simplify the presentation, we extract some common features and make the following assumptions.

Dominance of the first-order derivative is guaranteed by this assumption. Specifically, when estimate ?? in Proposition 2 holds, we can consequently derive the validity of ?? . In what follows, illustrative examples are provided to validate the assumption.

5.2 Estimation for \(\mathcal{K}_m^{\hat{p}}\)↩︎

Having validated the Assumption 4 through some examples, we now turn to estimate the equivalent terms in \(\mathcal{K}_m^{\hat{p}}\). We begin by analyzing the derivatives of \(u^h(t)\).

- *$m=2:$
  $\|A_h^{s/2}v^h(t)\| \leq C(t^{-1+\gamma/2-\sigma(s)/2} + 1)$ for
  $-2\leq s<\rho_1(\gamma)$.*

- *$m\geq 3:$
  $\delta^m\|A_h^{s/2}v^h(t)\| \leq C(\delta^3 t^{-2+\gamma/2-\sigma(s)/2} + 1)$
  for $s\geq -2$ with $\sigma(s)<\rho_1(\gamma)$.*

 - *$m=2,\:3:$
   $\|A_h^{s/2}v^h(t)\| \leq C(t^{-m+1+\gamma/2-\sigma(s)/2} + 1)$
   for $s\geq -2$ with $\sigma(s)<\rho_1(\gamma)$.*

 - *$m\geq 4:$
   $\delta^m\|A_h^{s/2}v^h(t)\| \leq C(\delta^4 t^{-3+\gamma/2-\sigma(s)/2} + 1)$
   for $s\geq -2$ with $\sigma(\sigma(s))<\rho_1(\gamma)$.*

5.3 Main result↩︎

Now, the main result of this paper is formulated in the following theorem.

6 Numerical Experiments↩︎

In this section, we present numerical tests to support the theoretical analysis in Section 5. We consider the semilinear parabolic problem 1 in the domain \(\Omega = (0,1)\times (0,1)\) up to \(T = 1\). The sectorial operator \(A\) is realization of \(\Delta-I\) in \(L^2\) under the homogeneous Neumann boundary condition on \(\partial \Omega\). The spatial discretization is performed using the linear Galerkin finite element method 9 with an adequately small spatial mesh \(h=1/64\), while the temporal discretization relies on the third-order EERK method 16 .

The main focus is to compute the convergence order of the third-order EERK method at \(T\) for semilinear parabolic problems 9 with various initial values and nonlinear terms. We consider the following initial values.

(i) \(u_0(x_1,x_2)=0.5\:\text{sign}(x_2-0.5)+1.3\in D(A^{1/4-\varepsilon})\cap L^\infty\).

(ii) \(u_0(x_1,x_2)=(x_1^2+x_2^2)^{-1/4}\in D(A^{1/4-\varepsilon})\).

(iii) \(u_0(x_1,x_2)=0.5(x_1^2+x_2^2)+1\in D(A^{3/4-\varepsilon})\).

(iv) \(u_0(x_1,x_2) = (2x_1^{3/2}-x_1^3)(2x_2^{3/2}-x_2^3) +1 \in D(A^{1-\varepsilon})\).

Substitute these initial values into the semilinear parabolic problem with the following nonlinear terms.

(1) \(f(u)=-(u+1)(u-1.5)+u\) with \(\rho_1(\gamma)=\gamma\).

(2) \(f(u)= -1/8\:(u+1)^2(u-1.5)+u\) with \(\rho_1(\gamma)=2\gamma-1\).

(3) \(f(u)=-u^4/(1+u^2)+u+81/52\) with \(\rho_1(\gamma)=\gamma\).

It can be verified through Examples 1-4 that Assumptions 1-4 hold. The numerical orders of convergence are computed by \[\log\left(\frac{\Vert U_{T,N}-U_{T,2N} \Vert }{\Vert U_{T,2N}-U_{T,4N} \Vert}\right)/\log(2),\] where \(U_{T,N}\) represents the value at \(T\), obtained by applying the third-order EERK method to the discretized problem 9 with a constant step size of \(T/N\). For \(\gamma<1\), Theorem 2 yields a convergence order of \(1+{\gamma}/{2}+\rho_1(\gamma)/{2}-\varepsilon\). The numerical results for initial values () and () in Tables 1-3 align with the theoretical predictions. Notably, initial value () belongs to \(L^{\infty}\), while for the nonlinear term in case \((2)\), the solution \(\Vert u^h(t)\Vert_{\infty}\) remains uniformly bounded with respect to the \(h\). As shown in Example 1, we have \(\rho_1(1/2)=1/2\) in this case. For \(\gamma>1\), Theorem 2 yields a convergence order of \(1+\gamma-\varepsilon\). Numerical results for initial values () and () in Tables 1-3 validate this theoretical analysis.

7 Conclusions↩︎

In this paper, we have developed a comprehensive and rigorous numerical analysis framework for a class of semilinear parabolic problems subjected to nonsmooth initial data. By employing a linear Galerkin finite element method in space and a high-order explicit exponential Runge-Kutta (EERK) method in time, we systematically investigated the severe order reduction phenomenon inherent in nonsmooth settings. The primary mathematical difficulty of bounding the higher-order Fréchet derivatives of the Nemytskii operator was rigorously resolved through a combination of fractional power space techniques and analytic semigroup estimates. Ultimately, our theoretical analysis establishes a sharp temporal convergence rate of \(\min(1 + \gamma/2 + \rho_1(\gamma)/2, \:p)\), which strictly adapts to the exact regularity of the initial data, thereby significantly improving upon suboptimal estimates frequently encountered in the existing literature.

The theoretical framework established herein not only elucidates the intricate interplay between strong nonlinearity and strong stiffness but also provides a robust and extensible foundation for future investigations. Natural continuations of this work include adapting the current analysis to phase-field models governed by the Allen-Cahn and Cahn-Hilliard equations, where handling non-linear stiffness and phase separation under low-regularity conditions remains profoundly challenging. Furthermore, investigating the convergence behavior of fully implicit or exponential Rosenbrock schemes within this fractional Sobolev framework constitutes another important direction for our forthcoming research.

Funding↩︎

This work was supported by the Guangdong Basic and Applied Basic Research Foundation of China under Grand No.2026A1515012144.

Author Contributions↩︎

S. Yang performed conceptualization, methodology, software development and wrote the original draft. R. Zhang carried out validation, formal analysis and revised the manuscript. Z. Yu completed formal analysis and validation. J. Fang contributed to conceptualization, methodology and revised the manuscript. All authors reviewed the manuscript.

Declarations↩︎

Data Availability

No datasets were generated or analysed during the study.

Conflicts of Interest

The authors declare that they have no conflict of interest.

Competing interests

The authors declare no competing interests.

References↩︎