Structural Obstruction in the Romanovski $\beta$ = 1 Compact-Interval Case and Domain-Decomposed Collocation for the Gompertz–Burgers Tumor Invasion Model

The Romanovski weight family $\rho(x ; \alpha, \beta)=\left(1+x^2\right)^{\beta-1} e^{\alpha \arctan (x)}$ generates orthogonal polynomials that appear in quantum-mechanical potential wells, quark models, and random matrix theory [1-3]. The classical theory is well-developed for $\beta < 0$ on the real line [4-5]. The case $\beta = 1$, $w(x)=\pi^{\arctan (x)}$ on the compact interval $[-1,1]$, has a natural motivation (bounded and strictly positive weight on a compact domain arising from transport-diffusion models [6]) but falls outside the standard regime. The precise structural reason for this has not been characterized explicitly in the literature.

Two specific gaps motivated the present study. Analytically, although the Romanovski family is well understood for $\beta < 0$ on the real line, the literature did not state explicitly why the $\beta = 1$ compact-interval case admits no Rodrigues or Sturm–Liouville representation, nor does it supply the corresponding orthogonal polynomials by an alternative route. Computationally, the numerical consequences of employing such a weight in a nonlinear transport--diffusion model, particularly whether the polynomial family matters at all for the attainable accuracy, have not been examined. The present work addressed both gaps directly.

This paper made two contributions. First, we proved that the Rodrigues and Sturm–Liouville constructions fail for $\beta = 1$ on any compact interval, identified this consequent to $\sigma(x) = 1 + x^2$ having only complex roots, and constructed the orthogonal polynomials directly from the three-term recurrence. Second, we applied these polynomials to the Gompertz–Burgers tumor invasion model [7-8] and found that the basis choice was irrelevant to numerical accuracy (Section 5). We identified the dominant practical obstacle: a transition layer of width $O(\nu)$ that rendered all single-domain collocation schemes unreliable at physiological parameters, regardless of polynomial basis. We characterized and remedied this failure through domain decomposition (Section 6), to produce the major numerical result.

The remainder of the paper is organised as follows. Section 2 defines the $\beta = 1$ weight and proves the structural obstruction (Theorem 1). Section 3 constructs the Keşan polynomials from the three-term recurrence. Section 4 sets up the Gompertz–Burgers tumor invasion model and its reference solution. Section 5 establishes that the single-domain collocation failure is independent of the polynomial basis. Section 6 introduces the domain-decomposition remedy, analyses interface placement, and reports verified negative results. Section 7 compares the findings with recent literature, and Section 8 concludes.

Let $w(x)=\pi^{\arctan (x)}$ on $[-1,1]$. Writing $\alpha = \ln\pi$, this is $w(x)=e^{\alpha \arctan (x)}$, the $\beta = 1$ case of the Romanovski family. Differentiating, $(1 + x^2)w'(x) = \alpha w(x)$. The Pearson equation $\frac{\mathrm{d}}{\mathrm{d} x}[\sigma(x) w(x)]=\tau(x) w(x)$, with $\deg(\sigma) \leq 2$, $\deg(\tau) = 1$, yields the unique solution $\sigma(x) = 1 + x^2$ and $\tau(x) = 2x + \alpha$. Uniqueness follows because any deviation in the coefficient of 1 in $\sigma$ leaves a genuine rational remainder.

Theorem 1. For the weight $w(x)=\pi^{\arctan (x)}$ and Pearson polynomial $\sigma(x) = 1 + x^2$, the boundary term $[\sigma(x) w(x)]_{-1}^1=2 \pi^{\pi / 4}-2 \pi^{-\pi / 4} \approx 4.101$, strictly nonzero. Consequently: (i) no iterated integration-by-parts in the Rodrigues formula reduced to a boundary-free expression; (ii) the differential operator $L$ = $\sigma(x) \frac{\mathrm{d}^2}{\mathrm{~d} x^2}+\tau(x) \frac{\mathrm{d}}{\mathrm{d} x}$ was not self-adjoint on $L^2([-1,1],w)$; and (iii) no second-order differential equation with a constant eigenvalue $\lambda_n$ per degree can hold for the monic orthogonal polynomials of this weight.

