ETDRK4–Chebyshev Collocation for the Generalized Burgers–Huxley Equation: Machine-Precision Benchmarks and a Corrected Exact Solution

1. Introduction

The propagation of excitation waves in one-dimensional excitable and bistable media — action potentials in nerve fibres, invasion fronts in ecological populations, and nonlinear dissipative waves in gases and dispersive optical media — is, in its simplest continuum form, described by a scalar reaction–diffusion–advection equation combining a Burgers-type convective nonlinearity with a Huxley-type cubic (or polynomial) reaction term [18,19,20,21,22]. The generalized Burgers–Huxley (gBH) equation

$$u_t+\alpha u^{\delta }{u}_x-u_{xx}=\beta u(1-u^{\delta })(u^{\delta }-\gamma ),(x,t)\in \Omega \times (0,T],$$

(1)

with $\alpha ,\beta \ge 0$, $\delta >0$, and typically $\gamma \in (0,1)$, is the canonical one-dimensional model unifying these three classes of phenomena, interpolating between the Burgers ($\beta =0$), Huxley ($\alpha =0$) and Newell–Whitehead ($\alpha =0,\gamma =-1,\delta =1$) equations [1,2]. When $\alpha =0$, equation (1) reduces to the Huxley component, which is the simplest continuum model of the Hodgkin–Huxley action-potential dynamics developed to explain nerve-impulse propagation in the squid giant axon [18], and is closely related to the FitzHugh–Nagumo excitable-media model used in computational neuroscience and cardiac electrophysiology [19,20]. The reaction term $\beta u(1-u^{\delta })(u^{\delta }-\gamma )$ is the bistable (Allee-threshold) kinetics that also arises in population-dynamics models with a threshold density below which a population collapses and above which it grows to carrying capacity [21], while the advection term $\alpha u^{\delta }{u}_x$ introduces the Burgers-type nonlinearity governing shock formation in fluid flow and nonlinear optical pulse propagation [22]. The travelling-wave solution of (1) discussed below therefore has a direct physical interpretation: its wave speed $A_2$ is the propagation velocity of an action potential in the nerve-pulse interpretation, of an invasion front in the population-dynamics interpretation, and of a shock wave in the nonlinear-wave interpretation. Quantitative prediction in any of these domains requires benchmark solutions accurate well beyond the precision of the physical measurements they model, which is the numerical gap the present work addresses.

The gBH equation admits the travelling-wave exact solution

$$u(x,t)={\left[\frac{\gamma }{2}+\frac{\gamma }{2}tanh\left(A_1(x-A_{2}t)\right)\right]}^{1/\delta },$$

(2)

with amplitude

$$A_1={\displaystyle \frac{\gamma \delta \left(-\alpha +\sqrt{{\alpha }^2+4\beta (1+\delta )}\right)}{4(1+\delta )}}.$$

(3)

For the wave speed $A_2$, two formulas appear in the literature. The Wang–Zhu–Lu formula [1], obtained by substituting the ansatz into (1) and balancing terms,

$$A_2^{\mathrm{W}}={\displaystyle \frac{\alpha \gamma }{1+\delta }}-{\displaystyle \frac{(1+\delta -\gamma )(-\alpha +\sqrt{{\alpha }^2+4\beta (1+\delta )})}{2(1+\delta )}},$$

(4)

is the formula most frequently cited in numerical-methods papers that use (2) as a reference solution. The alternative formula

$$A_2={\displaystyle \frac{\alpha \gamma }{1+\delta }}+{\displaystyle \frac{(1+\delta -\gamma )(\alpha +\sqrt{{\alpha }^2+4\beta (1+\delta )})}{2(1+\delta )}},$$

(5)

was derived by Deng [8] using the first-integral method and later verified symbolically by Appadu and Tijani [9].

