1 Introduction↩︎

Many physical, biological, and engineering systems are governed by partial differential equations (PDEs) whose coefficients vary on spatial or temporal scales much smaller than the scale at which the system is observed. Examples include flow through porous media, heat conduction in composite materials, elastic deformation of heterogeneous media, wave propagation in complex environments, and transport in random or disordered structures. Homogenization theory provides a systematic framework for replacing the original heterogeneous problem by an effective macroscopic model that captures the large-scale behavior of the solution [1]–[3]. They also guide the design of multiscale numerical methods for PDEs with highly oscillatory coefficients [4]–[6].

In the classical periodic setting, the microscopic medium is modeled by a repeating unit cell, allowing effective equations to be derived through asymptotic expansions, two-scale convergence, or variational methods. The homogenized coefficients describe the macroscopic response of the oscillatory medium. However, many realistic media are nonperiodic and exhibit disorder, randomness, intermittency, or long-range correlations. This motivates stochastic homogenization, where coefficients are modeled as stationary random fields rather than periodic functions [7]–[9], and studies when equations with random microscopic coefficients converge, at large scales, to deterministic effective equations. Under assumptions such as stationarity, ergodicity, coercivity, and suitable compactness or mixing conditions, solutions often converge to a deterministic homogenized limit. Effective coefficients are obtained not from periodic cell problems, but from corrector problems on probability spaces or large finite domains [8]–[10]. Nonlinearity further complicates the theory, since effective behavior is not determined by simple coefficient averaging. Oscillations in both the coefficients and solution gradients may interact, and weak convergence is generally insufficient for passing to nonlinear limits. Thus, nonlinear homogenization must account for microscopic oscillations, loss of compactness, and nonlinear material response [11]–[13].

Several tools support nonlinear and stochastic homogenization, including variational methods, \(\Gamma\)-convergence, compensated compactness, monotonicity methods, subadditive ergodic theorems, corrector estimates, concentration inequalities, and quantitative homogenization [7], [9], [10], [12], [14], [15]. In this work, we emphasize Young measures as a central framework for describing nonlinear oscillations and randomness-induced limiting distributions that are not captured by weak convergence alone but are essential for passing nonlinear quantities to the limit [16]–[18]. The goal is to connect variational structure, stochastic heterogeneity, Young measure descriptions, and numerical approximation in order to identify effective macroscopic laws without fully resolving every microscopic realization or scale. By retaining distributional information beyond averages, Young measure formulations provide a bridge between rigorous nonlinear and stochastic homogenization theory and computational approaches such as discretization, relaxation, linear programming, adaptive atomic approximations, and neural-network parameterizations [19], [20]. A key computational challenge is that Young measures are probability-valued maps encoding distributions over gradients, phases, microscopic states, and random environments, so their approximation is intrinsically high-dimensional across physical, phase, and probability spaces.

Quantum computing for differential equations aims to reduce the cost of high-dimensional PDE computation by exploiting techniques such as quantum linear algebra and Hamiltonian simulation. For linear PDEs, standard approaches either discretize the equation and apply quantum linear-system solvers, or embed the resulting possibly nonunitary dynamics into Hamiltonian simulation frameworks, for example through Schrödingerization [21]–[24]. For linear multiscale problems with scale separation, Hu et. al developed quantum algorithms for elliptic, parabolic, and hyperbolic PDEs based on either direct oscillatory models or homogenized equations, combined with higher-dimensional liftings and Schrödingerization [25]. Balazi et. al proposed a quantum-enhanced numerical homogenization method based on the Localized Orthogonal Decomposition, achieving logarithmic scaling in fine-scale resolution via quantum local solvers [26]. Nonlinear PDEs are more challenging because quantum evolution is intrinsically linear. This motivates a Young-measure formulation [27], in which oscillatory, singular, or uncertain nonlinear solutions are represented by probability-valued measures. Under this representation, certain nonlinear PDE problems can be reformulated as linear programming (LP) problems with linear objectives and constraints. Since Young measures are defined on high-dimensional physical, state, and possibly random-parameter domains, their full computation suffers from the curse of dimensionality. Quantum linear-programming methods may therefore offer polynomial advantages for recovering the full Young measure, although such advantages may be less relevant when only the expected weak solution is required [27]–[30].

In this paper, we first introduce nonlinear heterogeneous and random homogenization models and formulate them using Young measures. We then show how the measure-valued formulation transforms the nonlinear PDE into linear constraints involving normalization, barycenter consistency, compatibility, and weak force balance. Next, we describe the resulting discrete Young-measure linear programs for deterministic periodic and random media. We then discuss quantum algorithms for solving these LPs, focusing on the quantum central path method and comparing its complexity with that of direct classical solvers and classical Young-measure LP approaches.

The quantum advantage analysis is detailed in Sections 4.2 and 4.3. The baseline is a direct classical solver (DCS), which must resolve all \(\varepsilon\)-scale oscillations at cost \(\widetilde{O}(\varepsilon^{-d})\) with \(d\) being the spatial dimension. The Young-measure LP, by contrast, treats the microscale \(y\), the gradient \(\xi\), and the random environment \(\omega\) as independent discretization variables, and its size grows with the target accuracy \(\delta\) rather than with \(\varepsilon\). Under first-order algebraic discretization, which requires minimal regularity assumptions, two rigorous regimes of quantum advantage over DCS are identified (see Table 1):

1. Deterministic nonlinear homogenization. Recovering the homogenized solution \(u_0\) to accuracy \(\delta=\varepsilon^\alpha\), where \(\alpha\in(0,1)\) reflects the homogenization convergence rate, yields a polynomial QCP speedup over DCS whenever \(\alpha<2d/(3d+4)\).

2. Stochastic homogenization at fine-scale accuracy. With \(r\) random variables and \(N_\omega^r\) total random discretization points treated as a free parameter (independent of \(\varepsilon\)), QCP achieves a square-root reduction in stochastic sampling cost whenever \(N_\omega^r\gtrsim\varepsilon^{-(d+2)}\); the speedup grows with \(N_\omega^r\).

We further propose in Conjecture 1 that if the Young-measure density is analytic in \(x\) and \(y\), spectral discretization in those variables extends quantum advantage to fine-scale accuracy \(\delta=\varepsilon\) for all \(d\ge2\) in the deterministic setting, with an analogous improvement in the stochastic case.

Organizatin: The rest of the paper is organized as follows. Section 2 reviews the deterministic periodic and stochastic homogenization settings for nonlinear elliptic equations. Section 3 develops the Young-measure formulation and its finite-dimensional LP discretization for both periodic and random media. Section 4 discusses the quantum central path algorithm for the resulting Young-measure LPs, derives the associated query and gate complexity bounds, and compares them with direct classical multiscale solvers and classical LP methods. Section 5 presents numerical experiments for one- and two-dimensional deterministic and stochastic benchmark problems. Section 6 concludes with a summary of the proposed framework and outlines several directions for future work.