Proof. Direct computation at $x= \pm 1$, using $w(1)=\pi^{\pi / 4}$, $w(-1)=\pi^{-\pi / 4}$, $\sigma( \pm 1)=2$, gives the stated nonzero value. Part (i)–(iii) follow from standard Sturm–Liouville theory [9].

Remark 1. The obstruction is generic: when $\sigma$ has only complex roots ($\sigma=1 + x^2$, roots $\pm i$), integration produces the arctangent, yielding a transcendental weight that cannot vanish at any pair of real points. The Romanovski $\beta = 1$ case is the canonical instance on a compact interval.

Remark 2. We verified for $\beta=-2, \alpha=1$ on the unbounded interval that the Rodrigues formula produces the correct polynomials and the Sturm–Liouville equation holds with eigenvalue $\lambda_1=-4$, residual $<10^{-12}$. The Romanovski theory is sound; $\beta = 1$ on a compact interval lies outside its domain of validity.

Since Theorem 1 foreclosed the Rodrigues route, we constructed the orthogonal polynomials $\left\{K_{\mathrm{n}}\right\}$ directly from their moments $\mu_k=\int_{-1}{ }^1 x^k w(x) d x$ (all finite by $w>0$ bounded on $[-1,1]$) via Gram–Schmidt orthogonalisation [10], giving the three-term recurrence $K_{\mathrm{n}+1}(x)=\left(x-\beta_{\mathrm{n}}\right) K_{\mathrm{n}}(x)-\gamma_{\mathrm{n}} K_{\mathrm{n}-1}(x)$, with $\beta_n=\frac{\left\langle x K_n, K_n\right\rangle_w}{\left\langle K_n, K_n\right\rangle_w}$ and $\gamma_n=\frac{\left\langle K_n, K_n\right\rangle_w}{\left\langle K_{n-1}, K_{n-1}\right\rangle_w}$. By Favard’s theorem, the resulting polynomials were orthogonal with respect to $w$, regardless of the unavailability of a differential representation [11]. The resulting recurrence coefficients $\beta_{\mathrm{n}}$ and $\gamma_{\mathrm{n}}$, together with the squared norms $h_n=\left\|K_n\right\|_w^2$ and the numerically verified cross-orthogonality $\left|\left\langle K_n, K_m\right\rangle_w\right| \quad(m \neq n)$, are listed in Table 1 for degrees $n=0$ to $4$.

The asymmetry in $\beta_n$ reflects the asymmetry of $w$: $\frac{w(1)}{w(-1)} \approx 6.04$. The sequence $\gamma_{\mathrm{n}} \rightarrow 0.25$ is consistent with the spectral radius of the Jacobi matrix on this interval [12], confirming numerical stability of the recurrence.

We applied the collocation machinery to the nonlinear boundary value problem [6-8]:

The parameter values represent normalized and biologically representative conditions rather than a fit to a specific dataset. The reaction term $-\alpha u \ln u$ is the Gompertz growth law [13] with relative growth rate $\alpha=0.5$, an $O(1)$ normalized value consistent with the Gompertz-type tumor-growth descriptions of [7-8] [14]. The dependent variable u is a normalized cell density: the boundary value $u(-1)=1$ corresponds to the maximal (carrying-capacity) density, while $u(1)=0.36 \approx e^{-1}$ is the natural Gompertz reference level at which the relative-growth factor $\ln u$ equals $-1$. The diffusion coefficient $v=10^{-3}$ places the model in the strongly convection-dominated regime [15] characteristic of invasive-front dynamics [6], in which the transport term $\left(1+x^2\right) u^{\prime}$ dominates diffusion and a thin transition layer. This small value was chosen deliberately to expose the numerical difficulty that was the subject of the present study.

