Inhomogeneous Coefficient Equations

• Inhomogeneous coefficient equations are differential, integro-differential, stochastic, or difference equations with coefficients that depend on time, space, state, or spectral variables, altering their canonical structure.
• Transformations like the Liouville and phase-function methods reallocate variable inhomogeneity to reveal simplified, tractable forms, enabling sharper energy estimates and stability analysis.
• These equations are pivotal in modeling complex physical phenomena in wave propagation, quantum systems, fluid dynamics, and solid-state physics, offering insights into medium-induced effects.

An inhomogeneous coefficient equation is a differential, integro-differential, stochastic, or difference equation in which one or more coefficients depend nontrivially on time, space, state, or spectral variables, so that the governing operator departs from a constant-coefficient form. In the cited literature, the term covers second-order wave equations with singular time-dependent coefficients, second-order inhomogeneous linear ordinary differential equations with variable coefficients, spatial nonlinear Schrödinger equations with topography-dependent coefficients, kinetic and phase-field models with spatially varying transport or potential weights, and discrete recurrences with variable stencil coefficients (Discacciati et al., 2021, Serkh et al., 2022, Karjanto et al., 2017, Tristani, 2013, Elbar et al., 17 Feb 2025, Mane, 25 Mar 2025). In one recent solid-state setting, “Inhomogeneous Coefficient Equation” is also the explicit name of a $k$-space amplitude equation for harmonic radiation under spatially inhomogeneous driving fields (Zhang et al., 1 Aug 2025). This suggests a structural class rather than a single canonical equation.

1. Scope and representative forms

Across the literature, the class is defined less by a single normal form than by the way coefficient inhomogeneity enters the principal part, lower-order terms, source coupling, or constitutive law. Representative examples include hyperbolic, parabolic, dispersive, kinetic, stochastic, and discrete models.

Setting	Representative equation	Inhomogeneity source
Wave equation	$$\partial_t^2 u(t,x)-\sum_{i=1}^n a_i(t)\partial_{x_i}^2u+l(t,\partial_t,\partial_x)u=f(t,x)$$	singular time-dependent $a_i,c_i,d,e$
Second-order ODE	$y''(x)+p(x)y'(x)+q(x)y(x)=r(x)$	variable $p,q,r$
Difference equation	$c_0(n)f(n)+c_1(n)f(n-1)+\dots+c_d(n)f(n-d)=r(n)$	variable coefficients $c_i(n)$
Spatial NLS	$$i\left(\partial_xA+c_g(x)^{-1}\partial_tA+\mu(x)A\right)+\alpha(x)\partial_t^2A+\beta(x)A|A|^2=0$$	topography-induced $c_g,\mu,\alpha,\beta$
Anisotropic Cahn–Hilliard	$$\partial_tu=\operatorname{Div}(u\nabla\mu),\ \mu=-\operatorname{Div}(A_p(x,\nabla u))+\frac1\varepsilon K(x)W'(u)$$	anisotropy $$\partial_t^2 u(t,x)-\sum_{i=1}^n a_i(t)\partial_{x_i}^2u+l(t,\partial_t,\partial_x)u=f(t,x)$$0 and inhomogeneous $$\partial_t^2 u(t,x)-\sum_{i=1}^n a_i(t)\partial_{x_i}^2u+l(t,\partial_t,\partial_x)u=f(t,x)$$1

