A Hybrid Legendre–Laguerre Operational Matrix Method for the Variable-Coefficient Fifth-Order KdV Equation in Non-Homogeneous Optical Media

The phenomenon of trapping light within microscopic domains has emerged as a cornerstone of modern photonics, driving advancements in optical storage, quantum information processing, and slow-light switches [1], [2]. Traditional approaches rely primarily on static geometric modifications, such as engineered photonic crystal defects or metamaterial boundaries. While effective under rigid constraints, these mechanisms lack the dynamic adaptability required for real-time environmental perturbations in complex non-homogeneous optical media.

To overcome these limitations, the optical Kerr effect has been exploited to generate stable solitons [3], [4], [5], [6], [7]. While the nonlinear Schrödinger equation and the standard Korteweg–de Vries (KdV) equation are widely employed as classical integrable models [8], [9], they overlook high-order spatial dispersion [10], [11]. Fifth-order KdV equations address this by introducing $\partial^{5}u/\partial x^{5}$ terms. Spectral methods achieve exponential convergence for smooth problems [12], [13] with far fewer degrees of freedom than the finite-difference time-domain method or the finite element method.

The central innovation of this study is a hybrid shifted Legendre–Laguerre operational matrix approach combining two orthogonal polynomial families into a single tensor-product framework, extending established operational-matrix spectral techniques [14], [15], [16], [17], [18], [19], [20]. This study is structured as follows: Section 2 develops the mathematical framework. Section 3 describes the numerical solution procedure. Section 4 presents results and validation. Section 5 concludes.

This study considers the initial–boundary value problem governed by the variable-coefficient fifth-order KdV equation in Eq. (1), together with the initial condition in Eq. (2) and the boundary conditions in Eq. (3):

where, $A$ = 2.0, $x_{0}$ = 5, $w$ = 1.2,

Since the spatial operator is of fifth order, Eq.(3) alone does not fully determine the problem. In the numerical scheme, the boundary closure is completed by additionally enforcing vanishing first spatial derivatives and, at one boundary, a vanishing second spatial derivative, consistent with the rapidly decaying pulse tails. This treatment provides a well-posed discrete formulation on the bounded spatial domain.

The variable coefficients are defined as follows:

The domain is defined by $L=10~\mathrm{mm}$ and $t\in[ 0,5]$. The equation couples nonlinear advection ($\alpha u\partial u/\partial x$), cubic dispersion ($\partial^3u/\partial x^3$), and quintic dispersion ($\partial^5 u/\partial x^5$).

For the finite spatial domain $[ 0,L]$, this study employs shifted Legendre polynomials [15], [20]:

The orthogonality property is as follows:

The corresponding recurrence relation is given by:

The spatial approximation is represented as:

The spatial derivatives are expressed as:

where, operational matrix $D_{x}$ is sparse and upper-triangular.

For the semi-infinite temporal domain $[ 0,\infty)$, classical Laguerre Polynomials satisfy the following, which is the natural basis for spectral discretization on semi-infinite domains [21], [22],

The recurrence relation is as follows:

The temporal approximation is represented as:

The temporal derivatives are expressed as:

where, $(D_t)_{ji}=-1$ if $i>j$; otherwise, $(D_t)_{ji}=0$.

The combined spatio-temporal approximation is represented as:

The corresponding Kronecker-product form is:

where, $\mathbf{U}\in\mathbb{R}^{mn}$ is the vectorized coefficient matrix. For $m=10$ and $n=8$, the resulting system contains $mn=80$ unknowns.

The extended operational matrices are defined as:

Substitution of approximation into the partial differential equation at Chebyshev-Gauss-Lobatto collocation points yields [13]:

Use of the method proposed by Newton leads to:

where, $\mathrm{J}(\mathrm{U})=\partial \mathrm{F} / \partial \mathrm{U}$ is the Jacobian. Convergence is realized when $\|\mathrm{F}(\mathrm{U})\|<10^{-12}$.

The robust solver using the adaptive trust-region strategy is as follows [23]:

• Quadratic model: \(m_k(s)=F(U^{(k)})+J(U^{(k)})^{T}s+0.5s^{T}B_k s\)

• Step accepted if actual reduction \(\rho_k>0.01\), trust region \(\Delta_k\) expanded by 2

• Otherwise, \(\Delta_k\) contracted by 0.25

The parameters are as follows:

• Maximum iterations: 100,000

• Tolerance: \(\|F(U)\|<10^{-12}\)

• Convergence: 15–25 iterations from zero initial guess

The convergence pattern is: $4.2 \rightarrow 2.1 \rightarrow 8.7 \times 10^{-1} \rightarrow 1.2 \times 10^{-3} \rightarrow 3.4 \times 10^{-9}$ (quadratic in final iterations).

The local convergence governed by Jacobian spectral radius is given by:

By the Banach fixed-point theorem, \(\rho(J)<1\) guarantees convergence of the fixed-point iteration; the quadratic rate observed near the solution stems from the Newton/trust-region update rather than from \(\rho(J)\) itself. The condition number is \(\kappa(J)\approx 28\), indicating excellent numerical conditioning. No iterative refinement is needed.

The relative Cauchy error is as follows:

This exhibits exponential convergence [12]. Table 1 summarizes the convergence results for the constant-coefficient benchmark at $t=0.1$.

The convergence behaviour of the spectral discretization is summarized in Table 1. As the number of degrees of freedom increases, the $L^\infty$ Cauchy error decreases consistently, while the $L^2$ quantity remains conserved to machine precision after the pulse is sufficiently resolved.

Figure 1 compares the convergence of the proposed spectral discretization with that of a fourth-order finite-difference scheme for the constant-coefficient benchmark at $t=0.1$. The vertical axis represents the $L^\infty$ error on a logarithmic scale, while the horizontal axis denotes the number of degrees of freedom. The spectral error falls by several orders of magnitude over a narrow range of degrees of freedom, whereas the finite-difference error decreases algebraically, so that the spectral method reaches a given accuracy with roughly an order of magnitude fewer unknowns. At the finest resolutions the spectral error approaches machine precision, so further refinement yields diminishing returns.

The steepness of the curve illustrates the exponential convergence rate characteristic of spectral methods applied to smooth solutions. Each additional polynomial degree in both $m$ and $n$ contributes approximately 1.5–2.0 decimal digits of accuracy, a stark contrast to finite difference methods which gain only $O(1/m^4)$ convergence. The curve plateaus slightly above (12,10) due to accumulated roundoff error at machine precision, indicating that further truncation increases provide diminishing returns.

This section presents a controlled comparison (Table 2) with a constant-coefficient fifth-order KdV test case ($\alpha_0=1.5$, $\beta_0=0.1$, $\gamma=0.01$, no spatial modulation) at time $t=0.1$, before the gradient catastrophe ($t^\ast\approx0.5$) that transforms this solution into a dispersive shock. A spectrally converged solution at machine precision is used as the reference.

The hybrid spectral method achieves high efficiency through exponential convergence. For the same target accuracy, the fourth-order finite-difference scheme requires roughly an order of magnitude more unknowns than the spectral discretization, while its error is considerably larger at equal degrees of freedom, reflecting algebraic versus exponential convergence. Under the identical integrator, the CPU times are comparable across methods (within a few seconds); the advantage lies in the smaller number of unknowns required to reach a given accuracy.

The three-dimensional surface plot in Figure 2 visualizes the spatio-temporal evolution of the optical field $u(x,t)$ for the polynomial pair $(m,n)=(8,6)$. The horizontal axes represent the spatial position $x\in[ 0,10]$ and time $t\in[ 0,5]$, while the vertical axis indicates the field amplitude $u(x,t)\in[-0.5,2.2]$. Under the variable coefficients, the pulse drifts in the positive $x$ direction and its amplitude gradually decreases due to dispersive radiation. Therefore, the solution does not remain a stationary trapped soliton, while the numerical scheme still resolves the evolution stably.