Direct substitution of (2)–(3) with $A_2=A_2^{\mathrm{W}}$ into (1) produces a non-vanishing residual; the wave profile generated by $A_2^{\mathrm{W}}$ is not a solution of the gBH equation. The difference between the two formulas is not cosmetic. At the canonical parameter set $\alpha =\beta =\gamma =\delta =1$, (4) returns $A_2^{\mathrm{W}}=0$—a “travelling wave” that does not travel—while (5) returns $A_2=3/2$. At the Ismail–Raslan–Rabboh benchmark [3] ($\alpha =\beta =\delta =1$, $\gamma ={10}^{-3}$), (4) returns $A_2^{\mathrm{W}}\approx -0.999$, a wave travelling in the opposite direction from the true solution with $A_2\approx 1.9995$.

Despite the published corrections [8,9], (4) continues to be used as the reference solution in numerical studies, and reported “errors” in the range ${10}^{-4}$–${10}^{-8}$ are routinely interpreted as measures of the schemes’ numerical accuracy. A representative timeline of such studies includes: the Adomian-decomposition benchmark of Ismail, Raslan and Rabboh [3] reporting $3.75\times {10}^{-7}$; the Adomian work of Hashim et al. [4] reporting similar values; the B-spline collocation studies of Çelik [5] and Mohammadi [6] in the ${10}^{-8}$ range; and, most recently, the Haar-wavelet study of Shukla and Kumar [7] with errors of order ${10}^{-7}$. In each case the quoted error coincides to several significant digits with the analytical wave-profile gap $\parallel u_{A_2}-u_{A_2^{\mathrm{W}}}{\parallel }_{\infty }$ rather than with any genuine numerical defect.

This paper makes two contributions. First, we develop an ETDRK4–Chebyshev pseudospectral scheme for (1) with non-homogeneous Dirichlet data. The scheme combines Chebyshev–Gauss–Lobatto collocation in space with linear boundary lifting to homogenize the boundary data, and the Cox–Matthews/Kassam–Trefethen ETDRK4 exponential integrator in time [13,14], with matrix $\phi$-functions evaluated by the Trefethen–Weideman–Schmelzer Talbot contour quadrature [15]. On the Ismail–Raslan–Rabboh benchmark the scheme attains floating-point round-off accuracy with $N=2$ collocation points and a single time step of size $\Delta t=1.0$, a combination of accuracy, spatial economy and temporal stride that exceeds the performance of previously reported schemes for this equation by at least twelve orders of magnitude. The method complements but does not duplicate existing exponential-integrator work on the gBH equation: ETDRK4 with finite differences was applied by Owolabi [10]; bivariate Chebyshev collocation with quasilinearization rather than exponential time stepping was developed by Motsa et al. [11]; and barycentric Jacobi collocation with an ETD predictor–corrector was proposed by Pindza, Owolabi and Patidar [12]. None of these prior works combines the ETDRK4 time stepper with one-dimensional Chebyshev collocation on a bounded domain with linear boundary lifting.

Second, the extreme accuracy of the scheme turns it into an analytical probe. Measured against the corrected exact solution (5), the scheme’s error at the Ismail benchmark is exact to floating-point round-off at $t=1$ and at $t=10$. Measured against Wang’s (4), the same scheme plateaus at a constant independent of both N and $\Delta t$. We derive this plateau in closed form in Section 2: the residual produced by Wang’s formula factors as $R=\frac{\gamma A_1}{2}(A_2-A_2^{\mathrm{W}})(1-v^2)$, whose $L^{\infty }$ norm equals the analytical distance between the two wave profiles. At the Ismail benchmark this distance is $3.748\times {10}^{-7}$, matching the observed plateau to eight significant digits. The published ${10}^{-7}$-range errors of [3,4,7] are consequently dominated by this analytical gap rather than by any numerical defect in those schemes.