The wave problem in $$\partial_t^2 u(t,x)-\sum_{i=1}^n a_i(t)\partial_{x_i}^2u+l(t,\partial_t,\partial_x)u=f(t,x)$$2 studied in (Discacciati et al., 2021) allows principal coefficients $$\partial_t^2 u(t,x)-\sum_{i=1}^n a_i(t)\partial_{x_i}^2u+l(t,\partial_t,\partial_x)u=f(t,x)$$3 and lower-order coefficients that may be compactly supported distributions in time, including Heaviside jumps and delta-type singularities. The ODE setting in (Serkh et al., 2022) treats “slowly varying coefficients” $$\partial_t^2 u(t,x)-\sum_{i=1}^n a_i(t)\partial_{x_i}^2u+l(t,\partial_t,\partial_x)u=f(t,x)$$4 on an interval, even when the magnitude of $$\partial_t^2 u(t,x)-\sum_{i=1}^n a_i(t)\partial_{x_i}^2u+l(t,\partial_t,\partial_x)u=f(t,x)$$5 is large. The discrete analogue in (Mane, 25 Mar 2025) treats a $$\partial_t^2 u(t,x)-\sum_{i=1}^n a_i(t)\partial_{x_i}^2u+l(t,\partial_t,\partial_x)u=f(t,x)$$6th-order backward recurrence on $$\partial_t^2 u(t,x)-\sum_{i=1}^n a_i(t)\partial_{x_i}^2u+l(t,\partial_t,\partial_x)u=f(t,x)$$7 with nonconstant stencil coefficients. The spatially inhomogeneous NLS in (Karjanto et al., 2017) derives variable dispersion and cubic nonlinearity from slowly varying depth $$\partial_t^2 u(t,x)-\sum_{i=1}^n a_i(t)\partial_{x_i}^2u+l(t,\partial_t,\partial_x)u=f(t,x)$$8 through the local dispersion relation $$\partial_t^2 u(t,x)-\sum_{i=1}^n a_i(t)\partial_{x_i}^2u+l(t,\partial_t,\partial_x)u=f(t,x)$$9. The phase-field model in (Elbar et al., 17 Feb 2025) places the inhomogeneity in the double-well coefficient $a_i,c_i,d,e$0, so that $a_i,c_i,d,e$1.

A broader consequence is that “inhomogeneous coefficient” can refer either to an externally prescribed medium profile, as in $a_i,c_i,d,e$2 or $a_i,c_i,d,e$3, or to an effective coefficient generated by reduction, homogenization, or gauge transformation. The solid-state ICE of (Zhang et al., 1 Aug 2025) makes this explicit by writing field inhomogeneity directly as coefficient terms in the amplitude evolution.

2. Transformations and canonical reductions

A recurring strategy is to transform an inhomogeneous-coefficient equation into a form with simpler principal structure while keeping the coefficient dependence in auxiliary functions. In second-order ODEs, the Liouville transformation

$a_i,c_i,d,e$4

eliminates the first-derivative term and yields

$a_i,c_i,d,e$5

thereby isolating a self-adjoint homogeneous part $a_i,c_i,d,e$6 (Serkh et al., 2022).

In the spatially inhomogeneous NLS for surface gravity waves, the convective and dissipative terms are removed by the co-moving time

$a_i,c_i,d,e$7

together with the amplitude scaling by $a_i,c_i,d,e$8. The dissipative coefficient is

$a_i,c_i,d,e$9

and the transformed equation becomes

$y''(x)+p(x)y'(x)+q(x)y(x)=r(x)$0

with the same physical content up to phase conjugation (Karjanto et al., 2017).

Variable-coefficient integrable models admit more elaborate reductions. The inhomogeneous variable-coefficient Hirota equation is mapped to the constant-coefficient Hirota equation through

$y''(x)+p(x)y'(x)+q(x)y(x)=r(x)$1

with explicit $y''(x)+p(x)y'(x)+q(x)y(x)=r(x)$2, $y''(x)+p(x)y'(x)+q(x)y(x)=r(x)$3, and gauge phase $y''(x)+p(x)y'(x)+q(x)y(x)=r(x)$4 chosen so that the variable-coefficient equation inherits the RH and Darboux structure of the homogeneous model (Zhou et al., 2022).

For inhomogeneous diffusion and Fokker–Planck models, (Costa et al., 2018) introduces the deformed measure and derivative

$y''(x)+p(x)y'(x)+q(x)y(x)=r(x)$5

so that the inhomogeneous FPE becomes

$y''(x)+p(x)y'(x)+q(x)y(x)=r(x)$6

The effect is to replace a position-dependent diffusion coefficient in $y''(x)+p(x)y'(x)+q(x)y(x)=r(x)$7 by constant diffusion in a deformed coordinate $y''(x)+p(x)y'(x)+q(x)y(x)=r(x)$8.

The hyperbolic theory in (Discacciati et al., 2021) uses a different reduction: after Fourier transform in $y''(x)+p(x)y'(x)+q(x)y(x)=r(x)$9, the second-order wave equation is rewritten as a first-order system

$p,q,r$0

