Pith. sign in

REVIEW 4 major objections 4 minor 57 references

Extended Walk-on-Spheres Algorithm for Linear and Nonlinear Elliptic Problems of Divergence-type

T0 review · 4 major / 4 minor · reviewed 2026-08-10 · deepseek-v4-flash

Pith's one-line read The paper proposes a reusable Monte Carlo toolbox that turns divergence-type elliptic Dirichlet problems into small algebraic systems.

desk verdict The LOT precomputation is a genuinely new idea for WoS-type solvers, and the experiments are encouraging, but the missing error analysis for the algebraic solve step leaves the 'general solver' claim heuristic. read the letter →

arxiv 2608.07017 v1 pith:EQEQD22O submitted 2026-08-07 math.NA cs.NAmath.APmath.PR

classification math.NAcs.NAmath.APmath.PR MSC 65N7535J2565C0565N80
keywords Walk-on-SpheresMonteCarloellipticPDEDirichletproblemharmonicextensionGreenfunctionVoronoiprobabilisticnumericalmethod
verification ladder T0 review T1 audit T2 compute T3 formal

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The reading

The paper claims that the classical Walk-on-Spheres method can be extended from the Laplace equation to general divergence-type elliptic problems, linear and nonlinear, without deriving a new stochastic representation for each operator. Its strategy is to precompute, once per domain, random matrices that approximate the inverse Dirichlet Laplacian, the harmonic extension operator, and their gradients, forming a Laplacian Operator Toolbox. Any Dirichlet boundary value problem is then approximated by writing the solution as the toolbox applied to an unknown source term plus the harmonic extension of the boundary data, and solving a finite-dimensional algebraic system for that source term. If this reduction is correct, the same precomputed toolbox solves many different PDEs on the same domain, and for linear operators it also supplies estimators of the Green integral operator and the elliptic measure operator.

What carries the argument

The central object is the Laplacian Operator Toolbox: random matrices that approximate $U$, $\nabla U$, $H$, and $\nabla H$ after a Voronoi discretization of the domain and boundary, defined in (2.27). The matrices are Monte Carlo estimates of occupation times and exit locations of Walk-on-Spheres chains; entries of $A_{M,N}$ and $B^{(j)}_{M,N}$ accumulate weighted visits to spatial cells, while entries of $H_{M,N}$ and $K^{(j)}_{M,N}$ record where chains first approach the boundary. The argument then uses the representation $u\approx Uq+Hg$ to reduce the linear PDE to the algebraic system (3.5), built from the toolbox matrices and the coefficients, and reduces Euler-Lagrange problems to minimizing or finding saddle points of the discrete functional (3.12).

What would settle it

Take a domain with a known solution, such as the Poisson problem on the unit disk, precompute the LOT matrices with large Monte Carlo sample counts, and solve (3.5) for increasing numbers of interior and boundary points. If the approximate solution does not converge to the exact one as the Voronoi cell diameters tend to zero, or if ever stronger truncation of the singular values is needed as the discretization refines, then the point-evaluation quadrature or the invertibility step behind (3.5) is not sound.

Watch

Extended reading notes

Core claim

The central claim is that a precomputed Laplacian Operator Toolbox (LOT), consisting of the random matrices $A_{M,N}$, $B^{(j)}_{M,N}$, $H_{M,N}$, and $K^{(j)}_{M,N}$ that approximate $U=(-\frac{1}{2}\Delta_0)^{-1}$, its gradient, the harmonic extension $H$, and its gradient, as collected in (2.27), is consistent (Proposition 2.13) and can be used to solve arbitrary Dirichlet boundary value problems of divergence type. For linear operators, the solution is represented as $u\approx Uq+Hg$ with an unknown source $q$; testing the PDE against functions $U\varphi$ and using the toolbox matrices turns the weak form into the linear algebraic system $\Lambda\bar q = -A^\lambda_{M,N}\tilde f - \Sigma\tilde g$ (3.5), whose solution yields the approximate solution at the discretization points. The same formula also gives Monte Carlo estimators of the Green integral operator $G_L$ and the elliptic measure operator $H_L$. For nonlinear problems arising as Euler-Lagrange equations, the toolbox converts the energy functional into a finite-dimensional functional $J_{n,m,p,M,N}$ on $\mathbb{R}^m$, which is minimized or whose saddle points are searched. The numerical experiments cover non-symmetric and anisotropic linear equations, $p$-Laplace, minimal surface, and semilinear equations with multiple solutions.

Load-bearing premise

The reduction of the PDE to the algebraic system assumes that evaluating functions at representative points inside Voronoi cells gives accurate integrals, and that the resulting matrix $\Lambda$ is invertible or can be safely regularized; the paper proves consistency only for the Monte Carlo matrices, not for this reduction, and its examples require truncated SVD because $\Lambda$ is ill-conditioned.