Notation: \(\widetilde{O}(f(n))\) denotes complexity scaling as \(O(f(n))\) up to logarithmic factors.

2 Preliminaries on elliptic homogenization↩︎

We first recall the homogenization framework used throughout the paper. The deterministic setting describes periodic nonlinear elliptic equations with rapidly oscillating coefficients, while the random setting describes media whose constitutive laws depend on finitely many probabilistic states. In both cases, the goal is to replace the heterogeneous problem by an effective deterministic equation whose coefficient, or constitutive law, captures the macroscopic response of the medium.

2.1 Nonlinear periodic homogenization↩︎

We first consider deterministic periodic homogenization for nonlinear monotone elliptic equations. Let \(\Omega\subset\mathbb{R}^d\) be a bounded domain and let \(Y=(0,1)^d\) denote the periodicity cell. For \(\varepsilon>0\), we study \[\label{eq:pde95eps} -\nabla\cdot a\!\left(\frac{x}{\varepsilon},\nabla u_\varepsilon(x)\right) = f(x) \quad\text{in }\Omega, \qquad u_\varepsilon=0 \quad\text{on }\partial\Omega,\tag{1}\] where \(a(y,\xi)\) is \(Y\)-periodic in \(y\). We assume that \(a\) satisfies the standard monotone \(L^p\) conditions: for some \(p>1\) and constants \(c_0,c_1>0\), \[\begin{align} a(y,\xi)\cdot \xi &\ge c_0|\xi|^p-c_1, \tag{2} \\ |a(y,\xi)| &\le c_1(1+|\xi|^{p-1}), \tag{3} \\ \big(a(y,\xi)-a(y,\eta)\big)\cdot(\xi-\eta) &\ge 0. \tag{4} \end{align}\] Under these assumptions, \((u_\varepsilon)\) is uniformly bounded in \(W_0^{1,p}(\Omega)\). Hence, along a subsequence, \[u_\varepsilon \rightharpoonup u_0 \quad\text{weakly in }W_0^{1,p}(\Omega).\] The homogenized limit \(u_0\) solves \[\label{eq:hom95pde} -\nabla\cdot a_{\mathrm{hom}}(\nabla u_0) = f \quad\text{in }\Omega, \qquad u_0=0 \quad\text{on }\partial\Omega.\tag{5}\]