which is then analyzed by quasi-symmetriser methods. In each case, transformation is not merely a change of variables; it is the mechanism by which coefficient inhomogeneity is reallocated to a form compatible with energy estimates, oscillatory quadrature, RH factorization, or deformed conservation laws.

3. Existence theories, weak formulations, and estimates

The analytical treatment depends on how singular the coefficients are. For the inhomogeneous wave equation with singular time-dependent coefficients, (Discacciati et al., 2021) constructs very weak solutions in the sense of Garetto–Ruzhansky by regularizing coefficients and data with mollifiers and defining the solution as a moderate net of classical solutions. For distributional coefficients $p,q,r$1 with $p,q,r$2, existence of a very weak solution of order $p,q,r$3 is obtained under $p,q,r$4-level Levi conditions with

$p,q,r$5

and logarithmic regularization scales. When the coefficients are classical, the very weak solution converges to the unique classical solution, so the theory is both stable and consistent (Discacciati et al., 2021).

In kinetic theory, the spatially inhomogeneous diffusively driven inelastic Boltzmann equation

$p,q,r$6

admits global solutions in the close-to-equilibrium regime, together with exponential convergence to the stationary homogeneous equilibrium $p,q,r$7. The linearized operator has a spectral gap in weighted Sobolev spaces, with

$p,q,r$8

and the nonlinear theory closes by new bilinear estimates on the collision operator (Tristani, 2013).

The anisotropic Cahn–Hilliard equation with disparate mobility and inhomogeneous potential coefficient $p,q,r$9 has global weak solutions on $c_0(n)f(n)+c_1(n)f(n-1)+\dots+c_d(n)f(n-d)=r(n)$0, $c_0(n)f(n)+c_1(n)f(n-1)+\dots+c_d(n)f(n-d)=r(n)$1, under strong convexity of $c_0(n)f(n)+c_1(n)f(n-1)+\dots+c_d(n)f(n-d)=r(n)$2 and $c_0(n)f(n)+c_1(n)f(n-1)+\dots+c_d(n)f(n-d)=r(n)$3 with $c_0(n)f(n)+c_1(n)f(n-1)+\dots+c_d(n)f(n-d)=r(n)$4. Under an energy convergence assumption, weak solutions converge as $c_0(n)f(n)+c_1(n)f(n-1)+\dots+c_d(n)f(n-d)=r(n)$5 to BV solutions of a weighted anisotropic Hele–Shaw flow, with interfacial energy

$c_0(n)f(n)+c_1(n)f(n-1)+\dots+c_d(n)f(n-d)=r(n)$6

and Gibbs–Thomson-type boundary condition

$c_0(n)f(n)+c_1(n)f(n-1)+\dots+c_d(n)f(n-d)=r(n)$7

(Elbar et al., 17 Feb 2025).

For the two-dimensional stationary incompressible inhomogeneous Navier–Stokes system with variable viscosity coefficient $c_0(n)f(n)+c_1(n)f(n-1)+\dots+c_d(n)f(n-d)=r(n)$8, (He et al., 2020) proves existence of weak solutions by passing to a fourth-order nonlinear elliptic equation for the stream function $c_0(n)f(n)+c_1(n)f(n-1)+\dots+c_d(n)f(n-d)=r(n)$9, with

$c_i(n)$0

The density and viscosity may have large variations, including piecewise-constant viscosity. The solution framework is $c_i(n)$1, $c_i(n)$2, with no smallness assumption on the coefficient variation.

Dispersive equations exhibit a different well-posedness pattern. For the inhomogeneous nonlinear biharmonic Schrödinger equation

$c_i(n)$3

with $c_i(n)$4, $c_i(n)$5, and $c_i(n)$6, local well-posedness is proved in the whole $c_i(n)$7-subcritical range $c_i(n)$8 under

$c_i(n)$9

and global existence is obtained in the $$i\left(\partial_xA+c_g(x)^{-1}\partial_tA+\mu(x)A\right)+\alpha(x)\partial_t^2A+\beta(x)A|A|^2=0$$0-subcritical regime under structural assumptions on $$i\left(\partial_xA+c_g(x)^{-1}\partial_tA+\mu(x)A\right)+\alpha(x)\partial_t^2A+\beta(x)A|A|^2=0$$1 and $$i\left(\partial_xA+c_g(x)^{-1}\partial_tA+\mu(x)A\right)+\alpha(x)\partial_t^2A+\beta(x)A|A|^2=0$$2 (Liu et al., 2021).