An independent reference solution was computed via $scipy.integrate.solve_bvp$ with continuation in $v$ from 0.5 to $10^{-3}$ (six stages). The solution was $u \approx 1$ for $x$ up to $\approx 0.998$, then fell sharply to $u=0.36$ within $[ 0.999,1.0]$, a transition layer of width $\approx 0.001$ under 0.2\% of the domain, scaling as $O(v)$ (convective layer balance). This extreme sharpness was the source of all subsequent numerical difficulties.

Newton–Raphson collocation with $v$-continuation was implemented for the Keşan, Legendre, and Chebyshev bases [16] under identical conditions. All three bases produced indistinguishable results: weighted mean-squared error (WMSE) against the reference solution remained in the range 0.37–0.73 across $N=16$ to $N=48$ for all three bases, and did not decrease monotonically with $N$. The Keşan and Legendre bases agreed to be within the third decimal place at every degree. The failure was intrinsic to single-domain polynomial collocation [17] on this problem as it was not basis-dependent. Figure 1 (panel a) shows the non-physical zigzag oscillation of the single-domain $N = 16$ solution across the full domain.

We split $[-1,1]$ into $\left[-1, x_{\text {split }}\right]$ (bulk) and $\left[x_{\text {split }}, 1\right]$ (layer), each with an independent Chebyshev–Gauss–Lobatto grid [18], coupled by $C^1$ continuity (shared value and matched derivative) at $x_{\text {split }}=0.99$, solved by full-Jacobian Newton–Raphson with identical continuation. This single change reduced the WMSE from 0.37 (best single-domain, $N=48$) to $0.0195$ ($N_{\text {total }}=23$), a factor of 19–21 improvement. Table 2 summarizes all tested configurations. Figure 1 (panel b, zoom into the transition layer) shows that the two-domain solution closely tracked the reference solution in the sharpest part of the domain, while the single-domain scheme collapsed.

Two mechanisms explained this behaviour. First, the single-domain scheme should represent both the flat bulk, where $u \approx 1$ over more than 99% of the domain, and the near-discontinuous dropped within $[ 0.999,1]$ using a single global polynomial. The fixed global node distribution placed far too few points inside the layer and the interpolant developed the non-physical zigzag of panel (a) as it attempted to fit an under-resolved sharp gradient. Second, the two-domain strategy assigned a dedicated Chebyshev–Gauss–Lobatto grid to the narrow layer subdomain $\left[x_{\text {split }}, 1\right]$, where Gauss–Lobatto clustering nodes concentrated exactly where the gradient was largest, while the bulk subdomain required only a low-order representation of an essentially flat profile. The interface $x_{\text {split }}=0.99$ was effective because it lied just outside the layer, which occupied $[ 0.999,1]$. It kept the entire steep transition inside the refined layer subdomain, while leaving the bulk subdomain free of any sharp features. Placing the interface inside the layer split the gradient across the $C^1$ coupling, whereas placing it too far from the layer wasted resolution on the flat region. Both degraded the result, which was precisely the sensitivity quantified by the Shishkin-parameter analysis in Section 6.2.

The classical Shishkin formula [15] for convection-dominated layers, $\tau=\sigma v \ln \left(N_2\right)$, $x_{\text {split }}=1-\tau$, with the literature standard $\sigma=1$ (derived for linear problems) gave $x_{\text {split }}\approx 0.997$, underperforming hand-tuning by $\approx 40 \%$ (WMSE 0.027 vs.0.0195). Calibrating $\sigma \approx 4$ by a one-time search matched the hand-tuned result and generalised more consistently across $N$. The recalibration reflects this problem’s nonlinear reaction term $\frac{\partial f}{\partial u}<0$ near $u=1$, which violated the sign assumption of standard semi-linear Shishkin theory [19].