Editorial extensions

If this is right

  • For any single domain, one precomputed toolbox can be reused to solve many different linear and nonlinear Dirichlet problems without rebuilding a mesh or deriving new exit laws.
  • For a general linear second-order elliptic operator, the same matrices provide estimators of the Green integral operator and the elliptic measure operator through formula (3.6).
  • Nonlinear Euler-Lagrange problems, including $p$-Laplace and minimal surface equations, are reduced to finite-dimensional optimization problems whose solutions are then mapped back to approximate PDE solutions.
  • The method inherits the mesh-free and parallel-friendly character of Walk-on-Spheres, so the toolbox construction scales favorably with dimension compared to mesh-based solvers.

Reading between the lines

Editorial extensions of the paper, not claims the author makes directly.

  • The paper leaves open a convergence proof for the algebraic reduction itself: Proposition 2.13 proves consistency of the LOT matrices, but not that solutions of (3.5) converge to the PDE solution as the Voronoi cells refine, a gap the ill-conditioning of $\Lambda$ in the examples makes concrete.
  • If the reduction step can be rigorously stabilized, the same toolbox could plausibly be applied in higher dimensions, where mesh-based solvers degrade, but whether the conditioning of $\Lambda$ worsens with dimension is untested.
  • The Green-operator estimator in (3.6) hints at a possible one-step propagator for parabolic or time-dependent problems, though the paper does not pursue that direction.
Share X Bluesky LinkedIn Reddit HN

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, and a circularity audit.

Referee Report

4 major / 4 minor

Summary. The paper proposes a Walk-on-Spheres-based Laplacian Operator Toolbox (LOT): random matrices A_{M,N}, B_{M,N}, H_{M,N}, K_{M,N} that approximate the inverse Dirichlet Laplacian U, its gradient, the harmonic extension H, and its gradient. Section 2 develops probabilistic representations, Monte Carlo estimators, and consistency results (Propositions 2.7 and 2.13). Section 3 uses the LOT to reduce linear divergence-type Dirichlet problems to the finite-dimensional system (3.5), and Euler-Lagrange/energy minimization problems to the optimization of J_{n,m,p,M,N} in (3.14). Numerical experiments on Laplace/Poisson problems, non-symmetric/anisotropic linear equations, p-Laplace, minimal surface, and semilinear equations are reported. The central claim is that the precomputed LOT transforms arbitrary divergence-type Dirichlet problems into an algebraic system whose solution yields an approximate PDE solution, and for linear operators gives estimators of the Green and elliptic-measure operators.

Significance. The probabilistic representations and LOT consistency proofs are nontrivial and appear correct; the paper supplies full proofs in Appendix A.2, and the numerical benchmark suite is broader than in most WoS papers. If the algebraic reduction were backed by convergence and stability analysis, the method would be a significant contribution: it offers a mesh-free, parallel building block for variable-coefficient and nonlinear elliptic problems without custom stochastic representations. The downloadable matrices and pseudocode strengthen reproducibility. However, the current manuscript does not establish the central approximation claim; the missing step is exactly the passage from consistent LOT matrices to controlled errors in the solution of (3.5) or of the minimization problem (3.14).