Linear parabolic theory supplies scale-consistent mixed-norm estimates rather than new weak solution notions. For

$$i\left(\partial_xA+c_g(x)^{-1}\partial_tA+\mu(x)A\right)+\alpha(x)\partial_t^2A+\beta(x)A|A|^2=0$$3

the homogeneous and Duhamel components satisfy optimal anisotropic mixed Lebesgue estimates, with admissibility relations such as

$$i\left(\partial_xA+c_g(x)^{-1}\partial_tA+\mu(x)A\right)+\alpha(x)\partial_t^2A+\beta(x)A|A|^2=0$$4

for the inhomogeneous term (Ostrovsky et al., 2014). This establishes a quantitative template for coefficient-sensitive well-posedness and regularity.

4. Structural effects of coefficient inhomogeneity

Coefficient inhomogeneity changes the equation not only analytically but dynamically. In the singular-coefficient wave equation, the lower-order terms are controlled by Levi-type bounds

$$i\left(\partial_xA+c_g(x)^{-1}\partial_tA+\mu(x)A\right)+\alpha(x)\partial_t^2A+\beta(x)A|A|^2=0$$5

and the quasi-symmetriser estimate succeeds precisely because the principal part remains nonnegative and the Kinoshita–Spagnolo condition is trivially satisfied for $$i\left(\partial_xA+c_g(x)^{-1}\partial_tA+\mu(x)A\right)+\alpha(x)\partial_t^2A+\beta(x)A|A|^2=0$$6 (Discacciati et al., 2021). Here coefficient regularity, multiplicity, and lower-order growth directly determine whether energy moderateness survives regularization.

In the spatially inhomogeneous NLS for water waves, the sign of the product $$i\left(\partial_xA+c_g(x)^{-1}\partial_tA+\mu(x)A\right)+\alpha(x)\partial_t^2A+\beta(x)A|A|^2=0$$7 governs focusing or defocusing behavior: $$i\left(\partial_xA+c_g(x)^{-1}\partial_tA+\mu(x)A\right)+\alpha(x)\partial_t^2A+\beta(x)A|A|^2=0$$8 Using the corrected nonlinear coefficient, the sign change occurs when $$i\left(\partial_xA+c_g(x)^{-1}\partial_tA+\mu(x)A\right)+\alpha(x)\partial_t^2A+\beta(x)A|A|^2=0$$9 crosses the Benjamin–Feir threshold

$c_g,\mu,\alpha,\beta$0

For the example $c_g,\mu,\alpha,\beta$1 with $c_g,\mu,\alpha,\beta$2 and $c_g,\mu,\alpha,\beta$3, the equation alternates between focusing and defocusing along $c_g,\mu,\alpha,\beta$4, which prevents sustained classical bright or dark solitons over extended intervals (Karjanto et al., 2017).

The solid-state ICE makes the symmetry-breaking effect of coefficient inhomogeneity particularly explicit. For a driving field

$c_g,\mu,\alpha,\beta$5

the $c_g,\mu,\alpha,\beta$6-space amplitude equation acquires the additional terms

$c_g,\mu,\alpha,\beta$7

which are absent in the standard homogeneous-field semiconductor Bloch equation. In centrosymmetric solids, these $c_g,\mu,\alpha,\beta$8-terms break the cancellation that suppresses even-order harmonics, and even-order harmonic radiation grows as the field inhomogeneity increases; the second-harmonic intensity follows the perturbative $c_g,\mu,\alpha,\beta$9 law with respect to incident intensity (Zhang et al., 1 Aug 2025).

In the focusing inhomogeneous mass-critical half-wave equation

$$\partial_tu=\operatorname{Div}(u\nabla\mu),\ \mu=-\operatorname{Div}(A_p(x,\nabla u))+\frac1\varepsilon K(x)W'(u)$$0