A critical test of numerical method correctness is conservation of physical invariants, a property emphasized in structure-preserving spectral schemes for KdV-type equations [24]. For the variable-coefficient fifth-order KdV equation, the $L^2$ norm (monitored here as a numerical diagnostic; for the variable-coefficient equation it is not a rigorously proven physical invariant) (under periodic or vanishing boundary conditions) is:

The computed solution exhibits a normalized $L^2$ quantity $E(t)/E(0)$ that decreases to about 0.18 by $t=5$, as amplitude is lost to dispersive radiation. This decrease is a property of the solution under the variable coefficients, not an artefact of the discretization; the remaining numerical error sources are:

• Exponential truncation error from finite polynomial degree ($\sim10^{-7}$ per operator application)

• Gaussian quadrature error in computing $\int u^2\,\mathrm{d}x$

• Accumulated discretization error over entire time interval

The decrease is smooth and monotone, without oscillatory growth or blow-up, which indicates that it reflects genuine radiation loss rather than a numerical instability or artificial damping introduced by the scheme.

Figure 3 displays the normalized $L^2$ quantity $E(t)/E(0)$ as a function of time $t\in[ 0,5]$. The ratio decreases from 1.0 to about 0.18 over this interval, reflecting physical energy loss caused by dispersive radiation under the variable coefficients rather than numerical instability.

The following observations are made:

(i) The normalized $L^2$ quantity decreases monotonically from 1.0 to about 0.18 over $t\in[ 0,5]$. This decline is physical: under the variable coefficients the pulse sheds energy to dispersive radiation as it evolves. The absence of oscillatory growth or blow-up shows that the decrease is genuine radiation loss rather than a numerical instability, and that the scheme remains stable throughout.

(ii) No spurious energy generation: there is no upward trend in $E(t)/E(0)$; the monotone decrease is physical and the scheme does not artificially create energy.

(iii) Smoothness: the curve decreases smoothly, without high-frequency oscillations that would signal an unresolved discretization or the onset of instability.

(iv) Amplitude decay: the normalized $L^2$ quantity falls steadily throughout the integration, showing that error does not accumulate.

The behaviour of $E(t)/E(0)$ therefore serves as a diagnostic confirming that the discretization is stable and resolves the physical energy loss, rather than as evidence of an exactly conserved invariant.

The two-dimensional contour heatmap in Figure 4 provides a top-down view of the spatio-temporal evolution, with spatial coordinate $x$ (horizontal axis) and time $t$ (vertical axis). The colour intensity encodes field amplitude via a Viridis colormap. The high-amplitude region does not form a stationary vertical stripe: it migrates toward larger $x$ as time increases and its intensity fades, so the pulse drifts and loses amplitude rather than remaining confined near $x \approx 5$.

The critical observations are as follows:

(i) The pulse exhibits a non-zero drift velocity and lateral migration; it does not remain fixed at $x \approx 5$.

(ii) The contour shows a pulse that migrates in $x$ and fades in intensity with time, rather than a stationary stripe: the peak drifts toward larger $x$, the pulse broadens as it disperses, and the peak amplitude decreases, consistent with radiation loss. The spectral scheme resolves this evolution without spurious oscillations.

(iii) Fading intensity: the colour within the high-amplitude region becomes progressively weaker with time, signifying a decreasing peak amplitude consistent with radiation loss.

(iv) Edge structure: the boundary between the pulse and the background broadens with time as dispersive radiation is shed; the spectral scheme resolves this structure without spurious oscillations.

Taken together, the contour data show drift and amplitude decay under the selected variable coefficients; they do not support a claim of stationary pulse confinement, and no claim of physical light trapping in a real photonic device is made here.

Cross-sectional profiles in Figure 5 show the spatial evolution of the optical field $u(x)$ at four different times: $t=0$, $t=1$, $t=2$, and $t=4$. The profiles demonstrate that the pulse shifts toward larger $x$ values while its peak amplitude decreases and its width increases over time. The deformation of the waveform indicates dispersive radiation and evolution under variable coefficients rather than a stationary trapped state.