major comments (4)
  1. [§3.1, Eq. (3.5)] The system (3.5) is the load-bearing step of the paper, yet no theorem states that its solution \bar q produces \bar u close to the weak solution u. Proposition 2.13 proves only entrywise limits for A_{M,N}, B^{(j)}_{M,N}, H_{M,N}, K^{(j)}_{M,N} in the triple limit m→∞, M→∞, N→∞ for fixed data and fixed output points; it does not provide norm convergence of these matrices, and it says nothing about the spectrum of Λ, which is a quadratic function of the random matrices. Entrywise consistency of the building blocks does not imply that the pseudoinverse solution of (3.5) is close to any q for which Uq+Hg approximates u. A convergence theorem for \bar q and for \bar u, under explicit conditions on Λ, is required for the advertised claim.
  2. [§3.1, Remark 3.3 and Eq. (3.2)] The derivation of (3.2) replaces integrals over W_i by λ(W_i) times integrand values at the centers x_i, including f(x_i), b(x_i), c(x_i). The PDE assumptions in the same section only require bounded measurable coefficients and f∈L^2(D); pointwise evaluation is not defined for such data. No quadrature error estimate is given, and the regularity needed to justify the quadrature is not stated as a hypothesis of any theorem. Since (3.2) determines Λ and Σ in (3.5), this unquantified approximation is load-bearing; moreover, the limit n→∞ of the output partition W_n is not analyzed anywhere in the consistency results.
  3. [§A.1, Table 7 and Eq. (3.6)] The reported spectra σ_max(Λ)≈10^-2, σ_min(Λ)≈10^-7 (Example 3) and σ_max≈10^-3, σ_min≈10^-7 (Example 4) show that Λ is severely ill-conditioned. The paper uses TSVD with relative thresholds rcond∈{10^-3,10^-2} chosen per example, but provides no error bound for the effect of truncation or of Monte Carlo noise on \bar q, nor on the final estimator (3.6) of G_L f and H_L g. Without such a bound, the interpretation of \bar u as an approximate PDE solution, and of the matrices in (3.6) as estimators of the Green and elliptic-measure operators, remains heuristic.
  4. [§3.2, Eq. (3.14)] For the nonlinear problems the paper solves the minimization of J_{n,m,p,M,N} or its critical-point equations with Adam or p-HiSD; Table 7 lists learning rates, initial guesses, and saddle indices. The reported quantity ∥∇J(\bar f_opt)∥_{2,V} is the gradient norm of the discretized algebraic objective, not the PDE residual, and no theorem relates stationary points of J_{n,m,p,M,N} to solutions of the Euler-Lagrange equation (3.9) as n,m,M,N grow. This is an additional unquantified layer between the LOT and the claimed approximate solution.
minor comments (4)
  1. [§2.3, Proposition 2.13] In the last displayed limit, [K^{(j)}_{M,N} \bar f]_i should be [K^{(j)}_{M,N} \tilde g]_i, since K acts on boundary values; as written it conflicts with the definition (2.26).
  2. [§2.4, Algorithms 4 and 5] K^{(j)}_{M,N} is initialized as an n×m zero matrix in both pseudocode listings, but from (2.26) and the subsequent loops it should be n×p; both occurrences should be corrected.
  3. [§2.1, heading] The heading 'T aming the singular terms' contains a spacing typo and should read 'Taming the singular terms'.
  4. [§2.3, Remark 2.15] Remark 2.15 refers to 'A_{m,n}' while the toolbox definition (2.27) and Algorithms use A_{M,N}; please harmonize the notation.

Circularity Check

0 steps flagged · score 0.0 of 10

No significant circularity: LOT estimators are derived from first principles and the PDE reduction is a Galerkin-type scheme.

full rationale

The paper's central derivation chain is self-contained. Proposition 2.2 and Proposition 2.5 state probabilistic representations for U, ∇U, H, and ∇H, and proofs are provided in Appendix A.2, so the LOT matrices (2.23)-(2.26) are not defined in terms of the target PDEs. Proposition 2.13 proves entrywise consistency via the law of large numbers and partition refinement, using Propositions 2.7 and 2.11, again with proofs supplied in the paper. The reduction of a divergence-type BVP to the algebraic system (3.5) is a Galerkin-type method: the ansatz u ≈ Uq + Hg is inserted into the weak form and tested against Uφ; the resulting system is solved for q, and the LOT matrices are used only as numerical approximations of the building-block operators. Nothing in this chain fits a parameter to the target solution or defines the target in terms of the output. The Green operator and elliptic measure estimators in Section 3.1 are the numerical solution operator expressed in terms of the same precomputed matrices; this is an approximation claim, not a definitional identity, since G_L and μ_L are defined independently as the true solution operators. Self-citations [10], [16], and [17] appear for auxiliary facts such as distance-approximation guarantees, boundary Hölder estimates, and harmonic measure absolute continuity, but the corresponding statements are either proved in the paper, supported by external references, or not load-bearing for the main consistency result. The principal weakness is the lack of a convergence or stability analysis for the algebraic reduction (3.5), particularly the Voronoi quadrature of Remark 3.3 and the ill-conditioned inversion of Λ; that is a correctness risk, not circularity.

Assumptions & free parameters 3 free parameters · 4 assumptions · 0 invented entities

The method's central claim rests on standard PDE existence theory, on the probabilistic representations of the LOT (proven in Appendix A.2), and on two unproven numerical steps: the Voronoi cell point-evaluation quadrature and the stability and regularization of the reduced algebraic system. No new physical entities are introduced; the free parameters are hand-chosen regularization and optimization hyperparameters.