the spatial coefficient $$\partial_tu=\operatorname{Div}(u\nabla\mu),\ \mu=-\operatorname{Div}(A_p(x,\nabla u))+\frac1\varepsilon K(x)W'(u)$$1 breaks the symmetry group of the homogeneous case $$\partial_tu=\operatorname{Div}(u\nabla\mu),\ \mu=-\operatorname{Div}(A_p(x,\nabla u))+\frac1\varepsilon K(x)W'(u)$$2 and destroys momentum conservation for general $$\partial_tu=\operatorname{Div}(u\nabla\mu),\ \mu=-\operatorname{Div}(A_p(x,\nabla u))+\frac1\varepsilon K(x)W'(u)$$3. Under the conditions $$\partial_tu=\operatorname{Div}(u\nabla\mu),\ \mu=-\operatorname{Div}(A_p(x,\nabla u))+\frac1\varepsilon K(x)W'(u)$$4, $$\partial_tu=\operatorname{Div}(u\nabla\mu),\ \mu=-\operatorname{Div}(A_p(x,\nabla u))+\frac1\varepsilon K(x)W'(u)$$5, $$\partial_tu=\operatorname{Div}(u\nabla\mu),\ \mu=-\operatorname{Div}(A_p(x,\nabla u))+\frac1\varepsilon K(x)W'(u)$$6 even, with $$\partial_tu=\operatorname{Div}(u\nabla\mu),\ \mu=-\operatorname{Div}(A_p(x,\nabla u))+\frac1\varepsilon K(x)W'(u)$$7, $$\partial_tu=\operatorname{Div}(u\nabla\mu),\ \mu=-\operatorname{Div}(A_p(x,\nabla u))+\frac1\varepsilon K(x)W'(u)$$8, and $$\partial_tu=\operatorname{Div}(u\nabla\mu),\ \mu=-\operatorname{Div}(A_p(x,\nabla u))+\frac1\varepsilon K(x)W'(u)$$9, blowup solutions still exist at ground state mass, with

$$\partial_t^2 u(t,x)-\sum_{i=1}^n a_i(t)\partial_{x_i}^2u+l(t,\partial_t,\partial_x)u=f(t,x)$$00

The leading blowup rate is preserved, but the coefficient inhomogeneity changes the modulation analysis and removes the homogeneous symmetry-based shortcuts (Li, 2022).

A related ODE phenomenon is the turning-point obstruction in phase-function methods: when $$\partial_t^2 u(t,x)-\sum_{i=1}^n a_i(t)\partial_{x_i}^2u+l(t,\partial_t,\partial_x)u=f(t,x)$$01 changes sign or vanishes, the nonoscillatory phase-function representation can become delicate or numerically unstable, especially if one homogeneous solution grows rapidly (Serkh et al., 2022). Across models, the decisive issue is not merely whether coefficients vary, but how their variation interacts with hyperbolicity, symmetry, coercivity, and transport.

5. Discrete and computational treatments

Computational methods for inhomogeneous coefficient equations tend to mirror the analytic reductions. In the ODE setting, phase-function methods combine the numerical computation of a nonoscillatory phase $$\partial_t^2 u(t,x)-\sum_{i=1}^n a_i(t)\partial_{x_i}^2u+l(t,\partial_t,\partial_x)u=f(t,x)$$02 for

$$\partial_t^2 u(t,x)-\sum_{i=1}^n a_i(t)\partial_{x_i}^2u+l(t,\partial_t,\partial_x)u=f(t,x)$$03

with an adaptive Levin method for the oscillatory integral

$$\partial_t^2 u(t,x)-\sum_{i=1}^n a_i(t)\partial_{x_i}^2u+l(t,\partial_t,\partial_x)u=f(t,x)$$04

The resulting algorithm works in time essentially independent of the magnitude of $$\partial_t^2 u(t,x)-\sum_{i=1}^n a_i(t)\partial_{x_i}^2u+l(t,\partial_t,\partial_x)u=f(t,x)$$05 when $$\partial_t^2 u(t,x)-\sum_{i=1}^n a_i(t)\partial_{x_i}^2u+l(t,\partial_t,\partial_x)u=f(t,x)$$06 is slowly varying, and the reported examples exhibit near machine-precision accuracy over $$\partial_t^2 u(t,x)-\sum_{i=1}^n a_i(t)\partial_{x_i}^2u+l(t,\partial_t,\partial_x)u=f(t,x)$$07 (Serkh et al., 2022).