The remainder of the paper is organized as follows. Section 2 verifies the correctness of (5) by full coefficient-matching and derives the closed-form residual associated with Wang’s formula. Section 3 presents the ETDRK4–Chebyshev scheme with boundary lifting. Section 4 reports the numerical experiments—including three new convergence studies at different nonlinearity regimes, a visualization of the wave-profile gap, and a comparison with published benchmarks. Section 5 concludes with implications for future benchmarking work on the gBH equation.

2. Verification of the Travelling-Wave Solution

We verify (5) directly by substituting the ansatz (2) into the gBH equation (1) and requiring that each monomial in the auxiliary variable vanish. Set

$$U(x,t):=u^{\delta }(x,t)={\displaystyle \frac{\gamma }{2}}(1+v),v:=tanh\left(A_1(x-A_{2}t)\right).$$

(6)

Then $u=U^{1/\delta }$, and differentiation gives

$$U_t=-{\displaystyle \frac{\gamma A_1{A}_2}{2}}(1-v^2),U_x={\displaystyle \frac{\gamma A_1}{2}}(1-v^2),U_{xx}=-\gamma A_1^{2}v(1-v^2).$$

(7)

Using $U=u^{\delta }$, the chain rule gives $U_t=\delta u^{\delta -1}{u}_t$ so $u_t=U_t/\left(\delta u^{\delta -1}\right)$, and similarly for $u_x$. Substituting these into (1), multiplying through by $\delta u^{\delta -1}$, and simplifying (including the identity $u_{xx}=U_{xx}/\left(\delta u^{\delta -1}\right)-(\delta -1)U_x^2/\left(\delta u^{2\delta -1}\right)$), after algebraic manipulation the equation becomes, in the variable v,

$$-{\displaystyle \frac{\gamma A_2}{2}}(1-v^2)+\alpha {\displaystyle \frac{\gamma }{2}}(1+v){\displaystyle \frac{\gamma A_1}{2}}(1-v^2)·{\displaystyle \frac{1}{\delta }}-{\displaystyle \frac{U_{xx}}{\delta }}=\beta U(1-U)(U-\gamma )·{\displaystyle \frac{1}{\delta u^{\delta -1}}}·\delta .$$

(8)

After normalizing by $\gamma /\left(8\delta \right)$ and collecting powers of v, the residual reduces to a polynomial $P\left(v\right)$ of degree at most three. For $\delta =1$ only the coefficients of $v^0$ and $v^1$ survive, and vanishing of these two coefficients yields two algebraic conditions: the $v^1$ coefficient fixes $A_1$ as in (3), while the $v^0$ coefficient determines $A_2$. Carrying out the algebra, the $v^0$ condition reads

$$-A_2+{\displaystyle \frac{\alpha \gamma }{1+\delta }}+{\displaystyle \frac{(1+\delta -\gamma )\left(\alpha +\sqrt{{\alpha }^2+4\beta (1+\delta )}\right)}{2(1+\delta )}}=0,$$

(9)

which gives (5). For general $\delta >0$ the same argument applies after accounting for the extra factor of $u^{\delta -1}$; the final form of $A_2$ is unchanged.

The Wang residual. 

If the same ansatz is substituted with $A_2=A_2^{\mathrm{W}}$ as given in (4) rather than (5), only the $v^0$ coefficient fails to vanish. Tracking this coefficient through the algebra, the full residual takes the Remarkably simple closed form

$$R\left(v\right)={\displaystyle \frac{\gamma A_1}{2}}\left(A_2-A_2^{\mathrm{W}}\right)(1-v^2),$$

(10)

which is the product of the amplitude difference and the hyperbolic-secant-squared envelope. This explicit closed-form residual is a new observation, not previously reported in the literature.

Proposition 1. 1