free parameters (3)
  • rcond (truncated SVD relative threshold) = 1e-3 (Example 3), 1e-2 (Example 4)
    Used to regularize the ill-conditioned matrix Lambda when solving (3.5); the choice affects the numerical solution and no principled selection rule is given.
  • Adam learning rate alpha and initial guess f0 = alpha in {0.001, 0.01, 0.1, 1, 10}, f0 in {0.1, 1, 1000} depending on example
    Chosen per benchmark in Table 7; the number and type of solutions found in Example 8 depend on these and on restarts from nearby saddle points.
  • p-HiSD parameters (eta, tau, k) and f0 restarts = eta in {5, 20, 80}, tau = 10, k in {1,...,8}, f0 = 1000 or near previous saddle point
    Used to explore multiple solutions of the semilinear problem; the reported counts of 4 and 22 solutions are outcomes of hand-picked exploration strategies, not certified exhaustive enumerations.
assumptions (4)
  • standard math Existence and uniqueness of the weak solution to the generalized Dirichlet problem (3.1) on bounded Lipschitz domains for bounded measurable coefficients with A symmetric and strictly elliptic.
    Invoked in Section 3.1 to guarantee the target solution u in H1(D) exists; the paper cites Gilbarg and Trudinger [24, Theorem 8.3].
  • ad hoc to paper Quadrature approximation: integrands in the weak form are constant on each Voronoi cell and can be evaluated at the cell centers, requiring more regularity than the stated bounded measurable coefficients.
    Remark 3.3 explicitly assumes point values f(x_i), b(x_i), c(x_i) are well-defined and feasible; this is not implied by the hypotheses of Section 3.1 and the error is unquantified.
  • ad hoc to paper Stability and regularization of the algebraic system (3.5): the random matrix Lambda admits a pseudoinverse whose truncation yields a meaningful approximate source q_bar.
    The paper uses TSVD with hand-chosen thresholds (Table 7) because Lambda is ill-conditioned; no theory guarantees that the regularized solution converges to the true q or u.
  • domain assumption Harmless approximation of the distance-to-boundary function in WoS (Remark 2.16(i)).
    Relies on the authors' prior work [10] for quantifiable guarantees; the main body does not reprove this fact.

how reviews work

0 comments
Cite this review

Pith. "Pith review of Extended Walk-on-Spheres Algorithm for Linear and Nonlinear Elliptic Problems of Divergence-type." pith.science (2026). https://pith.science/paper/EQEQD22O

@misc{pith2026260807017,
  author       = {Pith},
  title        = {Pith review of: Extended Walk-on-Spheres Algorithm for Linear and Nonlinear Elliptic Problems of Divergence-type},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/EQEQD22O}},
  note         = {Machine review of arXiv:2608.07017}
}
read the original abstract

The Walk-on-Spheres algorithm, introduced by M. E. Muller in 1956, is a well known Monte Carlo method that leverages Brownian exit distributions from spheres to solve the Laplace equation with Dirichlet boundary conditions. Its mesh-free nature, robustness on complex geometries, favorable scaling with dimension, and intrinsic parallelism distinguish it from mesh-based solvers. However, its efficient applicability has been essentially limited to operators that admit explicit probabilistic exit laws, excluding most variable-coefficient and nonlinear elliptic operators. We propose a general framework that aims to overcome this limitation by using the classical Dirichlet Laplacian and harmonic extension as universal building blocks. Rather than seeking a custom stochastic representation for each operator, we employ Walk-on-Spheres to precompute a reusable numerical operator toolbox that approximates the inverse Dirichlet Laplacian, the harmonic extension operator, and their gradients. These precomputed operators are then used to represent candidate solutions and to transform arbitrary Dirichlet boundary value problems into finite-dimensional algebraic system for an unknown source term. Solving the resulting algebraic system and substituting back yields an approximate solution to the original PDE. Even more, for a general linear second order elliptic operator, the above mentioned precomputed toolbox can be directly used to obtain not just an approximation of a certain solution of the corresponding generalized Dirichlet problem, but an estimator of both the Green's integral operator and the elliptic measure operator. Numerical experiments on a range of benchmarks, including non-symmetric and anisotropic linear elliptic equations, semilinear and quasilinear problems, demonstrate the method's flexibility and efficiency.

Figures

Figures reproduced from arXiv: 2608.07017 by the authors.