For the 1-D wave equation in an inhomogeneous medium,

$$\partial_t^2 u(t,x)-\sum_{i=1}^n a_i(t)\partial_{x_i}^2u+l(t,\partial_t,\partial_x)u=f(t,x)$$08

the Krylov subspace spectral method achieves unconditional stability, spectral accuracy in space, and second-order accuracy in time in the case of constant $$\partial_t^2 u(t,x)-\sum_{i=1}^n a_i(t)\partial_{x_i}^2u+l(t,\partial_t,\partial_x)u=f(t,x)$$09 and bandlimited $$\partial_t^2 u(t,x)-\sum_{i=1}^n a_i(t)\partial_{x_i}^2u+l(t,\partial_t,\partial_x)u=f(t,x)$$10. The same paper gives the first stability analysis without the bandlimited assumption on $$\partial_t^2 u(t,x)-\sum_{i=1}^n a_i(t)\partial_{x_i}^2u+l(t,\partial_t,\partial_x)u=f(t,x)$$11, still for constant $$\partial_t^2 u(t,x)-\sum_{i=1}^n a_i(t)\partial_{x_i}^2u+l(t,\partial_t,\partial_x)u=f(t,x)$$12, and shows that variable $$\partial_t^2 u(t,x)-\sum_{i=1}^n a_i(t)\partial_{x_i}^2u+l(t,\partial_t,\partial_x)u=f(t,x)$$13 introduces $$\partial_t^2 u(t,x)-\sum_{i=1}^n a_i(t)\partial_{x_i}^2u+l(t,\partial_t,\partial_x)u=f(t,x)$$14-dependent growth that can impair unconditional stability (Rester et al., 2023).

The singular-coefficient wave theory of (Discacciati et al., 2021) is corroborated numerically by two toy models. For the Heaviside principal coefficient

$$\partial_t^2 u(t,x)-\sum_{i=1}^n a_i(t)\partial_{x_i}^2u+l(t,\partial_t,\partial_x)u=f(t,x)$$15

a Lax–Friedrichs scheme applied to the first-order system shows convergence of $$\partial_t^2 u(t,x)-\sum_{i=1}^n a_i(t)\partial_{x_i}^2u+l(t,\partial_t,\partial_x)u=f(t,x)$$16 to the piecewise distributional solution, independently of the mollifier. For the delta coefficient

$$\partial_t^2 u(t,x)-\sum_{i=1}^n a_i(t)\partial_{x_i}^2u+l(t,\partial_t,\partial_x)u=f(t,x)$$17

the relevant energy ratio is not bounded independently of $$\partial_t^2 u(t,x)-\sum_{i=1}^n a_i(t)\partial_{x_i}^2u+l(t,\partial_t,\partial_x)u=f(t,x)$$18, matching the theoretical failure of uniform energy bounds.

The discrete analogue of Green-function solution theory is developed in (Mane, 25 Mar 2025). For

$$\partial_t^2 u(t,x)-\sum_{i=1}^n a_i(t)\partial_{x_i}^2u+l(t,\partial_t,\partial_x)u=f(t,x)$$19

retarded and advanced Green’s functions are both required to construct the full particular solution over $$\partial_t^2 u(t,x)-\sum_{i=1}^n a_i(t)\partial_{x_i}^2u+l(t,\partial_t,\partial_x)u=f(t,x)$$20. The kernels are expressed by Casoratian ratios built from fundamental solutions of the homogeneous recurrence, making the method a direct discrete counterpart of variation of parameters.

These methods indicate that numerical treatment is most effective when coefficient inhomogeneity can be absorbed into spectral symbols, Green kernels, or nonoscillatory phases. Where the coefficient enters the principal part too singularly, as in the $$\partial_t^2 u(t,x)-\sum_{i=1}^n a_i(t)\partial_{x_i}^2u+l(t,\partial_t,\partial_x)u=f(t,x)$$21-wave toy model, numerical instability is not merely an implementation artifact but reflects a genuine analytical obstruction.

6. Physical realizations and model-specific phenomena