Let $A_1$ be given by (3) and let $A_2,A_2^{\mathrm{W}}$ be given by (5),(4). For the travelling-wave ansatz $u_W(x,t)={\left[\frac{\gamma }{2}(1+tanh{A}_1(x-A_2^{\mathrm{W}}t))\right]}^{1/\delta }$, the residual $R:={\partial }_t{u}_W+\alpha u_W^{\delta }{\partial }_x{u}_W-{\partial }_{xx}{u}_W-\beta u_W(1-u_W^{\delta })(u_W^{\delta }-\gamma )$ satisfies

$${\parallel R(·,t)\parallel }_{\infty }={\displaystyle \frac{\gamma A_1}{2}}\left|A_2-A_2^{\mathrm{W}}\right|,\mathrm{for}\mathrm{all}t\ge 0.$$

At the Ismail–Raslan–Rabboh benchmark ($\alpha =\beta =\delta =1$, $\gamma ={10}^{-3}$), we find $A_1=2.500\times {10}^{-4}$, $A_2=1.9995$, $A_2^{\mathrm{W}}=-0.9990$, and ${\parallel R\parallel }_{\infty }=\frac{\gamma A_1}{2}|A_2-A_2^{\mathrm{W}}|=3.748\times {10}^{-7}$. At leading order, the same prefactor controls the $L^{\infty }$ distance between the two travelling-wave profiles:

$$\parallel u_{A_2}-u_{A_2^{\mathrm{W}}}{\parallel }_{\infty }\approx {\displaystyle \frac{\gamma A_1}{2}}\left|A_2-A_2^{\mathrm{W}}\right|=3.748\times {10}^{-7}.$$

(11)

This is precisely the N- and $\Delta t$-independent plateau observed in Section 4 when the scheme is benchmarked against Wang’s profile, and it matches to eight significant digits the “errors” reported by [3,4,7] at the Ismail benchmark.

Remark 1. 

The closed-form residual (10) has a useful interpretation: because $R\left(v\right)\propto (1-v^2)$, the error concentrates around the centre of the front and vanishes exponentially in the tails. Any numerical method that benchmarks against $u_{A_2^{\mathrm{W}}}$ will therefore exhibit a characteristic bell-shaped error profile whose peak is $3.748\times {10}^{-7}$ independent of discretization. This is visible in Figure 1 below.

3. ETDRK4–Chebyshev Scheme

On the interval $[a,b]$ we use Chebyshev–Gauss–Lobatto nodes $x_k=\frac{a+b}{2}+\frac{b-a}{2}cos(k\pi /N)$, $k=0,\dots ,N$, and construct the associated differentiation matrices $D,D^{\left(2\right)}=D^2$ by the barycentric formula [17]. To homogenize the non-homogeneous Dirichlet conditions $u(a,t)=g_L\left(t\right)$, $u(b,t)=g_R\left(t\right)$, we introduce the linear boundary lift

$$\ell (x,t)={\displaystyle \frac{b-x}{b-a}}{g}_L\left(t\right)+{\displaystyle \frac{x-a}{b-a}}{g}_R\left(t\right),$$

(12)

and write $u=v+\ell$ with v vanishing at both endpoints. The motivation for the lifting approach, rather than eliminating the boundary rows from the pseudospectral system, is that the linear lift has a known time derivative and diffusive contribution, which can be evaluated analytically and added to the right-hand side at each stage. This is numerically more stable than boundary-row elimination, which introduces spurious eigenvalues at the corners of the domain.

Substituting $u=v+\ell$ into (1), restricting to the $N-1$ interior nodes, and using $L_{\mathrm{int}}:=D_{1:N-1,1:N-1}^{\left(2\right)}$ to denote the interior block of the Laplacian, the semi-discrete system is

$$\dot{\mathbf{v}}=L_{\mathrm{int}}\mathbf{v}+\tilde{F}(\mathbf{v},t),$$

(13)

where $\mathbf{v}\in {\mathbb{R}}^{N-1}$ collects the interior values of v and the nonlinear-plus-lift source is