Figure 1
Figure 1. for a visualization in a square domain): ◦ Let (xi)1≤i≤n ⊂ D be the prescribed centers for Wn and set Wi := {x ∈ D : ∥x − xi∥ ≤ ∥x − xj∥, ∀ 1 ≤ j ≤ n}, 1 ≤ i ≤ n. ◦ Let (yi)1≤i≤m ⊂ D be the prescribed centers for Vm and set Vi := {y ∈ D : ∥y − yi∥ ≤ ∥y − yj∥, ∀ 1 ≤ j ≤ m}, 1 ≤ i ≤ m. ◦ Let (zi)1≤i≤p ⊂ ∂D be the prescribed centers for Γp and set Γi := {z ∈ D : ∥z − zi∥ ≤ ∥z − zj∥, ∀ 1 ≤ j ≤ p}, 1 ≤ i ≤ p. As already … view at source ↗
Figure 2
Figure 2. Combined plots: Mean of pointwise absolute errors [PITH_FULL_IMAGE:figures/full_fig_p022_2.png] view at source ↗
Figure 3
Figure 3. Discretizations: (a) n = 20 equidistant interior points in (−1, 1) where the solution shall be ap￾proximated (orange) and p = 40 equidistant boundary points (blue); (b)-(d): FEM generated discretizations of D after fixing p ∈ {20, 30, 40} equidistant boundary points [PITH_FULL_IMAGE:figures/full_fig_p023_3.png] view at source ↗
Figures from the paper (17 more)
Figure 4
Figure 4. Figure 4: Comparison of pointwise absolute error |∂xuex−∂xufem| between the partial derivative of the exact solution and of the numerical solution obtained using FEM with P1 elements (top) and P2 elements (bottom), respectively (basic blue), and mean of pointwise absolute errors…
Figure 5
Figure 5. Figure 5: Exact solution uex (left), exact partial derivatives ∂xuex (middle), ∂yuex (right), plotted in (−1, 1)2 [PITH_FULL_IMAGE:figures/full_fig_p024_5.png]
Figure 6
Figure 6. Figure 6: Discretization of D with 106 points (left) as vertices of the triangulation (right) generated by FEM after fixing equidistant boundary points; points ordered by x, then y coordinate. 24 [PITH_FULL_IMAGE:figures/full_fig_p024_6.png]
Figure 7
Figure 7. Figure 7: Comparison of pointwise absolute error |uex −ufem| between the exact solution and the numerical solution obtained using FEM with P1 and P2 elements (light and basic blue, respectively), and mean of pointwise absolute errors 1 30 P30 i=1 |uex − ubi |, 30 estimates ubi o…
Figure 8
Figure 8. Figure 8: Numerical solution: 30 estimates ubi obtained with the proposed method (orange) together with their mean (black), using N = 103 Monte Carlo samples (top) and N = 104 (bottom), with M = 60 W oS steps. 25 [PITH_FULL_IMAGE:figures/full_fig_p025_8.png]
Figure 9
Figure 9. Figure 9: Comparison of pointwise absolute error |∂xuex−∂xufem| (top), |∂yuex−∂yufem| (bottom) between the partial derivatives of the exact solution and of the numerical solution obtained using FEM with P1 and P2 elements, respectively (light and basic blue), and pointwise absol…
Figure 10
Figure 10. Figure 10: Points (xi)1≤i≤n ⊂ D (a), (c) and points (zi)1≤i≤p ⊂ ∂D (b), (d), for D = B(0, 1) and D = (0, 1)2 . 3.3.2 Elliptic through ”solution operator” Example 3. Second order elliptic div(A(x)∇u) + b(x) · ∇u − cu = f. Consider ® div(A(x)∇u) + b(x) · ∇u − cu = f in D = B(0, 1)…
Figure 11
Figure 11. Figure 11: Source function f (left), exact solution uex (middle-left), exact partial derivatives ∂xuex and ∂yuex (middle-right and right), plotted in B(0, 1) [PITH_FULL_IMAGE:figures/full_fig_p027_11.png]
Figure 12
Figure 12. Figure 12: (Top) Exact solution uex (3 plots) and partial derivatives ∂xuex and ∂yuex; (Bottom) Approximate solution ub (3 plots) and approximate partial derivatives ∂dxu and ∂dyu, obtained with the proposed method. 27 [PITH_FULL_IMAGE:figures/full_fig_p027_12.png]
Figure 13
Figure 13. Figure 13: (Top) Exact solution uex (3 plots) and partial derivatives ∂xuex and ∂yuex; (Bottom) Approximate solution ub (3 plots) and approximate partial derivatives ∂dxu and ∂dyu, obtained with the proposed method [PITH_FULL_IMAGE:figures/full_fig_p028_13.png]
Figure 14
Figure 14. Figure 14: Exact solution uex (top) and approximate solution ub obtained with the proposed method (bottom) [PITH_FULL_IMAGE:figures/full_fig_p029_14.png]
Figure 15
Figure 15. Figure 15: Exact partial derivative ∂yuex of the solution (top) and its approximation ∂dyu obtained with the proposed method (bottom), for p−Laplace parameter p ∈ {1.5, 2, 5, 10, 50} (from left to right). 29 [PITH_FULL_IMAGE:figures/full_fig_p029_15.png]
Figure 16
Figure 16. Figure 16: Approximate solution ub obtained with the proposed method [PITH_FULL_IMAGE:figures/full_fig_p030_16.png]
Figure 17
Figure 17. Figure 17: Approximate of the partial derivative of the solution, [PITH_FULL_IMAGE:figures/full_fig_p030_17.png]
Figure 18
Figure 18. Figure 18: Approximate solution ub (3 plots) and approximate partial derivatives ∂dxu and ∂dyu, obtained with the proposed method. 3.3.6 Semilinear elliptic equation Example 8. Problem ∆u + u 2 = f studied by Breuer, McKenna, Plum [13]. Consider ® ∆u + u 2 = s sin(πx) sin(πy) in…
Figure 19
Figure 19. Figure 19: 4 essentially distinct approximate solutions of ∆ [PITH_FULL_IMAGE:figures/full_fig_p031_19.png]
Figure 20
Figure 20. Figure 20: 22 essentially distinct approximate solutions of ∆ [PITH_FULL_IMAGE:figures/full_fig_p032_20.png]