Three-domain decomposition consistently underperformed two-domain at matched $N$ (best WMSE 0.0276 at $N = 30$ vs. 0.0172 at $N = 29$). Legendre collocation within the two-domain framework was markedly less reliable at low $N_1$, attributed to Newton convergence to non-physical solutions more readily under Legendre nodes in the bulk subdomain. Hermite- and Laguerre-mapped collocation on $[-1,1]$ produced spurious eigenvalues at all tested $N$ due to node clustering [20]. A simplified two-point overlap method from [21] performed substantially worse (WMSE $\approx 0.46–0.48$). All were reported as verified negative results.

The computational workload of every scheme was dominated by the dense Newton solve, repeated over a small number of iterations within each of the six $v$-continuation stages. The linear systems involved were small: dimension $N+1$ ($N = 16–48$) for the single-domain schemes, and $N_1+N_2+2$ (23–29 unknowns, including the two interface-continuity constraints) for the two-domain scheme. The absolute cost of all configurations was therefore minor. The point relevant to the present comparison was relative: the two-domain scheme reached its target accuracy at a total of about $23$ unknowns, whereas the best single-domain configuration used $N=48$ and still did not attain comparable accuracy. Domain decomposition thus improved accuracy while operating on a smaller linear system, so its advantage carried no additional computational penalty.

Yang and Srinivasan [21] independently reported the same qualitative behaviour for the viscous Burgers equation at $v=5 \times 10^{-3}$. A single global Chebyshev collocation with $N = 300$ nodes places too few points in the region of sharpest gradient and attains only a peak error of order $10^{-4}$, whereas their overlapping domain-decomposition with nonlinear mapping (of the order of 150 nodes per subdomain) reduces the error by several orders of magnitude at a comparable total node count, and improves further as points are added to the overlap region. Although their smoother, larger-$v$ Burgers layer differs from the present convection-dominated case at $v=5 \times 10^{-3}$, the convergence of independent findings on related but distinct problems confirms that single-domain polynomial collocation is generically ill-suited to narrow-layer nonlinear problems of this class, and that domain decomposition is the effective remedy.

Two independent results have been established in this paper and they are: (1) Structural: the Romanovski $\beta=1$ weight on any compact interval cannot support a Rodrigues representation or constant-eigenvalue Sturm–Liouville equation. This is a consequence of $\sigma(x)=1+x^2$ having only complex roots (Theorem 1). The Keşan orthogonal polynomials are nevertheless well-defined and computable through the three-term recurrence (Table 1); and (2) Computational: for the Gompertz–Burgers tumor invasion model at $v=10^{-3}$, all single-domain collocation schemes fail irrespective of polynomial basis. Two-domain Chebyshev decomposition restores accuracy with a 19–21-fold reduction in WMSE (Table 2 and Figure 1). The optimal interface placement is most robustly specified by the calibrated Shishkin formula ($\sigma \approx 4$) rather than $\sigma=1$. The choice of orthogonal polynomial basis is immaterial to these conclusions.

Several limitations bounded these conclusions and indicated directions for future work. The structural theorem is specific to the $\beta=1$ weight with $\sigma(x)=1+x^2$ although the underlying mechanism, a Pearson polynomial with only complex roots is generic. The full $\beta$-range on compact intervals has not been characterized here. The computational study addressed a single one-dimensional, steady, and scalar Gompertz–Burgers model with one fixed set of normalized parameters, so the $\approx20$-fold gain from two-domain decomposition, the immateriality of the polynomial basis, and the recalibrated Shishkin constant ($\sigma \approx 4$) should be regarded as the representative of this setting rather than universal. Natural extensions include time-dependent invasion fronts, systems of coupled tumor–host equations, two spatial dimensions where the interface becomes a curve rather than a point, and adaptive rather than a priori interface selection. A rigorous error analysis of the two-domain scheme for this nonlinear reaction term, in the spirit of the semi-linear Shishkin theory likewise remains open.