The principal applications in the cited literature are medium-induced wave propagation, transport in heterogeneous continua, and coefficient-driven pattern formation. In wave propagation with abrupt temporal changes, the singular-coefficient wave equation models media in which the principal part switches on or experiences an impulsive event in time; the Heaviside and delta toy models are explicitly interpreted in terms of interface and impulse effects (Discacciati et al., 2021).

In water-wave theory, the spatial inhomogeneous NLS is the envelope equation for weakly nonlinear narrow-band gravity waves over slowly varying topography. Because $$\partial_t^2 u(t,x)-\sum_{i=1}^n a_i(t)\partial_{x_i}^2u+l(t,\partial_t,\partial_x)u=f(t,x)$$22 is fixed while $$\partial_t^2 u(t,x)-\sum_{i=1}^n a_i(t)\partial_{x_i}^2u+l(t,\partial_t,\partial_x)u=f(t,x)$$23 adjusts through $$\partial_t^2 u(t,x)-\sum_{i=1}^n a_i(t)\partial_{x_i}^2u+l(t,\partial_t,\partial_x)u=f(t,x)$$24, the coefficients $$\partial_t^2 u(t,x)-\sum_{i=1}^n a_i(t)\partial_{x_i}^2u+l(t,\partial_t,\partial_x)u=f(t,x)$$25, $$\partial_t^2 u(t,x)-\sum_{i=1}^n a_i(t)\partial_{x_i}^2u+l(t,\partial_t,\partial_x)u=f(t,x)$$26, and $$\partial_t^2 u(t,x)-\sum_{i=1}^n a_i(t)\partial_{x_i}^2u+l(t,\partial_t,\partial_x)u=f(t,x)$$27 inherit the depth variation. The corrected coefficient formulas imply that realistic topographies may force repeated transitions between modulational instability and stability (Karjanto et al., 2017).

In phase-field dynamics, the inhomogeneous coefficient $$\partial_t^2 u(t,x)-\sum_{i=1}^n a_i(t)\partial_{x_i}^2u+l(t,\partial_t,\partial_x)u=f(t,x)$$28 in $$\partial_t^2 u(t,x)-\sum_{i=1}^n a_i(t)\partial_{x_i}^2u+l(t,\partial_t,\partial_x)u=f(t,x)$$29 becomes, in the sharp-interface limit, a weighted anisotropic perimeter and a drift term in the Gibbs–Thomson relation. The coefficient landscape therefore acts directly on interface pinning, acceleration, and reorientation (Elbar et al., 17 Feb 2025).

In solid-state high-harmonic generation, the ICE of (Zhang et al., 1 Aug 2025) is designed precisely for spatially inhomogeneous driving fields, where the semiconductor Bloch equation fails. Using graphene as the example, the model predicts increasing even-order harmonics with increasing field inhomogeneity, second-harmonic intensity consistent with perturbation theory, and wavelength-dependent separation of interband and intraband contributions.

Integrable variable-coefficient dispersive systems generate another class of applications. For the inhomogeneous variable-coefficient Hirota equation, the special transformation from the constant-coefficient Hirota model yields multi-solitons and high-order solitons whose trajectories and amplitudes are modulated by $$\partial_t^2 u(t,x)-\sum_{i=1}^n a_i(t)\partial_{x_i}^2u+l(t,\partial_t,\partial_x)u=f(t,x)$$30 and $$\partial_t^2 u(t,x)-\sum_{i=1}^n a_i(t)\partial_{x_i}^2u+l(t,\partial_t,\partial_x)u=f(t,x)$$31. Periodic coefficient choices produce new morphologies, including heart-shaped periodic waves and O-shaped periodic waves (Zhou et al., 2022).

Kinetic and fluid applications are equally coefficient-sensitive. The inelastic Boltzmann equation uses a restitution coefficient that may depend on impact velocity and incorporates weak spatial inhomogeneity through the transport term on $$\partial_t^2 u(t,x)-\sum_{i=1}^n a_i(t)\partial_{x_i}^2u+l(t,\partial_t,\partial_x)u=f(t,x)$$32 (Tristani, 2013). The stationary incompressible inhomogeneous Navier–Stokes system uses variable viscosity $$\partial_t^2 u(t,x)-\sum_{i=1}^n a_i(t)\partial_{x_i}^2u+l(t,\partial_t,\partial_x)u=f(t,x)$$33 and admits explicit parallel, concentric, and radial solutions when $$\partial_t^2 u(t,x)-\sum_{i=1}^n a_i(t)\partial_{x_i}^2u+l(t,\partial_t,\partial_x)u=f(t,x)$$34 is piecewise constant (He et al., 2020). A plausible implication is that coefficient inhomogeneity functions as a compact encoding of geometry, constitutive response, and external control within a single operator framework.