Discussion (0). Continue with ORCID to comment.

Reference graph

Works this paper leans on

57 extracted references · 56 canonical work pages

  1. [1]

    H. Aikawa. H¨ older continuity of the Dirichlet solution for a general domain.Bull. London Math. Soc., 34(6):691–702, 2002. 22

  2. [2]

    Allgower, S.-G

    E. Allgower, S.-G. Cruceanu, and S. Tavener. Application of numerical continuation to compute all solutions of semilinear elliptic equations.Adv. Geom., 9(3):371–400, 2009. 32

  3. [3]

    Arag´ on, J

    A. Arag´ on, J. Fern´ andez Bonder, and D. Rubio. Effective numerical computation ofp(x)-Laplace equations in 2D.Int. J. Comput. Math., 100(11):2111–2123, 2023. 30

  4. [4]

    Barth, C

    A. Barth, C. Schwab, and N. Zollinger. Multi-level monte carlo finite element method for elliptic pdes with stochastic coefficients.Numerische Mathematik, 119(1):123–161, 2011. 2

  5. [5]

    Bathe.Finite element procedures

    K.-J. Bathe.Finite element procedures. Klaus-Jurgen Bathe, 2006. 2

  6. [6]

    C. Beck, S. Becker, P. Grohs, N. Jaafari, and A. Jentzen. Solving the kolmogorov pde by means of deep learning.Journal of Scientific Computing, 88(3):73, 2021. 3

  7. [7]

    C. Beck, M. Hutzenthaler, and A. Jentzen. On nonlinear feynman–kac formulas for viscosity solutions of semilinear parabolic partial differential equations.Stochastics and Dynamics, 21(08):2150048, 2021. 2

  8. [8]

    Berner, M

    J. Berner, M. Dablander, and P. Grohs. Numerically solving parametric families of high-dimensional kolmogorov partial differential equations via deep learning.Advances in neural information processing systems, 33:16615–16627, 2020. 3