The quantitative diagnostics are as follows:

• Peak amplitude (maximum of $u$): decreases from $\approx2.0$ at $t=0$ to $\approx0.7$ at $t=5$, a reduction of about 65%.

• Half-width at half-maximum: increases with time as the pulse broadens.

Quantitative diagnostics (peak amplitude, half-width, and centre-of-mass position) all change substantially over $t\in[ 0,5]$, confirming drift and decay rather than trapping.

The four profiles show a pulse that drifts and decays: the peak amplitude falls from $\approx2.0$ at $t=0$ to $\approx0.7$ at $t=5$; the profile broadens and develops a trailing oscillatory tail, departing from the initial $\mathrm{sech}^{2}$ shape; the peak position drifts from $x\approx5$ toward $x\approx8.9$; and the half-width grows as the pulse disperses. Under these variable coefficients the pulse is therefore neither trapped nor amplitude-preserving; the contribution of the study is the stable numerical resolution of this evolution.

(i) Amplitude decays: the peak falls from $\approx2.0$ at $t=0$ to $\approx0.7$ at $t=5$, indicating radiation loss rather than conservation.

(ii) Shape evolves: the profile broadens and develops a trailing oscillatory tail, departing from the initial $\mathrm{sech}^{2}$ form.

(iv) Spreading occurs: the half-width grows with time as the pulse disperses.

The most direct evidence that, under these variable coefficients, the pulse drifts and decays rather than being trapped; the value of the study is the stable numerical resolution of this evolution, not a demonstration of physical light trapping.

In engineered photonic waveguides with spatially periodic refractive index profiles [25], $n(x)=n_0+\Delta n\sin(2\pi x/\Lambda)$, the propagation coefficients $\alpha(x)$ and $\beta(x)$ vary spatially, producing a spatially modulated potential that shapes the propagation of the pulse. The proposed variable-coefficient model captures this dynamics with high accuracy using only a small number of spectral unknowns, enabling rapid numerical simulation and parameter studies.

Practical device applications include:

• Silicon-on-insulator waveguides with engineered sidewall roughness

• Photonic crystal microcavities with $Q$-factors $>10^4$

• Nonlinear metamaterials with spatially graded effective parameters

• Microresonators for quantum frequency comb generation [2], [5]

The small number of unknowns required by the spectral discretization makes parameter sweeps, sensitivity analysis, and inverse-design optimization practical on standard workstations.

This study developed and validated a hybrid shifted Legendre–classical Laguerre operational matrix method for the variable-coefficient fifth-order KdV equation, providing an accurate and computationally efficient spectral solver for variable-coefficient optical soliton dynamics.

The key achievements are as follows:

(i) Exponential convergence: the Cauchy error decreases rapidly with the polynomial degrees, reaching approximately $9.9\times10^{-6}$ at $(m,n)=(12,10)$ with 120 degrees of freedom.

(ii) Efficiency: the method reaches a given accuracy with about an order of magnitude fewer unknowns than fourth-order finite differences, and is orders of magnitude more accurate at equal degrees of freedom in the smooth-solution regime.

(iii) Physical behaviour: under the variable coefficients, the pulse drifts (peak from $x\approx5$ to $x\approx8.9$) and decays (amplitude $\approx2.0$ to $\approx0.7$), with the monitored $L^2$ quantity falling to $\approx0.18$ by $t=5$; the scheme resolves this evolution stably.

(iv) Spectral stability: $\rho(J)<1$ guarantees local convergence by Banach theorem.

(v) Computational cost: at matched accuracy the spectral solver is also faster in wall-clock time (96 versus 1024 unknowns; $\approx0.26$ s versus $\approx0.67$ s under the same time integrator).

The methodology therefore provides a stable and efficient framework for variable-coefficient soliton partial differential equations relevant to photonics. Demonstrating practical optical confinement in real-integrated-photonic systems lies beyond the scope of the present numerical study.