$$\tilde{F}(\mathbf{v},t)=-\alpha {(\mathbf{v}+\ell )}^{\delta }\odot {\mathbf{u}}_x+\beta (\mathbf{v}+\ell )\odot [1-{(\mathbf{v}+\ell )}^{\delta }]\odot [{(\mathbf{v}+\ell )}^{\delta }-\gamma ]-{\partial }_t\ell +(L\ell ){|}_{\mathrm{int}}.$$

(14)

Here ${\mathbf{u}}_x$ is the interior values of ${\partial }_{x}u$ computed from the full $(N+1)$-vector ${(g_R,\mathbf{v}+\ell ,g_L)}^{\top }$ via D, with the boundary columns of D contributing $g_R$ and $g_L$ directly (and similarly for ${(L\ell )|}_{\mathrm{int}}$ through the boundary columns of $D^{\left(2\right)}$), and ⊙ denotes Hadamard (elementwise) multiplication.

The system (13) is advanced by ETDRK4 [13,14]:

$$\begin{array}{cc}\hfill \mathbf{a}& ={\phi }_0\left(\frac{h}{2}{L}_{\mathrm{int}}\right){\mathbf{v}}_n+\frac{h}{2}{\phi }_1\left(\frac{h}{2}{L}_{\mathrm{int}}\right)\tilde{F}({\mathbf{v}}_n,t_n), \hfill & (15a)\\ \hfill \mathbf{b}& ={\phi }_0\left(\frac{h}{2}{L}_{\mathrm{int}}\right){\mathbf{v}}_n+\frac{h}{2}{\phi }_1\left(\frac{h}{2}{L}_{\mathrm{int}}\right)\tilde{F}(\mathbf{a},t_n+\frac{h}{2}),\hfill & (15b)\\ \hfill \mathbf{c}& ={\phi }_0\left(\frac{h}{2}{L}_{\mathrm{int}}\right)\mathbf{a}+\frac{h}{2}{\phi }_1\left(\frac{h}{2}{L}_{\mathrm{int}}\right)[2\tilde{F}(\mathbf{b},t_n+\frac{h}{2})-\tilde{F}({\mathbf{v}}_n,t_n)],\hfill & (15c)\\ \hfill {\mathbf{v}}_{n+1}& ={\phi }_0\left(h{L}_{\mathrm{int}}\right){\mathbf{v}}_n+h({\phi }_1-3{\phi }_2+4{\phi }_3)\left(h{L}_{\mathrm{int}}\right)\tilde{F}({\mathbf{v}}_n,t_n)\hfill \\ \hfill & +2h({\phi }_2-2{\phi }_3)\left(h{L}_{\mathrm{int}}\right)[\tilde{F}(\mathbf{a},t_n+\frac{h}{2})+\tilde{F}(\mathbf{b},t_n+\frac{h}{2})]\hfill \\ \hfill & +h(4{\phi }_3-{\phi }_2)\left(h{L}_{\mathrm{int}}\right)\tilde{F}(\mathbf{c},t_n+h),\hfill & (15d)\end{array}$$

where ${\phi }_0\left(z\right)=e^z$ and ${\phi }_k\left(z\right)={\int }_0^1{e}^{(1-s)z}{s}^{k-1}/(k-1)!ds$ for $k\ge 1$.

Evaluating the $\phi$-functions.

Because $L_{\mathrm{int}}$ is dense and non-normal, the scalar Kassam–Trefethen approach based on eigendecomposition is unsuitable. We instead use the modified Talbot contour quadrature of Trefethen, Weideman and Schmelzer [15], which expresses each matrix $\phi$-function as a Cauchy integral on a parabolic-hyperbolic contour in the resolvent half-plane:

$${\phi }_k\left(hA\right)={\displaystyle \frac{1}{2\pi i}}{\int }_{\Gamma }{e}^{hs}{s}^{-k}{(sI-A)}^{-1}ds,$$

(16)