Show all 57 references
  1. [9]

    L. Beznea. The stochastic solution of the dirichlet problem and controlled convergence.Lecture Notes of Seminario Interdisciplinare di Matematica, 10, 2011. 7

  2. [10]

    Beznea, I

    L. Beznea, I. C ˆ ımpean, O. Lupa¸ scu-Stamate, I. Popescu, and A. Z˘ arnescu. From Monte Carlo to neural networks approximations of boundary value problems.https://arxiv.org/abs/2209.01432, 2022. 3, 6, 14

  3. [11]

    Bossy, N

    M. Bossy, N. Champagnat, S. Maire, and D. Talay. Probabilistic interpretation and random walkon spheres algorithms for the poisson-boltzmann equationin molecular dynamics.ESAIM: Mathematical Modelling and Numerical Analysis, 44(5):997–1048, 2010. 3

  4. [12]

    Braess.Finite elements: Theory, fast solvers, and applications in solid mechanics

    D. Braess.Finite elements: Theory, fast solvers, and applications in solid mechanics. Cambridge University Press, 2001. 2

  5. [13]

    Breuer, P

    B. Breuer, P. J. McKenna, and M. Plum. Multiple solutions for a semilinear boundary value problem: a computational multiplicity proof.J. Differential Equations, 195(1):243–269, 2003. 31, 32

  6. [14]

    Cheridito, H

    P. Cheridito, H. M. Soner, N. Touzi, and N. Victoir. Second-order backward stochastic differential equations and fully nonlinear parabolic pdes.Communications on Pure and Applied Mathematics: A Journal Issued by the Courant Institute of Mathematical Sciences, 60(7):1081–1110, 2007. 2

  7. [15]

    C ˆ ımpean, T

    I. C ˆ ımpean, T. Do, L. Gonon, A. Jentzen, and I. Popescu. Error analysis for the deep kolmogorov method.arXiv preprint arXiv:2508.17167, 2025. 3

  8. [16]

    C ˆ ımpean, I

    I. C ˆ ımpean, I. Popescu, and A. Zarnescu. Quantitative boundary h\”{o}lder estimates for the inho- mogeneous poisson problem through a probabilistic approach.arXiv preprint arXiv:2510.06906, 2025. 6, 22, 24

  9. [17]

    C ˆ ımpean, A

    I. C ˆ ımpean, A. Grecu, and L. Marin. Numerical spectral analysis of Cauchy-type inverse problems: A probabilistic approach.IMAJNA, 2026. 3, 6

  10. [18]

    Deaconu, S

    M. Deaconu, S. Herrmann, and S. Maire. The walk on moving spheres: a new tool for simulating brownian motion’s exit time from a domain.Mathematics and Computers in Simulation, 135:28–38, 2017. 39

  11. [19]

    Deaconu and A

    M. Deaconu and A. Lejay. A random walk on rectangles algorithm.Methodology and Computing in Applied Probability, 8(1):135–151, 2006. 3

  12. [20]

    J. L. Doob and J. Doob.Classical potential theory and its probabilistic counterpart, volume 262. Springer, 1984. 7

  13. [21]

    Fahim, N

    A. Fahim, N. Touzi, and X. Warin. A probabilistic numerical method for fully nonlinear parabolic PDEs.Ann. Appl. Probab., 21(4):1322–1364, 2011. 2

  14. [22]

    J. H. Friedman, J. L. Bentley, and R. A. Finkel. An algorithm for finding best matches in logarithmic expected time.ACM Transactions on Mathematical Software (TOMS), 3(3):209–226, 1977. 15

  15. [23]

    W. D. Gerhard. The probabilistic solution of the Dirichlet problem for 1 2 ∆ +⟨a,∇⟩+bwith singular coefficients.J. Theoret. Probab., 5(3):503–520, 1992. 34, 36

  16. [24]

    Gilbarg and N

    D. Gilbarg and N. S. Trudinger.Elliptic partial differential equations of second order. Springer-Verlag, Berlin, 2001. 5, 6, 7, 16, 19

  17. [25]

    M. B. Giles. Multilevel monte carlo methods.Acta numerica, 24:259–328, 2015. 2

  18. [26]

    Graham and D

    C. Graham and D. Talay.Stochastic simulation and Monte Carlo methods: mathematical foundations of stochastic simulation, volume 68. Springer Science & Business Media, 2013. 2

  19. [27]

    Grisvard.Elliptic problems in nonsmooth domains, volume 24 ofMonographs and Studies in Math- ematics

    P. Grisvard.Elliptic problems in nonsmooth domains, volume 24 ofMonographs and Studies in Math- ematics. Pitman (Advanced Publishing Program), Boston, MA, 1985. 5

  20. [28]

    Grisvard.Elliptic problems in nonsmooth domains, volume 69 ofClassics in Applied Mathematics

    P. Grisvard.Elliptic problems in nonsmooth domains, volume 69 ofClassics in Applied Mathematics. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2011. Reprint of the 1985 original [MR0775683], With a foreword by Susanne C. Brenner. 24

  21. [29]

    Gr¨ uter and K.-O

    M. Gr¨ uter and K.-O. Widman. The Green function for uniformly elliptic equations.Manuscripta Math., 37(3):303–342, 1982. 38

  22. [30]

    P. C. Hansen. The truncated SVD as a method for regularization.BIT, 27(4):534–553, 1987. 18

  23. [31]

    Heinrich

    S. Heinrich. Multilevel monte carlo methods. InInternational Conference on Large-Scale Scientific Computing, pages 58–67. Springer, 2001. 2

  24. [32]

    Huang, H

    B. Huang, H. Su, L. Zhang, and J. Zhao. Computing saddle points in stiff problems via a preconditioned high-index saddle dynamics method, 2026. 21, 31, 32, 33

  25. [33]

    Johnson.Numerical solution of partial differential equations by the finite element method

    C. Johnson.Numerical solution of partial differential equations by the finite element method. Courier Corporation, 2009. 2

  26. [34]

    D. P. Kingma and J. Ba. Adam: A method for stochastic optimization.https://arxiv.org/abs/1412.6980,

  27. [35]

    D. P. Kroese, T. Taimre, and Z. I. Botev.Handbook of monte carlo methods. John Wiley & Sons, 2013. 2

  28. [36]

    A. E. Kyprianou, A. Osojnik, and T. Shardlow. Unbiased ‘walk-on-spheres’ monte carlo methods for the fractional laplacian.IMA Journal of Numerical Analysis, 38(3):1550–1578, 2018. 3

  29. [37]

    L. Li, Y. Ye, and H. Li. An adaptive orthogonal basis method for computing multiple solutions of differential equations with polynomial nonlinearities.J. Sci. Comput., 100(1):Paper No. 11, 28, 2024. 32

  30. [38]

    W. C. H. McLean.Strongly elliptic systems and boundary integral equations. Cambridge university press, 2000. 7

  31. [39]

    Miller, R

    B. Miller, R. Sawhney, K. Crane, and I. Gkioulekas. Boundary value caching for walk on spheres.ACM Trans. Graph., 42(4), 2023. 3 40

  32. [40]

    M. E. Muller. Some continuous Monte Carlo methods for the Dirichlet problem.Ann. Math. Statist., 27:569–589, 1956. 3

  33. [41]

    Quarteroni and A

    A. Quarteroni and A. Valli.Numerical approximation of partial differential equations. Springer, 1994. 2

  34. [42]

    Sabelfeld and A

    K. Sabelfeld and A. Kireeva. A global random walk on spheres algorithm for calculating the solution and its derivatives of the drift-diffusion-reaction equations.Mathematical Methods in the Applied Sciences, 45(3):1420–1431, 2022. 3

  35. [43]

    K. K. Sabelfeld and D. Talay. Integral formulation of the boundary value problems and the method of random walk on spheres.Monte Carlo Methods Appl., 1(1):1–34, 1995. 3

  36. [44]

    Sakakibara and Y

    K. Sakakibara and Y. Shimizu. Numerical analysis of the Plateau problem by the method of fundamental solutions.J. Sci. Comput., 100(1):Paper No. 2, 20, 2024. 31

  37. [45]

    Sawhney and K

    R. Sawhney and K. Crane. Monte carlo geometry processing: A grid-free approach to pde-based methods on volumetric domains.ACM Transactions on Graphics, 39(4), 2020. 3, 6, 15

  38. [46]

    Sawhney and K

    R. Sawhney and K. Crane. Monte carlo geometry processing: a grid-free approach to pde-based methods on volumetric domains.ACM Trans. Graph., 39(4), aug 2020. 14

  39. [47]

    Sawhney, D

    R. Sawhney, D. Seyb, W. Jarosz, and K. Crane. Grid-free monte carlo for pdes with spatially varying coefficients.ACM Transactions on Graphics (TOG), 41(4):1–17, 2022. 3

  40. [48]

    Shalimova and K

    I. Shalimova and K. K. Sabelfeld. Random walk on ellipsoids method for solving elliptic and parabolic equations.Monte Carlo Methods & Applications, 26(4), 2020. 3

  41. [49]

    Shalimova and K

    I. Shalimova and K. K. Sabelfeld. Development and implementation of branching random walk on spheres algorithms for solving the 2D elastostatics Lam´ e equation.Monte Carlo Methods Appl., 29(1):79–93, 2023. 3

  42. [50]

    G. D. Smith.Numerical solution of partial differential equations: finite difference methods. Oxford university press, 1985. 2

  43. [51]

    J. C. Strikwerda.Finite difference schemes and partial differential equations. SIAM, 2004. 2

  44. [52]

    Sugimoto, T

    R. Sugimoto, T. Chen, Y. Jiang, C. Batty, and T. Hachisuka. A practical walk-on-boundary method for boundary value problems.ACM Trans. Graph., 42(4), 2023. 3

  45. [53]

    Talay.Probabilistic numerical methods for partial differential equations: elements of analysis

    D. Talay.Probabilistic numerical methods for partial differential equations: elements of analysis. Springer, 2006. 2

  46. [54]

    Thom´ ee.Galerkin finite element methods for parabolic problems, volume 25

    V. Thom´ ee.Galerkin finite element methods for parabolic problems, volume 25. Springer Science & Business Media, 2007. 2

  47. [55]

    Viswanath, H

    H. Viswanath, H. C. Nam, X. Deng, J. Berner, A. Anandkumar, and A. Bera. Operator learning using weak supervision from walk-on-spheres.arXiv preprint arXiv:2603.01193, 2026. 3

  48. [56]

    N. Wiener. The Dirichlet problem.J. Math. Phys., 3(3):127–146, 1924. 7

  49. [57]

    Z. Zou, Z. Wang, and G. E. Karniadakis. Learning and discovering multiple solutions using physics- informed neural networks with random initialization and deep ensemble.Proc. A, 481(2325):Paper No. 20250205, 28, 2025. 32 41

Pith tools

Reviewed August 10, 2026 · model on record in the stance chip above.