7. Limitations, controversies, and open directions

The main limitations are highly model-dependent, but several recurring obstructions appear. In the singular wave theory, the principal coefficients must satisfy $$\partial_t^2 u(t,x)-\sum_{i=1}^n a_i(t)\partial_{x_i}^2u+l(t,\partial_t,\partial_x)u=f(t,x)$$35, the KS structure must remain available, and the Levi constants must grow no faster than $$\partial_t^2 u(t,x)-\sum_{i=1}^n a_i(t)\partial_{x_i}^2u+l(t,\partial_t,\partial_x)u=f(t,x)$$36 at the regularized level; stronger singularities, or lack of control on the lower-order terms, can destroy moderateness or uniform energy bounds (Discacciati et al., 2021). Extending the analysis to higher-order systems, combined time-space singularities, or refined microlocal propagation remains open.

In the phase-function ODE framework, the case $$\partial_t^2 u(t,x)-\sum_{i=1}^n a_i(t)\partial_{x_i}^2u+l(t,\partial_t,\partial_x)u=f(t,x)$$37 or $$\partial_t^2 u(t,x)-\sum_{i=1}^n a_i(t)\partial_{x_i}^2u+l(t,\partial_t,\partial_x)u=f(t,x)$$38 is delicate because the phase representation may become ill-conditioned when one homogeneous solution grows rapidly (Serkh et al., 2022). In the KSS analysis, unconditional stability is not guaranteed for variable leading coefficient $$\partial_t^2 u(t,x)-\sum_{i=1}^n a_i(t)\partial_{x_i}^2u+l(t,\partial_t,\partial_x)u=f(t,x)$$39, and sharper $$\partial_t^2 u(t,x)-\sum_{i=1}^n a_i(t)\partial_{x_i}^2u+l(t,\partial_t,\partial_x)u=f(t,x)$$40-independent bounds remain open (Rester et al., 2023). For the anisotropic Cahn–Hilliard problem, uniqueness of weak solutions in the disparate mobility case is generally open, and a varifold-type sharp-interface limit without the energy convergence assumption is also open (Elbar et al., 17 Feb 2025). In the granular Boltzmann setting, weakly inhomogeneous far-from-equilibrium theory for velocity-dependent restitution is not covered (Tristani, 2013).

A notable controversy concerns stochastic models with position-dependent mass and damping. The deformed Fokker–Planck formulation of (Costa et al., 2018) rewrites inhomogeneous diffusion through deformed derivatives and a deformed entropy balance, whereas the microscopic Caldeira–Leggett derivation in (Łuczka, 2021) argues that the corresponding phenomenological Langevin and Fokker–Planck equations are not equivalent. The microscopic derivation yields multiplicative noise, damping

$$\partial_t^2 u(t,x)-\sum_{i=1}^n a_i(t)\partial_{x_i}^2u+l(t,\partial_t,\partial_x)u=f(t,x)$$41

the Stratonovich interpretation, and a stationary phase-space density

$$\partial_t^2 u(t,x)-\sum_{i=1}^n a_i(t)\partial_{x_i}^2u+l(t,\partial_t,\partial_x)u=f(t,x)$$42

while the overdamped equilibrium position density is mass-independent. The disagreement is therefore not terminological but structural: it concerns which coefficient dependence is compatible with fluctuation–dissipation and detailed balance (Łuczka, 2021, Costa et al., 2018).

Taken together, these results show that the study of inhomogeneous coefficient equations is organized around three questions: which coefficient variations preserve analytic control, which transformations convert them to tractable canonical forms, and which physical effects are created precisely because the coefficients are not homogeneous. The modern literature treats those questions across continuous and discrete models, from very weak hyperbolic solutions and weighted sharp-interface limits to Green-function recurrences, deformed stochastic dynamics, and symmetry-breaking harmonic radiation.