evaluated by the midpoint rule on the Talbot parameterization $s=s\left(\theta \right)$. With $n_q=64$ quadrature nodes this yields all four functions ${\phi }_0,{\phi }_1,{\phi }_2,{\phi }_3$ to double precision in a single vectorized solve per node. The precomputation is performed once per distinct $(h,L_{\mathrm{int}})$ pair and reused for all time steps.

3.0.0.3. Computational cost.

The dominant cost per time step is six applications of the full ${(N-1)}^2$ dense $\phi$-operators. The precomputation cost is $\mathcal{O}\left(n_q{(N-1)}^3\right)$. For the benchmarks reported in Section 4 this is negligible; total runtime for the full Ismail benchmark with $N=20$ and $\Delta t={10}^{-2}$ over $[0,1]$ is under one second on a laptop.

4. Numerical Experiments

We test the scheme on two benchmark regimes from the literature:

• R1 (Ismail benchmark [3]): $\alpha =\beta =\delta =1$, $\gamma ={10}^{-3}$. This is the most frequently cited gBH benchmark.
• R2 (Hashim benchmark [4]): $\alpha =\beta =1$, $\delta =2$, $\gamma ={10}^{-2}$. This regime has higher nonlinearity in the $\delta$ parameter.

Both are solved on $\Omega =[0,1]$ with Dirichlet data read from the exact travelling-wave solution (2), (5), at $x=0$ and $x=1$.

4.1. Single-Step Machine Precision and Comparison with Published Errors

Figure 2(a) shows the $L^{\infty }$ error at $t=1$ as a function of N, using a single ETDRK4 step $\Delta t=1.0$. Measured against the correct solution (5), the error reaches the level of floating-point round-off ($\sim {10}^{-19}$ absolute, or $\sim {10}^{-15}$ relative since $u\sim {10}^{-4}$ on this benchmark) already at $N=2$. Measured against Wang’s (4), the error plateaus at $3.748\times {10}^{-7}$, exactly the analytical gap predicted in (11).

Table 1 compares the scheme’s accuracy with published results for R1 and R2. All published errors shown are the values reported in the respective papers at the same parameter sets, read at $t=1$ (R1, R2) and $t=10$ (R1 only). The separation between our accuracy and the Wang-plateau accuracy of prior schemes at R1 is twelve orders of magnitude.

4.2. Visualization of the Analytical Wave-Profile Gap

Figure 1 makes the diagnostic interpretation of the $3.748\times {10}^{-7}$ plateau visually explicit. The top panel overlays the two travelling-wave profiles $u_{A_2}(x,1)$ and $u_{A_2^{\mathrm{W}}}(x,1)$ on the Ismail benchmark; on the ${10}^{-4}$ scale of the solution, the two profiles are visually indistinguishable. The bottom panel plots the pointwise difference $|u_{A_2}(x,1)-u_{A_2^{\mathrm{W}}}(x,1)|$ on a log scale. The difference attains its maximum of $3.748\times {10}^{-7}$ near the centre of the front, matching the analytical prediction from (11) to eight significant digits.

4.3. Temporal Order Reduction Across Nonlinearity Regimes

Figure 2(b) shows a $\Delta t$-refinement study at $\gamma =0.5$ (a strongly nonlinear regime) with $N=30$ fixed. The observed order is approximately $\mathcal{O}\left(\Delta t^{2.5}\right)$, well below the classical fourth order expected of ETDRK4 in smooth-data settings. This is the classical Hochbruck–Ostermann order reduction [16]: exponential integrators applied to parabolic problems with non-homogeneous (in particular, time-dependent) Dirichlet boundary data generically exhibit order reduction from p to $p-3/2$ because the algebraic boundary conditions are only first-order accurate in the stage residuals.