The effective constitutive law \(a_{\mathrm{hom}}\) is defined through a periodic cell problem. For each \(\xi\in\mathbb{R}^d\), find \(w_\xi\in W^{1,p}_{\#}(Y)\), with \(\int_Y w_\xi\,dy=0\), such that \[\label{eq:cell95problem} -\nabla_y\cdot a\big(y,\xi+\nabla_y w_\xi(y)\big) = 0 \quad\text{in }Y.\tag{6}\] Then \[\label{eq:ahom95definition} a_{\mathrm{hom}}(\xi) = \int_Y a\big(y,\xi+\nabla_y w_\xi(y)\big)\,dy.\tag{7}\] Here \(W^{1,p}_{\#}(Y)\) denotes the Sobolev space of \(Y\)-periodic functions with \(Y\)-periodic weak derivatives. The zero-average condition fixes the additive constant in the cell corrector. This formulation is standard in periodic homogenization and two-scale convergence; see, for example, [1], [2], [31], [32].

2.2 Stochastic homogenization↩︎

We also consider homogenization in random media. Let \((\Theta,\mathcal{F},\mathbb{P})\) be a probability space equipped with a measure-preserving ergodic group action \((\tau_x)_{x\in\mathbb{R}^d}\). The constitutive law is described by a random monotone operator \[a:\Theta\times\mathbb{R}^d\to\mathbb{R}^d, \qquad (\omega,\xi)\mapsto a(\omega,\xi),\] which is assumed to be stationary in the sense that \[a(\tau_x\omega,\xi)\] represents the material response at the spatial position \(x\).

For \(\varepsilon>0\), we consider the stochastic boundary-value problem \[\label{eq:sto95pde} -\nabla\cdot a\!\left( \tau_{x/\varepsilon}\omega, \nabla u_\varepsilon(x,\omega) \right) = f(x) \quad\text{in }\Omega, \qquad u_\varepsilon=0 \quad\text{on }\partial\Omega.\tag{8}\] Here the factor \(x/\varepsilon\) describes rapidly oscillating random heterogeneities at the microscopic scale \(\varepsilon\).

Under the standard assumptions of stationarity, ergodicity, uniform monotonicity, and polynomial growth, the family \(u_\varepsilon(\cdot,\omega)\) converges, for \(\mathbb{P}\)-almost every \(\omega\), to a deterministic limit \(u_0\) solving the homogenized problem \[\label{eq:sto95hom95limit} -\nabla\cdot a_{\mathrm{hom}}(\nabla u_0) = f \quad\text{in }\Omega, \qquad u_0=0 \quad\text{on }\partial\Omega.\tag{9}\] The effective operator \(a_{\mathrm{hom}}\) is deterministic as a consequence of ergodicity; see [2], [7], [33], [34].

In the finite-state setting used later in this paper, the random parameter takes values \[\omega_i\in\Theta, \qquad i=1,\ldots,N,\] with probabilities \[p_i>0, \qquad \sum_{i=1}^N p_i=1.\] Each state \(\omega_i\) corresponds to a constitutive law \[a_i(\xi) := a(\omega_i,\xi).\] Accordingly, the random medium is represented by the finite ensemble \[\label{eq:finite95random95medium} \big\{(a_i,p_i)\big\}_{i=1}^N.\tag{10}\]

To approximate the stochastic homogenized operator, we introduce in the next section a Young-measure formulation together with a finite-dimensional linear-programming discretization adapted to the finite-state random setting.

3 Young-measure formulation and LP implementation↩︎

This section presents a measure-valued formulation for nonlinear homogenization problems in both periodic and stochastic settings, together with the finite-dimensional linear programming problems used in the numerical implementation. The central idea is to replace rapidly oscillatory gradients by parametrized Young measures that encode the local distribution of microscopic states [16]–[18]. In both the periodic and stochastic frameworks, the resulting constraints: normalization, barycenter consistency, weak equilibrium, and admissibility conditions, are linear in the measures. Consequently, after discretization of the spatial variables and gradient state space, the homogenization problem becomes a large but structured linear program. The formulation provides a unified framework for deterministic periodic media and stochastic heterogeneous materials while avoiding direct resolution of fine microscopic scales.

3.1 Nonlinear periodic homogenization↩︎

Let \(\Omega\subset\mathbb{R}^d\) be bounded and let \(Y=(0,1)^d\) denote the periodic reference cell. Consider the nonlinear elliptic problem \[-\nabla\cdot a\!\left(\frac{x}{\varepsilon},\nabla u^\varepsilon(x)\right)=f(x) \qquad \text{in }\Omega,\] subject to the boundary condition \[u^\varepsilon=g \qquad \text{on }\partial\Omega.\] Here \[a:Y\times\mathbb{R}^d\to\mathbb{R}^d\] is assumed to be \(Y\)-periodic in the first variable.

To describe microscopic oscillations in the gradients, we introduce a parametrized Young measure \[\mu_{x,y}\in\mathcal{P}(\mathbb{R}^d), \qquad (x,y)\in\Omega\times Y,\] which represents the local probability distribution of microscopic gradients at macroscopic position \(x\) and microscopic phase \(y\). The normalization condition is imposed pointwise: \[\int_{\mathbb{R}^d} d\mu_{x,y}(\xi)=1 \qquad \text{for a.e. }(x,y)\in\Omega\times Y.\]

The admissible unknowns are \[u\in W^{1,p}(\Omega), \qquad \mu_{x,y}\in\mathcal{P}(\mathbb{R}^d),\] subject to the following constraints: \[\begin{align} &\mu_{x,y}\ge0, \\ &\int_{\mathbb{R}^d} d\mu_{x,y}(\xi)=1 \qquad \text{for a.e. }(x,y), \\ &\nabla u(x) = \int_Y\int_{\mathbb{R}^d} \xi\, d\mu_{x,y}(\xi)\,dy \qquad \text{for a.e. }x, \\ &-\nabla\cdot \left( \int_Y \int_{\mathbb{R}^d} a(y,\xi)\, d\mu_{x,y}(\xi)\,dy \right) = f(x) \qquad \text{in } \Omega . \end{align}\]

Importantly, all terms involving \(\mu\) are linear in the measure. After discretization, the formulation therefore reduces to a linear programming problem in the discrete measure masses.

For presentation purposes, it is convenient to introduce the marginal measure \[\widehat\mu_x(\xi) := \int_Y \mu_{x,y}(\xi)\,dy,\] which represents the effective distribution of microscopic gradients at the macroscopic point \(x\).

3.2 Stochastic homogenization↩︎

We now consider the nonlinear random elliptic problem \[-\nabla\cdot a\!\left(\tau_{x/\varepsilon}\omega,\nabla u^\varepsilon(x)\right) = f(x) \qquad \text{in }\Omega,\] with boundary condition \[u^\varepsilon=g \qquad \text{on }\partial\Omega.\]

Here \((\Theta,\mathcal{F},\mathbb{P})\) is a probability space of environments, \[\tau_z:\Theta\to\Theta\] is a measure-preserving ergodic group action, and \[a:\Theta\times\mathbb{R}^d\to\mathbb{R}^d\] is the nonlinear random constitutive law.

For almost every macroscopic point \(x\in\Omega\), we introduce a probability measure \[\mu_x\in\mathcal{P}(\Theta\times\mathbb{R}^d),\] interpreted as the local joint distribution of the random environment and microscopic gradient state. The admissibility conditions become \[\begin{align} &\mu_x\ge0, \qquad \int_{\mathbb{R}^d} d\mu_x(\omega,\xi) = d\mathbb{P}(\omega) \quad \mathrm{for~a.e.}~x, \\ &\nabla u(x) = \int_{\Theta\times\mathbb{R}^d} \xi\,d\mu_x(\omega,\xi) \qquad \text{for a.e. }x, \\ &-\nabla\cdot \left( \int_{\Theta\times\mathbb{R}^d} a(\omega,\xi)\, d\mu_x(\omega,\xi) \right) = f \qquad \text{in } \Omega, \\ &u=g \qquad \text{on }\partial\Omega. \end{align}\]

As in the periodic case, the objective and all admissibility constraints are linear in the measures. For numerical implementation, the probability space is approximated by finitely many environments \[\omega_i, \qquad i=1,\dots,N_\omega,\] with associated probabilities \[p_i>0, \qquad \sum_{i=1}^{N_\omega}p_i=1.\] In contrast to the periodic setting, there is no explicit fast variable \(y\). The unknown measure family is therefore indexed only by the macroscopic point and the random state: \[\mu_{x,i}\ge0.\]

The constraints become \[\begin{align} &\int_{\mathbb{R}^d} d\mu_{x,i}(\xi)=p_i \qquad \text{for a.e. }x, \quad i=1,\dots,N_\omega, \\ &\nabla u(x) = \sum_{i=1}^{N_\omega} \int_{\mathbb{R}^d} \xi\,d\mu_{x,i}(\xi) \qquad \text{for a.e. }x, \\ &-\nabla\cdot \left( \sum_{i=1}^{N_\omega} \int_{\mathbb{R}^d} a_i(\xi)\, d\mu_{x,i}(\xi) \right) = f \qquad \text{in } \Omega, \end{align}\] where \[a_i(\xi)=a(\omega_i,\xi).\]

Equivalently, one may introduce normalized probability measures \[\nu_{x,i}:=\frac{\mu_{x,i}}{p_i} \in\mathcal{P}(\mathbb{R}^d),\] so that \[\mu_{x,i}=p_i\nu_{x,i}.\]

The global marginal distribution used for visualization is \[\widehat\mu_x(\xi) = \sum_{i=1}^{N_\omega}\mu_{x,i}(\xi) = \sum_{i=1}^{N_\omega} p_i\nu_{x,i}(\xi).\]

The discrete random objective functional becomes \[\mathcal{J}_{\mathrm{rand}}(u,\mu) = \int_\Omega \sum_{i=1}^{N_\omega} \int_{\mathbb{R}^d} W_i(\xi)\, d\mu_{x,i}(\xi)\,dx - \int_\Omega f(x)u(x)\,dx,\] where \[W_i(\xi)=W(\omega_i,\xi).\]

3.3 Discrete LP implementation↩︎

After discretization of the macroscopic variable and, when required, the microscopic phase, random state, and gradient variable, all unknowns are assembled into a vector \[z=[u;F], \qquad F\ge0,\] where \(F\) contains the discrete masses approximating the Young measures \(\mu_{x,y}\) in the periodic setting or \(\mu_{x,i}\) in the stochastic setting.

The resulting finite-dimensional optimization problem takes the form \[\min_z c^\top z \qquad \text{subject to} \qquad A_{\mathrm{eq}}z=b_{\mathrm{eq}}, \qquad F\ge0.\]

The equality constraints enforce: \[\begin{align} &\text{local normalization or state-wise mass conservation},\\ &\text{consistency between barycenters and macroscopic gradients},\\ &\text{weak equilibrium or force balance},\\ &\text{boundary conditions on }u. \end{align}\]

4 Quantum Algorithm for Young Measures↩︎

This section examines the application of the quantum central path (QCP) algorithm to the Young-measure formulation of the homogenization problem. After reformulating the discretized Young-measure model as a finite-dimensional linear program, we review the homogeneous self-dual embedding underlying QCP and derive the corresponding query and gate complexity bounds. These complexities are then compared with those of direct fine-scale PDE solvers and classical interior-point methods. The analysis identifies parameter regimes—including small \(\varepsilon\), high-dimensional measure discretizations, and sparse LP structure—in which QCP may provide computational advantages, particularly when the Young-measure discretization is independent of the microscopic scale and the resulting LP remains sparse and well-conditioned.

4.1 Complexity of the Quantum Central Path Algorithm↩︎

We consider the finite-dimensional Young-measure linear program obtained from the discretization introduced in the previous section. Let \[z=[u;F], \qquad F\ge0,\] where \(u\) denotes the discrete macroscopic variable and \(F\) collects the nonnegative discrete Young-measure masses. The discrete problem takes the form \[\label{eq:qcp-young-measure-lp} \min_{z} \quad c^\top z \qquad \text{subject to} \qquad Az=b, \qquad z\ge0.\tag{11}\] The matrix \(A\in\mathbb{R}^{m\times n}\) collects the discrete normalization conditions, barycenter consistency relations, weak equilibrium equations, and boundary constraints.

Following the homogeneous self-dual embedding, the LP is reformulated in variables \(w\in\mathbb{R}^d, d=m+n+2\). The embedding is written as \[\label{eq:qcp-self-dual-embedding} s(w)=Mw+q, \qquad w\ge0, \qquad s(w)\ge0,\tag{12}\] where \(M^\top=-M\). For a barrier parameter \(\mu>0\), the central path is characterized by \[\label{eq:qcp-central-path} w_i s_i(w)=\mu, \qquad i=1,\dots,d.\tag{13}\]

Equivalently, define the central-path potential \[\label{eq:qcp-potential} f_\mu(w) := \frac{1}{2} \left\| w\circ s(w)-\mu\mathbf{1} \right\|_2^2.\tag{14}\] The central-path solution \(w(\mu)\) is the global minimizer of \(f_\mu\), satisfying \[\label{eq:qcp-potential-minimum} f_\mu(w(\mu))=0.\tag{15}\]

Here we give a heuristic explanation for the theorem. The QCP algorithm associates each value of the barrier parameter \(\mu\) with a Schrödinger Hamiltonian \[\label{eq:qcp-hamiltonian} H_\mu = -\frac{h}{2}\Delta + f_\mu(z),\tag{16}\] where \(h>0\) is a semiclassical parameter. Since \(f_\mu\) is minimized at the central-path point \(z(\mu)\), the ground state of \(H_\mu\) is concentrated near \(z(\mu)\).

Near \(z(\mu)\), we have the Taylor expansion \[\label{eq:qcp-taylor} f_\mu(z) = \frac{1}{2} (z-z(\mu))^\top \nabla^2 f_\mu(z(\mu)) (z-z(\mu)) + O(\|z-z(\mu)\|^3),\tag{17}\] where we have used \(f_\mu(z(\mu))=0\) and \(\nabla f_\mu(z(\mu))=0\). Thus, locally, \(H_\mu\) is approximated by a multidimensional quantum harmonic oscillator. Its ground state is approximately Gaussian and localized around \(z(\mu)\).

To recover an \(\varepsilon\)-accurate classical point in dimension \(d=m+n+2\), the wavepacket must be localized at scale \(O(\varepsilon/\sqrt d)\) in each coordinate. Indeed, if the coordinatewise error is \(O(\eta)\), then the Euclidean error is \(O(\sqrt d\,\eta)\). Hence one needs \(\displaystyle \eta = O\!\left(\frac{\varepsilon}{\sqrt d}\right)\). The total length of the central path is controlled by the assumed bound \(\|z^\star\|_1 \le R_1\). Therefore, the effective number of localization-scale steps required to transport the quantum state along the path is \(\displaystyle \frac{R_1}{\varepsilon/\sqrt d} = \sqrt d\,\frac{R_1}{\varepsilon}\). Up to logarithmic factors hidden in \(\widetilde{O}(\cdot)\), this yields the query complexity \(\displaystyle \widetilde{O}\!\left( \sqrt d\,\frac{R_1}{\varepsilon} \right)\).

Substituting \(d=m+n+2\), and absorbing constants into the asymptotic notation, gives \[\widetilde{O}\!\left( \sqrt{m+n}\,\frac{R_1}{\varepsilon} \right).\]

Finally, in the standard gate model, each query to \(f_\mu\) requires evaluating \(s(z)=Mz+q\). Since \(M\) is constructed from the original constraint matrix \(A\), this evaluation costs \(\widetilde{O}(\mathrm{nnz}(A))\) elementary gates. Multiplying this cost by the query complexity gives the overall gate complexity \[\widetilde{O}\!\left( \sqrt{m+n}\,\mathrm{nnz}(A)\frac{R_1}{\varepsilon} \right).\]

4.2 Setup for the comparison with classical solvers↩︎

We compare the QCP complexity with direct classical solvers (DCS) for the multiscale PDE, using the Young-measure LP degrees of freedom and discretization error as the bridge between the two approaches.

4.2.0.1 Degrees of freedom.

We consider both deterministic and stochastic Young-measure formulations. In the deterministic periodic setting, let \(x\in\Omega\subset\mathbb{R}^d\), \(y\in Y=[0,1]^d\), and \(\xi\in\mathbb{R}^d\). We discretize \(x\), \(y\), and \(\xi\) using \(N_x^d\), \(N_y^d\), and \(N_\xi^d\) degrees of freedom respectively. The discrete Young-measure variable is \(\mu_{i,j,\ell}\ge0\), where \(i=1,\dots,N_x^d\), \(j=1,\dots,N_y^d\), and \(\ell=1,\dots,N_\xi^d\), giving \[n_{\rm det} \sim N_x^d N_y^d N_\xi^d.\]

In the stochastic setting, we additionally let \(\omega\in\Gamma\subset\mathbb{R}^r\) be a random parameter with \(r\) random variables, discretized on a tensor grid of \(N_\omega^r\) points. The discrete Young-measure variable is \(\mu_{i,j,k,\ell}\ge0\), where \(i=1,\dots,N_x^d\), \(j=1,\dots,N_y^d\), \(k=1,\dots,N_\omega^r\), and \(\ell=1,\dots,N_\xi^d\), giving \[n_{\rm sto} \sim N_x^d N_y^d N_\omega^r N_\xi^d.\]

In both cases the LP is written abstractly as \(\min_{z\in\mathbb{R}^n} c^\top z\) subject to \(A_{\rm LP}z=b\) and \(z\ge0\), with \(m\sim n\) for local Young-measure formulations. The QCP query complexity is \(\widetilde{O}(\sqrt{n}\,R_1/\delta)\).

4.2.0.2 Discretization error.

We use first-order (\(p=1\)) discretization in \(x\), \(y\), and \(\xi\), which requires minimal regularity assumptions on the Young-measure density: \[\mathrm{err}_{\rm LP} \;\sim\; N_x^{-1} + N_y^{-1} + N_\xi^{-1} + \varepsilon_{\rm alg}.\] Balancing at target accuracy \(\delta\) gives \(N_x\sim N_y\sim N_\xi\sim\delta^{-1}\) and \(\varepsilon_{\rm alg}\sim\delta\), so \[n_{\rm det}\sim\delta^{-3d},\qquad n_{\rm sto}\sim\delta^{-3d} N_\omega^r.\] The case of exponential (spectral) convergence in \(x\) and \(y\), which requires analytic regularity of the Young-measure density in those variables (see Section 4.3), is deferred to Conjecture 1.

4.2.0.3 Desired accuracy.

We distinguish two target accuracy regimes.

• Homogenized-solution regime. The desired output is the homogenized solution \(u_0\). The natural accuracy threshold is \[\delta_{\rm hom} = \varepsilon^\alpha, \quad 0<\alpha<1,\] since the homogenization error satisfies \(\|u_\varepsilon-u_0\|\sim\varepsilon^\alpha\). This threshold is coarser than resolving the full oscillatory field.

• Fine-solution regime. The desired output is the full fine-scale solution \(u_\varepsilon\), with threshold \(\delta_{\rm fine}=\varepsilon\). Under \(p=1\) algebraic discretization this regime does not yield quantum advantage (see Conjecture 1 for the spectral case).

4.3 Possible quantum advantage regimes↩︎

We now describe three regimes in which QCP may outperform DCS, using the setup from Section 4.2. Since \(m\sim n\) for the sparse local Young-measure LPs considered here, the QCP query complexity reduces to \(\widetilde{O}\!\bigl(\sqrt{n}\,R_1/\delta\bigr)\).

4.3.0.1 1. Nonlinear homogenization with homogenized accuracy (\(\delta=\varepsilon^\alpha\), \(0<\alpha<1\), \(p=1\)).

We seek the homogenized solution \(u_0\) to accuracy \(\delta=\varepsilon^\alpha\) with \(0<\alpha<1\). This is the natural target in nonlinear homogenization, where the homogenization error satisfies \(\|u_\varepsilon-u_0\|\sim\varepsilon^\alpha\) with \(\alpha<1\) in general norms (for instance \(\alpha=1/2\) in \(H^1\) for smooth linear periodic problems, and possibly smaller for nonlinear settings). A direct fine-scale solver must resolve all \(\varepsilon\)-scale oscillations regardless of the desired accuracy for \(u_0\), incurring cost \[\mathrm{DCS}_{\rm det}\sim\widetilde{O}(\varepsilon^{-d}).\]

For the Young-measure LP with first-order (\(p=1\)) algebraic convergence in \(x\), \(y\), and \(\xi\), balancing all discretization errors at \(\delta=\varepsilon^\alpha\) gives \(N_x\sim N_y\sim N_\xi\sim\varepsilon^{-\alpha}\), so \(n_{\rm det}\sim\varepsilon^{-3\alpha d}\), and \[\mathrm{QCP}_{\rm det} \sim \widetilde{O}\!\left( R_1\,\varepsilon^{-\alpha(2+3d/2)} \right) = \widetilde{O}\!\left( R_1\,\varepsilon^{-\alpha(3d+4)/2} \right).\] QCP beats DCS when \(\alpha(3d+4)/2 < d\), i.e., \[\boxed{\alpha < \frac{2d}{3d+4}.}\] For \(d=2\) this requires \(\alpha<4/10=2/5\); for \(d=3\), \(\alpha<6/13\); for \(d>>1\), \(\alpha\lesssim 2/3\). Since \(\alpha<1\) is the natural regime for nonlinear homogenization, this condition can be satisfied for sufficiently small \(\alpha\) in all dimensions \(d\ge1\).

4.3.0.2 2. Stochastic homogenization with fine-scale accuracy (\(\delta=\varepsilon\), \(p=1\)).

The second advantage regime arises in stochastic homogenization with \(r\) random variables \(\omega\in\mathbb{R}^r\), where we seek the fine-scale solution at accuracy \(\delta=\varepsilon\). We treat the total number of random discretization points \(N_\omega^r\) as a free parameter, independent of the spatial accuracy \(\varepsilon\), since the convergence rate in \(\omega\) depends on the smoothness of the solution in the random parameter and need not match the spatial rate. A direct solver must handle each of the \(N_\omega^r\) realizations at fine-scale cost \(\varepsilon^{-d}\): \[\mathrm{DCS}_{\rm sto}\sim\widetilde{O}(N_\omega^r\,\varepsilon^{-d}).\] The stochastic Young-measure LP encodes all random realizations simultaneously in a single LP. With \(p=1\) algebraic convergence in \(x\), \(y\), \(\xi\) and \(\delta=\varepsilon\), \(N_x\sim N_y\sim N_\xi\sim\varepsilon^{-1}\), so \(n_{\rm sto}\sim\varepsilon^{-3d}\,N_\omega^r\), and \[\mathrm{QCP}_{\rm sto} \sim \widetilde{O}\!\left( R_1\,\varepsilon^{-(3d/2+1)}\,N_\omega^{r/2} \right).\] QCP beats DCS when \(\varepsilon^{-(3d/2+1)}\,N_\omega^{r/2}\lesssim N_\omega^r\,\varepsilon^{-d}\), i.e., \[\boxed{N_\omega^r \gtrsim \varepsilon^{-(d+2)}.}\] This condition depends only on the total stochastic dimension \(N_\omega^r\) and not on \(r\) or \(N_\omega\) individually: whenever the random environment is discretized finely enough that \(N_\omega^r\gg\varepsilon^{-(d+2)}\), quantum advantage holds. The quantum speedup factor is \(N_\omega^{r/2}\) — a square root in the stochastic dimension — which grows without bound as \(N_\omega^r\) increases, independently of \(\varepsilon\). When the advantage condition is met, the DCS cost is at least \(\varepsilon^{-(2d+2)}\), while the QCP cost is \(\widetilde{O}(R_1\,\varepsilon^{-(3d/2+1)}\,N_\omega^{r/2})\), giving a genuine polynomial separation.

4.3.0.3 3. Conjecture: quantum advantage under spectral discretization.

Above discussions use first-order algebraic convergence in \(x\), \(y\), and \(\xi\). As noted in Remark 11, one may improve convergence in \(x\) and \(y\) while keeping first-order in \(\xi\), since regularity of \(\nu_{x,y}\) in \(\xi\) is not inherited from smoothness of \(a\). In the extreme case where \(\nu_{x,y}(\xi)\) is analytic in \(x\) and \(y\), spectral discretization gives \(N_x\sim N_y\sim\log(1/\varepsilon)\) while \(N_\xi\sim\varepsilon^{-1}\), reducing the LP size from \(n\sim\varepsilon^{-3d}\) (all algebraic) to \(n_{\rm det}\sim\varepsilon^{-d}|\log\varepsilon|^{2d}\) at fine-scale accuracy \(\delta=\varepsilon\).

5 Numerical Experiments↩︎

We validate the proposed Young-measure LP formulation on several benchmark homogenization problems, including 1D deterministic periodic media, 1D stochastic media, and 2D anisotropic periodic media. The primary purpose of these experiments is to verify that the LP formulation correctly reproduces the expected homogenized solutions in linear, nonlinear and stochastic settings.

The numerical error is measured using the relative \(L^\infty\)-error \[\label{eq:rel95err95def} \frac{ \|u_{\mathrm{num}}-u_{\mathrm{exact}}\|_\infty }{ \|u_{\mathrm{exact}}\|_\infty }.\tag{18}\]

Each benchmark uses a problem-dependent discretization consisting of:

• macroscopic spatial grids,

• microscopic periodic or random-state discretizations,

• gradient-state discretizations for the Young measures,

• and the associated constitutive law.

The corresponding parameters are summarized in the tables preceding each experiment.

The experiments compare the LP-computed homogenized solution with analytically known homogenized solutions. In all cases considered below, the numerical results exhibit good agreement with the exact homogenized profiles and small relative errors, supporting the correctness and consistency of the Young-measure LP formulation.

5.1 1D Deterministic Periodic Media↩︎

5.1.1 A linear case↩︎

The oscillatory problem is \(\displaystyle -\frac{d}{dx}\left(k\!\left(\frac{x}{\varepsilon}\right)u_\varepsilon'(x)\right)=1, u_\varepsilon(0)=u_\varepsilon(1)=0\), with periodic coefficient \(k(y)=2+\sin(2\pi y)\).

The effective coefficient is \(\displaystyle k_{\mathrm{hom}}=\left(\int_0^1\frac{1}{k(y)}\,dy\right)^{-1}=\sqrt{3}\), and the homogenized problem \(\displaystyle -\frac{d}{dx}\left(k_{\mathrm{hom}}u_0'(x)\right)=1, u_0(0)=u_0(1)=0, u_{\mathrm{exact}}(x)=\frac{x(1-x)}{2k_{\mathrm{hom}}}\). Numerically, the LP solution agrees closely with the parabola, with the pointwise relative error remaining uniformly at the \(10^{-5}\) level across the domain.

5.1.1.1 Discussion.

In the linear case, the Young measure is theoretically a Dirac mass \(\delta_{\xi^*(x,y)}\) at each \((x,y)\), where \(\xi^*(x,y) = u_0'(x) + w_y'(y,u_0'(x))\) is determined by the cell corrector \(w\). The marginal Young-measure distribution shown in Figure 2 therefore reflects the pushforward of Lebesgue measure on \(Y\) through the analytic map \(y\mapsto\xi^*(x,y)\): it is smooth and unimodal, concentrated near the macroscopic gradient \(u_0'(x)\) at each \(x\). Since \(k(y)=2+\sin(2\pi y)\) is analytic, both \(\xi^*(x,y)\) and the effective coefficient \(k_{\rm hom}=\sqrt{3}\) are determined by analytic corrector computations. The excellent agreement at the \(10^{-5}\) level confirms that the LP correctly recovers this analytic structure through the discrete gradient-state representation.

5.1.2 Nonlinear case↩︎

The deterministic nonlinear benchmark is \[-\frac{d}{dx}\left(k\!\left(\frac{x}{\varepsilon}\right)u_\varepsilon'(x)^3\right)=1, \qquad u_\varepsilon(0)=u_\varepsilon(1)=0,\] with the same \(k(y)=2+\sin(2\pi y)\). The effective coefficient is \[k_{\mathrm{hom}}=\left(\int_0^1k(y)^{-1/3}\,dy\right)^{-3},\] and the exact homogenized profile is \[u_{\mathrm{exact}}(x)=\frac{3}{4}k_{\mathrm{hom}}^{-1/3}\left(2^{-4/3}-\left|x-\tfrac12\right|^{4/3}\right).\] The computed solution captures the non-quadratic shape accurately.

5.1.2.1 Discussion.

The nonlinear cubic-flux operator \(a(y,\xi)=k(y)\xi^3\) satisfies the strict monotonicity condition \(\partial_\xi a>0\) for \(k>0\), so the cell corrector map \(y\mapsto\xi^*(x,y)\) is a diffeomorphism on \(Y\) for each fixed macroscopic gradient. Accordingly, the marginal Young measure (Figure 4) is again the pushforward of Lebesgue measure on \(Y\) through this map and appears as a smooth, spread distribution in \(\xi\). The non-quadratic profile \(u_{\rm exact}\propto|x-1/2|^{4/3}\) has a Hölder singularity at \(x=1/2\), where \(u_{\rm exact}'=0\); this causes a mild degradation in the spatial convergence rate near the center of the domain (visible in the error panel), consistent with the algebraic \(O(N_x^{-4/3})\) rate expected for \(W^{1,\infty}\)-approximation of a \(C^{4/3}\) function. The LP nonetheless captures the overall shape with errors below \(10^{-3}\), demonstrating robustness to mild solution singularities.

5.2 2D Deterministic Periodic Media↩︎

For both 2D tests, \[\Omega=(0,1)^2, \qquad u=0\text{ on }\partial\Omega,\] with anisotropic periodic coefficients \[k_1(s)=2+\sin(2\pi s),\qquad k_2(s)=3+\cos(2\pi s),\] and manufactured exact field \[u_{\mathrm{exact}}(x_1,x_2)=x_1(1-x_1)x_2(1-x_2).\]

5.2.1 Linear anisotropic case↩︎

The oscillatory equation is \[-\nabla\cdot\left(\begin{bmatrix}k_1(x_1/\varepsilon)&0\\0&k_2(x_2/\varepsilon)\end{bmatrix}\nabla u_\varepsilon\right)=f_{\mathrm{lin}}(x), \qquad u_\varepsilon|_{\partial\Omega}=0.\] The right-hand side is chosen as \[f_{\mathrm{lin}}(x_1,x_2)=2k_{1,\mathrm{eff}}x_2(1-x_2)+2k_{2,\mathrm{eff}}x_1(1-x_1),\] where \(k_{1,\mathrm{eff}}\) and \(k_{2,\mathrm{eff}}\) are the homogenized coefficients.

5.2.1.1 Discussion.

The anisotropic coefficient structure decouples the two spatial directions: the effective coefficients \(k_{1,\rm eff}=\sqrt{3}\) and \(k_{2,\rm eff}=1/(\int_0^1(3+\cos 2\pi s)^{-1}ds)\) are determined independently by 1D cell problems in \(y_1\) and \(y_2\) respectively. The marginal Young-measure distributions (Figure 6) reflect this decoupling: the support in the \((\xi_1,\xi_2)\)-plane is concentrated along a curve determined by the product of two independent corrector maps. The maximum relative error of \(1.561\times10^{-2}\) is consistent with the coarser 2D discretization (\(20\times20\) macro grid, \(5\times5\) micro grid, \(25\times25\) gradient grid), where the LP dimension \(n\sim(20)^2(5)^2(25)^2\approx 1.56\times10^7\) is already large for classical solvers. This illustrates the LP dimension scaling \(n\sim N_x^d N_y^d N_\xi^d\) in 2D and motivates the use of quantum LP solvers for larger problems.

5.2.2 Nonlinear anisotropic case↩︎

The oscillatory equation is \[-\partial_{x_1}\!\left(k_1(x_1/\varepsilon)(\partial_{x_1}u_\varepsilon)^3\right) -\partial_{x_2}\!\left(k_2(x_2/\varepsilon)(\partial_{x_2}u_\varepsilon)^3\right) =f_{\mathrm{nonlin}}(x), \qquad u_\varepsilon|_{\partial\Omega}=0.\]

The right-hand side is chosen as \[f_{\mathrm{nonlin}}(x_1,x_2)=6K_{1,\mathrm{eff}}x_2^3(1-x_2)^3(1-2x_1)^2+6K_{2,\mathrm{eff}}x_1^3(1-x_1)^3(1-2x_2)^2,\] where \(K_{1,\mathrm{eff}}\) and \(K_{2,\mathrm{eff}}\) are the homogenized coefficients.

5.2.2.1 Discussion.

The 2D nonlinear case combines the decoupled anisotropic structure with a cubic-flux nonlinearity in each direction. The effective coefficients \(K_{i,\rm eff}=(k_{i,\rm eff}^{-1/3})^{-3}\) are determined by the nonlinear 1D cell problems and appear as the independent cubic homogenization constants for each direction. The marginal Young-measure distributions (Figure 8) are more spread in the \(\xi\)-plane compared to the linear case, reflecting the stronger nonlinear distortion of the corrector map \(\xi^*(x,y)\) under the cubic flux. The maximum relative error of \(8.599\times10^{-3}\) is achieved on a \(15\times15\) macro grid with \(35\times35\) gradient states, and is smaller than the 2D linear error despite the coarser macro grid, because the nonlinear effective profile \(u_{\rm exact}=x_1(1-x_1)x_2(1-x_2)\) is smooth and well-resolved by the LP. The result confirms that the Young-measure LP handles nonlinear constitutive laws in multiple dimensions without structural modification.

5.3 1D Stochastic Media↩︎

5.3.1 Random linear case↩︎

We solve \[-\frac{d}{dx}\left(c\!\left(\frac{x}{\varepsilon},\omega\right)u_\varepsilon'(x,\omega)\right)=1, \qquad u_\varepsilon(0,\omega)=u_\varepsilon(1,\omega)=0,\] with \(c\in\{5,1\}\) and probabilities \([0.4,0.6]\). The effective coefficient is \[c_{\mathrm{hom}}=\left(\sum_i p_i/c_i\right)^{-1},\] and \(u_{\mathrm{exact}}(x)=x(1-x)/(2c_{\mathrm{hom}})\). Script output gives \[\max|u_{\mathrm{num}}-u_{\mathrm{exact}}|=7.36\times10^{-4}, \qquad \frac{\|u_{\mathrm{num}}-u_{\mathrm{exact}}\|_\infty}{\|u_{\mathrm{exact}}\|_\infty}=8.66\times10^{-3}.\]

5.3.1.1 Discussion.

With only two random states \(\omega\in\{c_1,c_2\}\), the stochastic Young measure \(\nu_{x,i}(\xi)\) for each realization \(i\) is a Dirac mass \(\delta_{\xi^*(x,i)}\) at the corrector gradient for that realization. The global marginal \(\hat{\mu}_x(\xi)=\sum_i p_i\delta_{\xi^*(x,i)}\) is therefore a weighted combination of two Dirac masses, visible as two discrete peaks in Figure 10. The LP correctly identifies both peaks and their weights, recovering the harmonic average effective coefficient \(c_{\rm hom}\) to within relative error \(8.66\times10^{-3}\). This stochastic case with \(N_\omega=2\) represents the simplest instance of the large-\(r\) advantage regime: the LP encodes both realizations simultaneously in a single solve, at half the DCS cost of solving two independent fine-scale problems.

5.3.2 Random nonlinear case↩︎

We solve \[-\frac{d}{dx}\left(c\!\left(\frac{x}{\varepsilon},\omega\right)|u_\varepsilon'(x,\omega)|u_\varepsilon'(x,\omega)\right)=1, \qquad u_\varepsilon(0,\omega)=u_\varepsilon(1,\omega)=0,\] with three random states \(c\in\{1,3,6\}\) and probabilities \([0.5,0.3,0.2]\). The effective coefficient is \[c_{\mathrm{hom}}=\left(\sum_i p_i c_i^{-1/2}\right)^{-2},\] and the exact profile is \(u_{\mathrm{exact}}(x)=\frac{2}{3}c_{\mathrm{hom}}^{-1/2}\left(2^{-3/2}-\left|x-\tfrac12\right|^{3/2}\right)\).

Script output gives \[\max|u_{\mathrm{num}}-u_{\mathrm{exact}}|=1.27\times10^{-4}, \qquad \frac{\|u_{\mathrm{num}}-u_{\mathrm{exact}}\|_\infty}{\|u_{\mathrm{exact}}\|_\infty}=7.11\times10^{-4}.\]

5.3.2.1 Discussion.

With three random states and a quadratic flux \(a(\omega)|\xi|\xi\), the LP must simultaneously represent three distinct corrector gradients weighted by their probabilities \([0.5,0.3,0.2]\). The marginal Young measure (Figure 12) shows three discrete peaks at \(\xi\)-values determined by \(\xi^*(x,i)\propto c_i^{-1/2}\,{\rm sgn}(u_0'(x))\), each with height proportional to \(p_i\). The effective coefficient \(c_{\rm hom}=(\sum_i p_i c_i^{-1/2})^{-2}\) blends the three realizations nonlinearly via the cubic homogenization formula. The very small relative error \(7.11\times10^{-4}\) demonstrates that the LP correctly resolves the nonlinear stochastic averaging, recovering a relative error nearly an order of magnitude smaller than the linear stochastic case. This improvement is consistent with the smooth \(|x-1/2|^{3/2}\) profile being slightly less singular than the linear parabola under the cubic flux energy, and with the finer effective resolution provided by three gradient peaks compared to two.

5.3.3 Linear non-variational anisotropic case↩︎

The oscillatory equation is given by \[-\nabla\cdot\left(A(x/\varepsilon)\nabla u_\varepsilon\right)=f_{\mathrm{nonvar}}(x), \qquad u_\varepsilon|_{\partial\Omega}=0,\] where the micro-scale coefficient matrix \(A(y)\) is explicitly asymmetric. Due to the presence of the skew-symmetric terms (\(a_{12}\) and \(a_{21}\)), the system lacks an underlying scalar energy potential \(W(y,\xi)\) such that \(\nabla_\xi W = A(y)\xi\). Consequently, the problem is strictly non-variational.

To uniquely determine the homogenized solution within the Young-measure LP framework, the objective function is set to zero (a pure feasibility problem), and the micro-scale divergence-free condition must be explicitly enforced as an equality constraint: \[-\nabla_y \cdot \left( A(y)\xi \right) = 0,\] along with the microscopic curl-free constraint \[\nabla_y \times \mathbb{E}[\xi] = 0.\] To avoid unnatural non-vanishing measure on the phase-space boundary, we impose an artificial energy regularization objective (\(c(\xi) = |\xi|^2\)) to enforce the collapse of the Young measure into a deterministic Dirac mass. Given the severe degeneracy caused by the measure collapse on the refined phase-space grid, the global LP system is efficiently solved using the interior-point algorithm of the linprog function in MATLAB. For the method of manufactured solutions, we design the same exact macroscopic solution as \(u_0(x) = x_1(1-x_1)x_2(1-x_2)\). The corresponding right-hand side is chosen as

\[f_{\mathrm{nonvar}}(x_1,x_2)=2\sqrt{3}x_2(1-x_2)+4x_1(1-x_1),\]

based on the analytically derived homogenized matrix \(A_{\mathrm{eff}} = \begin{bmatrix} \sqrt{3} & 1 \\ -1 & 2 \end{bmatrix}\).

5.3.3.1 Discussion.

The non-variational case is qualitatively different from the preceding examples: there is no underlying energy \(W(y,\xi)\) whose gradient equals \(A(y)\xi\), so the LP cannot be posed as an energy minimization. The feasibility-only formulation with zero objective identifies the unique homogenized solution via the combined microscopic divergence-free and macroscopic curl-free constraints. The energy regularization objective \(c(\xi)=|\xi|^2\) is introduced artificially to collapse the Young measure to a Dirac mass, preventing the LP from exploiting degenerate feasible directions in the gradient-state space. The marginal distributions (Figure 14) accordingly concentrate near a single point in \((\xi_1,\xi_2)\) at each macro-cell. The maximum relative error of \(1.2\times10^{-2}\) confirms that even in the non-variational setting, the LP feasibility approach correctly reproduces the analytically derived homogenized matrix \(A_{\rm eff}=\begin{bmatrix}\sqrt{3}&1\\-1&2\end{bmatrix}\), demonstrating the scope of the Young-measure LP framework beyond energy minimization.

5.4 Summary on Numerical Results↩︎

The numerical experiments show that the proposed Young-measure LP formulation accurately reproduces known homogenized solutions across deterministic periodic, stochastic, linear, nonlinear, and higher-dimensional elliptic homogenization problems. In the 1D deterministic and stochastic cases, the LP solutions closely match analytical homogenized profiles for both linear and cubic-flux models with very small relative errors. The 2D anisotropic periodic benchmarks further demonstrate agreement with manufactured exact solutions in both linear and nonlinear settings, yielding maximum relative errors of the order of \(10^{-2}\), which is quite acceptable given the relatively rough grids limited by the huge number of variables introduced by higher dimensions. Moreover, the non-variational case suggests that our framework also works for problems that do not admit an energy minimization formulation. Overall, these results confirm the correctness, consistency, and flexibility of the Young-measure LP framework for homogenization problems involving periodic media, random media, and nonlinearity.

6 Conclusion and Outlook↩︎

This work develops a Young-measure linear programming framework for nonlinear deterministic and stochastic homogenization and investigates its compatibility with quantum optimization algorithms. The Young-measure formulation lifts the nonlinear homogenization problem to a linear one in higher dimensions by treating the microscale \(y\), the gradient \(\xi\), and the random environment \(\omega\) as independent variables, thereby characterizing effective macroscopic quantities through a structured LP without directly resolving fine-scale oscillations.

The complexity analysis compares the quantum central path (QCP) algorithm against direct classical solvers (DCS), whose cost \(\widetilde{O}(\varepsilon^{-d})\) is fixed by the need to resolve all \(\varepsilon\)-scale oscillations regardless of the desired output accuracy. Using first-order (\(p=1\)) algebraic discretization, two regimes of quantum advantage are identified:

1. Deterministic nonlinear homogenization. Recovering the homogenized solution to accuracy \(\delta=\varepsilon^\alpha\) yields QCP cost \(\widetilde{O}(R_1\varepsilon^{-\alpha(3d+4)/2})\), improving over DCS cost \(\widetilde{O}(\varepsilon^{-d})\) whenever \(\alpha<2d/(3d+4)\). For \(d=2\) this requires \(\alpha<2/5\); for \(d=3\), \(\alpha<6/13\).

2. Stochastic homogenization at fine-scale accuracy. With \(N_\omega^r\) random discretization points, QCP achieves cost \(\widetilde{O}(R_1\,N_\omega^{r/2}\,\varepsilon^{-(3d/2+1)})\), giving a square-root reduction in stochastic sampling cost that grows with \(N_\omega^r\), whenever \(N_\omega^r\gtrsim\varepsilon^{-(d+2)}\).

We further conjecture that if the Young-measure density is analytic in \(x\) and \(y\), spectral discretization in those variables reduces the LP size to \(n_{\rm det}\sim\varepsilon^{-d}|\log\varepsilon|^{2d}\) at fine-scale accuracy \(\delta=\varepsilon\), yielding QCP cost \(\widetilde{O}(R_1\varepsilon^{-(1+d/2)}|\log\varepsilon|^d)\) and quantum advantage for all \(d\ge2\). In the stochastic case, the advantage condition under spectral discretization relaxes to \(N_\omega^r\gtrsim\varepsilon^{-(d-2)}\), trivially satisfied for \(d\ge2\).

Numerical experiments on 1D and 2D deterministic and stochastic benchmarks confirm the correctness and consistency of the Young-measure LP formulation.

Future work will focus on: rigorous convergence analysis for the Young-measure LP in appropriate function spaces and norms; verification of the spectral conjecture, including analysis of the central-path scale \(R_1\) under spectral discretization and of the analytic regularity of the Young-measure density; extensions to broader classes of nonlinear and stochastic PDEs; and comprehensive quantum complexity estimates including state preparation, oracle construction, and measurement costs beyond the query-complexity level. Overall, the proposed framework establishes a bridge among nonlinear homogenization, stochastic multiscale PDEs, linear programming, and quantum scientific computing, and points a promising direction for efficient quantum computation of nonlineaer multiscale and stochastic PDEs.

Acknowledgments↩︎

SJ and ZL were supported by NSFC grant No. 12341104, the Shanghai Pilot Program for Basic Research, the Shanghai Jiao Tong University 2030 Initiative, the Science and Technology Innovation Key R&D Program of Chongqing grant No. CSTB2024TIAD-STX0035, and the Fundamental Research Funds for the Central Universities.

References↩︎