To confirm that this is a robust feature of the gBH equation rather than a coincidence of one parameter set, Figure 3 reports the same refinement study at four gamma values spanning weakly to strongly nonlinear regimes: $\gamma \in \{0.1,0.3,0.5,0.9\}$. Using a linear fit to the three smallest $\Delta t$ values, the observed convergence orders are

$$p_{0.1}=2.46,p_{0.3}=2.44,p_{0.5}=2.45,p_{0.9}=2.47,$$

tightly clustered around $2.45$ and fully consistent with the Hochbruck–Ostermann prediction. The absolute error levels grow with $\gamma$ (as expected, since the solution amplitude scales as ${\gamma }^{1/\delta }$) but the convergence order is essentially $\gamma$-independent.

4.4. Long-Time Integration

The Ismail benchmark at $t=10$ in Table 1 shows no accumulation of error: the scheme’s $L^{\infty }$ error at $t=10$ is $5.4\times {10}^{-19}$, essentially identical to the value at $t=1$. This is because the Chebyshev spatial discretization resolves the analytic solution exactly (to round-off), so temporal error is the only source, and ETDRK4 at small enough $\Delta t$ gives errors that remain at round-off. In contrast, a published scheme whose error is dominated by the Wang-formula plateau will exhibit the same $3.748\times {10}^{-7}$ “error” at any time t, which is a tell-tale signature of the analytical gap rather than genuine numerical accumulation.

5. Conclusions and Future Directions

Coupling ETDRK4 with Chebyshev collocation and linear boundary lifting produces a scheme for the generalized Burgers–Huxley equation that is exact to floating-point round-off on the Ismail–Raslan–Rabboh benchmark with $N=2$ collocation points and a single time step of size $\Delta t=1.0$, and that exhibits $\mathcal{O}\left(\Delta t^{2.45}\right)$ Hochbruck–Ostermann-reduced temporal order across four $\gamma$-regimes. Beyond its value as a method, the scheme resolves a decades-old issue in the benchmarking literature for this equation: Wang’s wave-speed formula, still dominant in numerical-methods papers as recently as 2020, does not satisfy the PDE. The corrected formula of Deng (2008), verified in Appadu and Tijani (2022), is the unique value for which the travelling-wave ansatz is a solution, and the residual generated by Wang’s formula admits the closed form $\frac{\gamma A_1}{2}(A_2-A_2^{\mathrm{W}})(1-v^2)$, whose $L^{\infty }$ maximum equals the analytical wave-profile gap.

Consequently, any benchmark that uses $A_2^{\mathrm{W}}$ as the reference exhibits an N- and $\Delta t$-independent error floor equal to $\parallel u_{A_2}-u_{A_2^{\mathrm{W}}}{\parallel }_{\infty }$. Reported errors at or near this floor (including the widely-cited $3.748\times {10}^{-7}$ value from the Ismail–Raslan–Rabboh series) reflect the wave-profile gap rather than the accuracy of the numerical scheme. Method comparisons on the gBH equation should henceforth use the corrected formula (5) as the reference solution, and published “errors” against Wang’s formula at or above $3.748\times {10}^{-7}$ at the Ismail benchmark (or the corresponding gap at other parameter sets) should be reinterpreted as measurements of the analytical gap rather than of scheme accuracy.

Several directions for future work suggest themselves. First, the closed-form residual (10) depends only on the amplitude difference $|A_2-A_2^{\mathrm{W}}|$ and not on the spatial discretization: an analogous analysis for other nonlinear wave equations with disputed or variant exact solutions (the Fisher equation, the Fitzhugh–Nagumo equation) may yield similar diagnostics. Second, extending the ETDRK4–Chebyshev approach to two-dimensional generalizations of the gBH equation and to systems of coupled reaction-diffusion equations would be of interest for applications in excitable media and population dynamics. Third, the Hochbruck–Ostermann order reduction could be mitigated by using structure-preserving ETDRK variants that account for the boundary data at each stage; constructing such variants for the Chebyshev-collocation setting remains an open problem.