Random domain decomposition for parabolic PDEs on graphs

TL;DR

This paper introduces a stochastic domain decomposition method using the Random Batch Method for efficiently solving parabolic PDEs on complex graph structures, reducing computational resources while maintaining accuracy.

Contribution

The paper applies RBM directly at the PDE level on graphs without prior discretization, proving convergence and demonstrating significant computational and memory savings.

Findings

Convergence of the RBM scheme with first-order accuracy.

Substantial reductions in memory and computational time.

Method effectiveness across various time discretization schemes.

Abstract

The simulation of complex systems, such as gas transport in large pipeline networks, often involves solving PDEs posed on intricate graph structures. Such problems require considerable computational and memory resources. The Random Batch Method (RBM) has shown promise in addressing these challenges via stochastic decomposition techniques. In this paper, we apply the RBM at the PDE level for parabolic equations on graphs, without assuming any preliminary discretization in space or time. We consider a non-overlapping domain decomposition in which the PDE coefficients and source terms are randomized. We prove that the resulting RBM-based scheme converges, in the mean-square sense and uniformly in time, to the true PDE solution with first-order accuracy in the RBM step size. Numerical experiments confirm this convergence rate and demonstrate substantial reductions in both memory usage and computational time compared to solving on the full graph. Moreover, these advantages persist across different time discretization schemes.

Tables

	${ℬ}_1$	${ℬ}_2$	${ℬ}_3$	${ℬ}_4$	${ℬ}_5$	${ℬ}_6$
${𝒱}^j$	$\{v_1,v_2,v_3,v_4\}$	$\{v_4,v_5,v_7\}$	$\{v_4,v_6,v_7\}$	$\{v_7,v_8,v_9,v_{10}\}$	$\{v_1,v_2,v_3,v_4,v_5,v_6,v_7\}$	$\{v_4,v_5,v_6,v_7,v_8,v_9,v_{10}\}$
${𝒱}_0^j$	$\left\{v_3\right\}$	$\left\{v_5\right\}$	$\left\{v_6\right\}$	$\left\{v_8\right\}$	$\{v_3,v_4,v_5,v_6\}$	$\{v_4,v_5,v_6,v_7,v_8\}$
${𝒱}_b^j$	$\{v_1,v_2,v_4\}$	$\{v_4,v_7\}$	$\{v_4,v_7\}$	$\{v_7,v_9,v_{10}\}$	$\{v_1,v_2,v_7\}$	$\{v_3,v_9,v_{10}\}$
${𝒱}_b^j\cap {𝒱}_0$	$\left\{v_4\right\}$	$\{v_4,v_7\}$	$\{v_4,v_7\}$	$\left\{v_7\right\}$	$\left\{v_7\right\}$	$\left\{v_3\right\}$

Table 2. Table 1 : Coefficients α e \alpha_{e} and β e \beta_{e} for each edge.

Coefficients	$e_1$	$e_2$	$e_3$	$e_4$	$e_5$	$e_6$	$e_7$	$e_8$	$e_9$	$e_{10}$
${\alpha }_e$	10	-3	5	4	4	10	-3	5	-3	4
${\beta }_e$	-12	1	-7	-6	-6	-12	1	-7	1	-6

Table 3. Table 2 : Deterministic solver comparison on the full graph problem

Metric	IE	CN	$𝜽=\frac{\mathrm{𝟑}}{\mathrm{𝟒}}$	SIEM
Error	$3.764\times {10}^{-1}$	$7.327\times {10}^{-2}$	$3.747\times {10}^{-1}$	$3.689\times {10}^{-1}$
Execution Time (s)	$57.13$	$74.68$	$59.79$	$56.71$
Memory Usage (MB)	$200.14$	$208.00$	$200.36$	$199.93$

Table 4. Table 3 : Performance of RBM with IE and CN, using decomposition ( 39 ).

Metric	$h=0.002$		$h=0.004$		$h=0.006$
	IE	CN	IE	CN	IE	CN
${\mathrm{Error}}_1$	$3.896\times {10}^{-1}$	$3.866\times {10}^{-1}$	$7.904\times {10}^{-1}$	$7.910\times {10}^{-1}$	$1.105\times {10}^0$	$1.109\times {10}^0$
${\mathrm{Error}}_2$	$4.273\times {10}^{-2}$	$3.576\times {10}^{-2}$	$2.092\times {10}^{-1}$	$2.238\times {10}^{-1}$	$1.474\times {10}^{-1}$	$1.409\times {10}^{-1}$
Av. Execution Time (s)	$27.153$	$38.144$	$22.278$	$32.382$	$28.506$	$40.005$
Av. Memory Usage (MB)	$163.945$	$165.754$	$182.832$	$183.754$	$183.441$	$179.992$

Table 5. Table 4 : Performance RBM with θ \theta -method and SIEM, using decomposition ( 39 ).

Metric	$h=0.002$		$h=0.004$		$h=0.006$
	$𝜽$	SIEM	$𝜽$	SIEM	$𝜽$	SIEM
${\mathrm{Error}}_1$	$3.879\times {10}^{-1}$	$3.915\times {10}^{-1}$	$7.912\times {10}^{-1}$	$7.949\times {10}^{-1}$	$1.103\times {10}^0$	$1.107\times {10}^0$
${\mathrm{Error}}_2$	$3.741\times {10}^{-2}$	$4.162\times {10}^{-2}$	$2.092\times {10}^{-1}$	$2.113\times {10}^{-1}$	$1.408\times {10}^{-1}$	$1.422\times {10}^{-1}$
Av. Execution Time (s)	$27.559$	$24.273$	$22.854$	$22.378$	$29.251$	$27.849$
Av. Memory Usage (MB)	$162.992$	$162.570$	$182.121$	$185.500$	$183.512$	$183.043$

Table 6. Table 6 : Performance RBM with IE and CN, using decomposition ( 42 ).

Metric	$h=0.002$		$h=0.004$		$h=0.006$
	IE	CN	IE	CN	IE	CN
${\mathrm{Error}}_1$	$3.563\times {10}^{-1}$	$3.503\times {10}^{-1}$	$8.981\times {10}^{-1}$	$9.134\times {10}^{-1}$	$1.307\times {10}^0$	$1.309\times {10}^0$
${\mathrm{Error}}_2$	$7.238\times {10}^{-2}$	$4.725\times {10}^{-2}$	$1.796\times {10}^{-1}$	$1.892\times {10}^{-1}$	$2.656\times {10}^{-1}$	$2.462\times {10}^{-1}$
Av. Execution Time (s)	$23.142$	$33.798$	$23.672$	$35.799$	$23.068$	$34.339$
Av. Memory Usage (MB)	$106.582$	$107.004$	$106.672$	$102.113$	$101.703$	$100.848$

Table 7. Table 7 : Performance RBM with θ \theta -method and SIEM, using decomposition ( 42 ).

Metric	$h=0.002$		$h=0.004$		$h=0.006$
	$𝜽$	SIEM	$𝜽$	SIEM	$𝜽$	SIEM
${\mathrm{Error}}_1$	$3.551\times {10}^{-1}$	$3.589\times {10}^{-1}$	$8.952\times {10}^{-1}$	$8.997\times {10}^{-1}$	$1.305\times {10}^0$	$1.309\times {10}^0$
${\mathrm{Error}}_2$	$5.931\times {10}^{-2}$	$7.029\times {10}^{-2}$	$1.739\times {10}^{-1}$	$1.774\times {10}^{-1}$	$2.606\times {10}^{-1}$	$2.632\times {10}^{-1}$
Av. Execution Time (s)	$23.952$	$23.742$	$23.941$	$23.864$	$23.741$	$22.442$
Av. Memory Usage (MB)	$106.352$	$107.762$	$102.180$	$102.109$	$102.113$	$101.734$

Equations

$$\{\begin{array}{ll}{\partial }_t{y}_e-{\partial }_x(a_e{\partial }_x{y}_e)+b_e{\partial }_x{y}_e+p_e{y}_e=f_e,& (x,t)\in e\times (0,T),\\ y_{e^{\prime }}(v,t)=y_{e^{\prime \prime }}(v,t),& (v,t)\in {\mathcal{V}}_0\times (0,T),\forall e^{\prime },e^{\prime \prime }\in \mathcal{E}(v),\\ {\sum }_{e\in \mathcal{E}(v)}{a}_e(v){\partial }_x{y}_e(v,t)n_e(v)=0,& (v,t)\in {\mathcal{V}}_0\times (0,T),\\ y_e(v,t)=g_v(t),& (v,t)\in {\mathcal{V}}_b\times (0,T),e\in \mathcal{E}(v),\\ y_e(x,0)=y_e^0(x),& x\in e,\end{array}$$

$$n_e(v)=-1,\text{if }e=(v,\cdot ),n_e(v)=1,\text{if }e=(\cdot ,v),\text{and}{n}_e(v)=0\text{otherwise. }$$

$$\mathcal{E}=\bigcup _{i\in [M]}{\mathcal{E}}_s^i,\mathcal{V}=\bigcup _{i\in [M]}{\mathcal{V}}_s^i,\text{and}{\mathcal{E}}_s^i\cap {\mathcal{E}}_s^j=\varnothing ,\text{for }i\ne j.$$

$${\mathcal{E}}^j:=\bigcup _{i\in {\mathcal{B}}_j}{\mathcal{E}}_s^i,{\mathcal{V}}^j:=\bigcup _{i\in {\mathcal{B}}_j}{\mathcal{V}}_s^i,{\mathcal{V}}_0^j:=\{v\in {\mathcal{V}}^j:{\mathcal{E}}^j(v)=\mathcal{E}(v)\},{\mathcal{V}}_b^j:={\mathcal{V}}^j\setminus {\mathcal{V}}_0^j,$$

$${\chi }_{\psi }^{i,j}(x):=\psi (x)1_{S_j^i}(x),S_j^i:={\mathcal{G}}_s^i\setminus ({\mathcal{V}}_b^j\cap {\mathcal{V}}_0),$$

$${\zeta }_{\psi }^{{\omega }_k}(x):=\sum _{i\in {\mathcal{B}}_{{\omega }_k}}\frac{{\chi }_{\psi }^{i,{\omega }_k}(x)}{{\pi }_i},t\in I_k,\text{with}{\pi }_i=\sum _{j\in [N]:i\in {\mathcal{B}}_j}{p}_j,\text{for all }i\in [M],$$

$$\mathbb{E}[{\zeta }_{\psi }^{{\omega }_k}(x)]$$

$$=\sum _{i\in [M]}\psi (x)1_{{\mathcal{G}}_s^i}(x)=\psi (x).$$

$$k_t:=\min \{k\in [K]:t<kh\},\text{and}{\zeta }_{\psi }(x,t):=\sum _{i\in {\mathcal{B}}_{{\omega }_{k_t}}}\frac{{\chi }_{\psi }^{i,{\omega }_{k_t}}(x)}{{\pi }_i}.$$

$$\displaystyle\begin{cases}\partial_{t}z_{e}-\partial_{x}(\zeta_{a}\partial_{x}z_{e})+\zeta_{b}\partial_{x}z_{e}+\zeta_{p}z_{e}=\zeta_{1}f_{e},&(x,t)\in e\times(0,T),\\ z_{e^{\prime}}(v,t)=z_{e^{\prime\prime}}(v,t),&(v,t)\in\mathcal{V}_{0}\times(0,T),e^{\prime},e^{\prime\prime}\in\mathcal{E}(v),\\ \zeta_{1}(v,t)\bigl{(}z_{e}(v,t)-g_{v}(t)\bigr{)}=0,&(v,t)\in\mathcal{V}_{b}\times(0,T),e\in\mathcal{E}(v),\\ \sum_{e\in\mathcal{E}(v)}\zeta_{a}(v,t)\partial_{x}z_{e}(v,t)n_{e}(v)=0,&(v,t)\in\mathcal{V}_{0}\times(0,T),\\ z_{e}(x,0)=y_{e}^{0}(x),&x\in e.\end{cases}$$

$$a_e(x)>a_0>0\text{for a.a. }x\in e,\text{for every }e\in \mathcal{E}.$$

$$\sum _{e\in \mathcal{E}(v)}{b}_e(v)n_e(v)=0,\text{ for every }v\in {\mathcal{V}}_0,$$

$$\Vert y{\Vert }_{H^1(0,T;L^2(\mathcal{E}))}+$$

$$\Vert z{\Vert }_{H^1(0,T;L^2(\mathcal{E}))}+$$

$$\Lambda (t):=$$

$$+\Vert (p-{\zeta }_p)y{\Vert }_{L^2(\mathcal{E})}^2+\Vert (1-{\zeta }_1)f{\Vert }_{L^2(\mathcal{E})}^2],$$

$$\displaystyle\sup_{t\in[0,T]}\mathbb{E}\big{[}\|z(t)-y(t)\|^{2}_{L^{2}(\mathcal{E})}\big{]}\leq C\|\Lambda\|_{L^{1}(0,T)}^{1/2}h,$$

$$\underset{t\in [0,T]}{\sup }\Vert \mathbb{E}[z(t)]-y(t){\Vert }_{L^2(\mathcal{E})}^2\le C\Vert \Lambda {\Vert }_{L^1(0,T)}^{1/2}h.$$

$$\displaystyle\sup_{t\in[0,T]}\operatorname{Var}\big{[}\|z(t)-y(t)\|_{L^{2}(\mathcal{E})}\big{]}\leq C\|\Lambda\|_{L^{1}(0,T)}^{1/2}h.$$

$$\displaystyle\sup_{t\in[0,T]}\mathbb{P}\bigg{(}\|z(t)-y(t)\|_{L^{2}(\mathcal{E})}>\varepsilon\bigg{)}\leq\frac{C\|\Lambda\|_{L^{1}(0,T)}^{1/2}h}{\varepsilon^{2}}.$$

$${\zeta }_{\psi }^5(x)=\sum _{i\in {\mathcal{B}}_5}\frac{{\chi }_{\psi }^{i,5}}{{\pi }_i}=\frac{{\chi }_{\psi }^{1,5}(x)}{{\pi }_1}+\frac{{\chi }_{\psi }^{2,5}(x)}{{\pi }_2}+\frac{{\chi }_{\psi }^{3,5}(x)}{{\pi }_3},$$

$${\zeta }_{\psi }^6(x)=\sum _{i\in {\mathcal{B}}_6}\frac{{\chi }_{\psi }^{i,6}}{{\pi }_i}=\frac{{\chi }_{\psi }^{2,6}(x)}{{\pi }_2}+\frac{{\chi }_{\psi }^{3,6}(x)}{{\pi }_3}+\frac{{\chi }_{\psi }^{4,6}(x)}{{\pi }_4},$$

$${\zeta }_{\psi }^j(x)=\sum _{i\in {\mathcal{B}}_j}\frac{{\chi }_{\psi }^{i,j}}{{\pi }_i}=\frac{{\chi }_{\psi }^{j,j}(x)}{{\pi }_j},\text{for }j\in [4].$$

$${\pi }_1=\frac{1}{6}+\frac{1}{6}=\frac{1}{3},{\pi }_2=\frac{1}{6}+\frac{1}{6}+\frac{1}{6}=\frac{1}{2},{\pi }_3=\frac{1}{6}+\frac{1}{6}+\frac{1}{6}=\frac{1}{2},{\pi }_4=\frac{1}{6}+\frac{1}{6}=\frac{1}{3}.$$

$$\{\begin{array}{ll}{\partial }_t{\xi }_e^k=0,& (x,t)\in e\times I_k,e\in {\left({\mathcal{E}}^k\right)}^c,\\ {\xi }_{e^{\prime }}^k(v,t)={\xi }_{e^{\prime \prime }}^k(v,t),& v\in \left(({\mathcal{V}}^k{)}^c\setminus {\mathcal{V}}_b\right),\forall e^{\prime },e^{\prime \prime }\in \mathcal{E}(v),\\ {\xi }_e^k(x,t_{k-1})=z_e^{k-1}(x,t_{k-1}),& x\in e\in {\left({\mathcal{E}}^k\right)}^c,\end{array}$$

$$z_e^k=\{\begin{array}{ll}{\xi }_e^k(x,t_k),& \text{in }{\left({\mathcal{E}}^k\right)}^c,\\ {\eta }_e^k(x,t_k),& \text{in }{\mathcal{E}}^k,\end{array}$$

$$\{\begin{array}{ll}{\partial }_t{\eta }_e^k-{\partial }_x({\zeta }_a^k{\partial }_x{\eta }_e^k)+{\zeta }_b^k{\partial }_x{\eta }_e^k+{\zeta }_p^k{\eta }_e^k={\zeta }_1^k{f}_e,& (x,t)\in e\times I_k,e\in {\mathcal{E}}^k,\\ {\eta }_{e^{\prime }}^k(v,t)={\eta }_{e^{\prime \prime }}^k(v,t),& v\in {\mathcal{V}}_0^k,\forall e^{\prime },e^{\prime \prime }\in \mathcal{E}(v),\\ {\eta }_{e^{\prime }}^k(v,t)={\xi }_{e^{\prime \prime }}^k(v,t),& v\in ({\mathcal{V}}_b^k\cap {\mathcal{V}}_0),\forall e^{\prime },e^{\prime \prime }\in \mathcal{E}(v),\\ {\eta }_e^k(v,t)=g_v(t),& v\in ({\mathcal{V}}_b^k\cap {\mathcal{V}}_b),e\in \mathcal{E}(v),\\ {\sum }_{e\in \mathcal{E}(v)}{\zeta }_a^k(v){\partial }_x{\eta }_e^k(v)n_e(v)=0,& v\in {\mathcal{V}}_0^k,\\ {\eta }_e^k(x,t_{k-1})=z^{k-1}(x,t_{k-1}),& x\in e\in {\mathcal{E}}^k,\end{array}$$

$$\{\begin{array}{ll}{\partial }_t{u}_e-{\partial }_x(a_e{\partial }_x{u}_e)+b_e{\partial }_x{u}_e+p_e{u}_e={\vartheta }_e,& (x,t)\in e\times (0,T),e\in \mathcal{E},\\ u_{e^{\prime }}(v,t)=u_{e^{\prime \prime }}(v,t),& v\in {\mathcal{V}}_0,\forall e^{\prime },e^{\prime \prime }\in \mathcal{E}(v),\\ u_e(v,t)=0,& v\in {\mathcal{V}}_b,e\in \mathcal{E}(v),\\ {\sum }_{e\in \mathcal{E}(v)}{a}_e(v){\partial }_x{u}_e(v,t)n_e(v)=0,& v\in {\mathcal{V}}_0,\\ u_e(x,0)=u_e^0(x),& x\in e,\end{array}$$

$$\Vert u{\Vert }_{H^1(0,T;L^2(\mathcal{E}))}+\Vert u{\Vert }_{L^2(0,T;H^2(\mathcal{E}))}\le C\left(\Vert \vartheta {\Vert }_{L^2(0,T;L^2(\mathcal{E}))}+\Vert u^0{\Vert }_{H_0^1(\mathcal{E})}\right).$$

$$\displaystyle\begin{cases}Au:=-\partial_{x}(a_{e}\partial_{x}u_{e})+b_{e}\partial_{x}u_{e}+p_{e}u_{e},\quad\text{for every }e\in\mathcal{E},\\ D(A):=\left\{u\in H^{1}_{0}(\mathcal{E})\cap H_{pw}^{2}(\mathcal{E})\,\Big{|}\,\sum_{e\in\mathcal{E}(v)}a_{e}(v)\partial_{x}u_{e}(v)n_{e}(v)=0,\,\forall v\in\mathcal{V}_{0}\right\}.\end{cases}$$

$$u\in H^1(0,T;L^2(\mathcal{E}))\cap L^2(0,T;H^2(\mathcal{E}))\cap C([0,T];H_0^1(\mathcal{E})),$$

$$u\in H^1(0,T;L^2(\mathcal{E}))\cap L^2(0,T;H^2(\mathcal{E}))\cap C([0,T];H^1(\mathcal{E})),$$

Full text

Keywords: Partial differential equations on graphs, Random batch method, Domain decomposition.

1 Introduction

1.1 Motivation

Modeling and simulating physical phenomena such as gas transport or heat transfer in large, complex networks often leads to partial differential equations (PDEs) on graph structures [16, 24, 27, 14]. PDEs on graphs have been extensively studied for diffusion and transport over networked domains [24, 27, 14, 22, 6], and in singular-limit regimes [12, 2] and for control of parabolic and hyperbolic equations [8, 1, 11]. Nonetheless, the high complexity induced by the network size and topology can make classical solvers prohibitively expensive in both time and memory [30]. A common strategy to reduce memory usage is domain decomposition [23]. Within the context of graphs, both overlapping (Schwarz-type) and non-overlapping (Schur-complement) methods have been developed to reduce the solver cost; see [23, Chapter III] and [31] for overviews. However, these methods rely on an outer iteration in which one alternates between solving the PDE on each subdomain (subgraph) with Dirichlet (or Neumann) data borrowed from its neighbors and then updating those interface data before the next pass. Although each local solve is smaller, repeated synchronization of boundary data across all subdomains until global convergence can dominate runtime, especially as the number of subdomains grows.

On the other hand, Random Batch Method (RBM)—inspired by minibatch gradient descent—was introduced in [19, 20] to reduce the cost of simulating large-scale interacting particle systems by randomly grouping particles into small batches over short temporal intervals. RBM shares with domain decomposition the idea of splitting the problem into subproblems, but rather than solving every subproblem at each step, RBM randomly selects a single subproblem in each time window. Subsequent extensions of RBM encompass finite-dimensional dynamical systems [33, 21, 32], time-discrete PDEs and gradient flows [13, 10], and spatially discrete PDEs [17]. However, each of these variants remains tied to a particular mesh or time-discretization scheme, which restricts their applicability. In addition, the hallmark savings in memory and CPU time should emerge from the RBM principle itself, rather than from particular discretization choices. To address these issues, we introduce a continuous RBM at the PDE level, decoupling the RBM from any discretization to enable mesh-independent convergence proofs and intrinsic computational savings.

1.2 Contribution

Consider a connected finite metric graph $\mathcal{G}:=(\mathcal{E},\mathcal{V})$ with vertices $v\in\mathcal{V}$ and edges $e\in\mathcal{E}$. Denote by $\mathcal{V}_{b}$ and $\mathcal{V}_{0}$ the boundary and interior vertices of $\mathcal{G}$, and $\mathcal{E}(v)$ the edges adjacent to the vertex $v$. Let $T>0$, and over each edge $e\in\mathcal{E}$, we consider the parabolic equation

[TABLE]

where $y_{e}^{0}$ is the initial condition, $f_{e}$ the source term, $g_{v}$ the boundary data, and $a_{e}$, $b_{e}$, and $p_{e}$ are the diffusion, convection, and potential coefficients, respectively. Here $n_{e}:\mathcal{V}\rightarrow\{-1,0,1\}$ corresponds to the incidence vector, defined as

[TABLE]

Conditions $\eqref{eq:heat1}_{2}$–$\eqref{eq:heat1}_{3}$ correspond to continuity and the flux (Kirchhoff) condition at the interior vertices $\mathcal{V}_{0}$, while $\eqref{eq:heat1}_{4}$–$\eqref{eq:heat1}_{5}$ correspond to the boundary condition at $\mathcal{V}_{b}$ and the initial condition in $\mathcal{E}$.

Our main contributions are:

Continuous RBM on graphs. We introduce an RBM for PDEs defined on graphs, specifically for (1). Our method relies on the following steps: First, we introduce a non-overlapping decomposition of the graph into small subgraphs. Then, on each temporal window of length $h$, we randomly select one subgraph to evolve a localized version of (1) while keeping the solution constant on the rest of the graph (freezing). This freeze-and-evolve mechanism yields a randomized domain decomposition of the full-graph dynamics. We prove the resulting randomized system is well-posed and admits a unique solution in the appropriate function class. 2. 2.

Convergence estimates. We prove that the mean-square error between the RBM solution and the true solution decays like $O(h)$, and we derive probability bounds ensuring each realization of the random system stays close to the exact solution of (1) for sufficiently small $h$. 3. 3.

Numerical examples. We present numerical simulations to demonstrate the $O(h)$ convergence rate. The results show that the Continuous RBM achieves accuracy comparable to solving on the full graph. Furthermore, we observe reductions in CPU time and memory usage compared to full-graph solvers. These results are illustrated using the implicit Euler scheme, the $\theta$-method, the Crank–Nicolson scheme, and the semi-implicit method.

1.3 Related work

RBM was first proposed in [19, 20] for large-scale interacting particle systems, grouping particles into random batches to reduce complexity. In [33], it was extended to finite-dimensional linear systems—via decomposition of the system matrix—and to the linear-quadratic regulator. Subsequent work [21, 32] applies continuous RBM to model predictive control in finite dimensions.

For PDEs, existing RBM schemes typically rely on a preliminary time discretization. In [13], an RBM framework for nonlinear PDEs is introduced, viewed as an overlapping random domain decomposition. However, the authors first discretize the PDE in time using an implicit Euler scheme and then apply RBM, which forces the RBM parameter $h$ to coincide with the time discretization parameter $\Delta t$. They prove that the expected error between the randomized and exact solutions is $O(h)=O(\Delta t)$. Using similar ideas, [7] introduces a time-discretized scheme for stochastic evolution equations and establishes convergence. In [10], the same $O(h)$ rate is proved for finite-dimensional gradient flows and extended to infinite-dimensional flows, again, after temporal discretization via the minimizing movement scheme. In [17], a “discretize-RBM” approach is proposed for PDEs on graphs: this scheme, based on spatial discretization, yields expected-error convergence under requirements on both $h$ and the spatial mesh size $\Delta x$. RBM has also been extended to the mean-field limit of particle systems, notably to approximate kinetic equations such as the Landau and Fokker–Planck equations defined on the whole space [18, 9, 29].

It is important to note that the convergence estimates and frameworks from [13], [10], and [17] do not imply a continuous RBM setting, since their formulations do not decouple the RBM parameter $h$ from the respective discretization parameters. Consequently, our main result cannot be deduced from any prior work.

1.4 Outline

The remainder of this article is organized as follows. In section˜2, we present the domain decomposition setting and the continuous RBM dynamics. Subsequently, we state the main results. In section˜3, we show a decomposition example for a particular graph. In section˜4, we derive the continuous RBM system. Sections 5 and 6 contain proofs of the well-posedness and convergence results, respectively. Finally, in section˜7, we present numerical experiments that corroborate our theoretical convergence rate and the computational savings provided by the RBM.

2 Overview of the decomposition scheme and main results

In this section, we present the RBM formulation, which leads to the definition of the random system. Then, we present the well-posedness and convergence of the continuous RBM.

2.1 Continuous random batch method

We begin by introducing the graph decomposition, then the probabilistic setting, and finally, the random dynamics. Denote by $[M]=\{1,\dots,M\}$. An example of a decomposition satisfying the framework below is provided in Section˜3.

Graph decomposition. Let $M\in\mathbb{N}$ and consider a decomposition of $\mathcal{G}$ via a family of subgraphs $\{\mathcal{G}_{s}^{i}\}_{i\in[M]}:=\{(\mathcal{E}_{s}^{i},\mathcal{V}_{s}^{i})\}_{i\in[M]}$ that is pairwise edge-disjoint and covers $\mathcal{G}$, that is,

[TABLE]

In the following, we denote by $\mathcal{V}_{s,0}^{i}$ and $\mathcal{V}_{s,b}^{i}$ the interior and boundary nodes of $\mathcal{V}_{s}^{i}$, respectively. For every $v\in\mathcal{V}_{s,0}^{i}$, we denote by $\mathcal{E}_{s}^{i}(v):=\mathcal{E}_{s}^{i}\cap\mathcal{E}(v)$. We assume that $\mathcal{E}^{i}_{s}(v)=\mathcal{E}(v)$ for every $v\in\mathcal{V}_{s,0}^{i}$, that is, interior vertices in the subgraphs keep all adjacent edges.

Denote by $\mathcal{P}([M])$ the power set of $[M]$. Let $1\leq N\leq 2^{M}$ and $\mathcal{B}_{1},\dots,\mathcal{B}_{N}\subset\mathcal{P}([M])$ be nonempty sets such that $\bigcup_{j\in[N]}\mathcal{B}_{j}=[M].$ Then, for every $j\in[N]$, consider the set of edges and vertices:

[TABLE]

where $\mathcal{E}^{j}(v):=\mathcal{E}^{j}\cap\mathcal{E}(v)$. Hence, $\mathcal{V}_{0}^{j}$ collects the vertices whose incident edges are all in $\mathcal{E}^{j}$, while $\mathcal{V}_{b}^{j}$ contains the remaining (interface and boundary) vertices; see Section˜3 and Figure˜3. Then, we have the following structural assumption for the batches.

Assumption A1**.**

For every $v\in\mathcal{V}_{0}$ there exists a batch $\mathcal{B}_{\star}$ such that $v\in\mathcal{V}^{\star}_{0}$.

Localization function. Given $\psi:\mathcal{G}\rightarrow\mathbb{R}$ set

[TABLE]

for every $i\in[M]$ and $j\in[N]$, where $\mathbf{1}_{S^{i}_{j}}$ denotes the indicator function on $S^{i}_{j}$. By definition $\operatorname{supp}\chi_{\psi}^{i,j}\subset\mathcal{G}_{s}^{i}$ for every $i\in[M]$, and $\sum_{i\in[M]}\chi_{\psi}^{i,j}=\psi$ for every $j\in[N]$ and on each $e\in\mathcal{E}$.

Probabilistic framework: Let $K\in\mathbb{N}$ and consider the RBM parameter $h=T/(K+1)>0$. Define the time grid $t_{k}:=kh$ for $k\in\{0,\dots,K\}$, and the intervals $I_{k}:=(t_{k-1},t_{k})$ for each $k\in[K]$. Now, we introduce a family of i.i.d. random variables $\{\omega_{k}\}_{k\in[K]}$, each taking values in $[N]$. We denote by $p_{j}\in(0,1)$ the probability that $\omega_{k}$ equals $j\in[N]$, for each $k\in[K]$, which satisfy $\sum_{j\in[N]}p_{j}=1$.

Then, for a given function $\psi:\mathcal{G}\rightarrow\mathbb{R}$, we consider the piecewise constant in time random function

[TABLE]

for every $k\in[K]$, where $\pi_{i}$ is a normalization constant. We assume that $\pi_{i}>0$ for every $i\in[M]$. It is noteworthy that by construction, for every $x\in e\in\mathcal{E}$ and $k\in[K]$ we have that

[TABLE]

Therefore, the random function $\zeta_{\psi}^{\omega_{k}}$ is an unbiased estimator of $\psi$ on the edges for every $k\in[K]$. For $t\in(0,T)$ let us introduce

[TABLE]

Note that $\zeta_{\psi}(\cdot,t)=\zeta_{\psi}^{\omega_{k}}$ for every $t\in I_{k}$, thus $\zeta_{\psi}(\cdot,t)$ is piecewise constant in time.

Random system: Now, for every $e\in\mathcal{E}$, the continuous RBM system reads as:

[TABLE]

Note that the vertex–continuity condition in (7) is enforced on the whole graph, whereas the boundary and flux terms are switched on or off by the functions $\zeta_{1}$ and $\zeta_{a}$, depending on the subgraph that is active. A rigorous derivation of (7)—obtained by alternately freezing the state on an inactive subgraph and evolving it on the active subgraph at each interval $I_{k}=(t_{k-1},t_{k})$—is given in Section 4.

Remark 2.1** (On ˜A1).**

˜A1* ensures that every interior vertex in $\mathcal{V}_{0}$ has a positive probability of appearing as an interior vertex of some randomly chosen subgraph. Consequently, although the Kirchhoff flux condition is imposed only at the active interior vertices $\mathcal{V}_{0}^{k}$ within each interval, it is enforced in expectation across all $v\in\mathcal{V}_{0}$, thereby guaranteeing the consistency of the continuous RBM scheme.*

2.2 Main results

In the following, we consider the standard Lebesgue and Sobolev spaces on the graph; see Appendix˜A. Regarding the well-posedness of (1), we have the following theorem, whose proof is given in Section˜5.

Theorem 2.1** (Maximal regularity).**

Let the coefficients $(a,b,p)\in W^{1,\infty}_{pw}(\mathcal{E})\times L^{\infty}(\mathcal{E})\times L^{\infty}(\mathcal{E})$ with

[TABLE]

Consider $y^{0}\in L^{2}(\mathcal{E})$, $g\in H^{1}(0,T;l^{2}(\mathcal{V}_{b}))$, $f\in L^{2}(0,T;L^{2}(\mathcal{E}))$. Then, there exists a unique solution $y\in C([0,T];L^{2}(\mathcal{E}))\cap L^{2}(0,T;H^{1}(\mathcal{E}))$ of (1). Moreover, if $y^{0}\in H^{1}(\mathcal{E})$ and $g\in H^{1}(0,T;l^{2}(\mathcal{V}_{b}))$ satisfies the compatibility condition $y_{e}^{0}(v)=g_{v}(0)$ for each $v\in\mathcal{V}_{b}$, and the coefficient $b$ satisfies

[TABLE]

then the unique solution of (1) belongs to $C([0,T];H^{1}(\mathcal{E}))\cap L^{2}(0,T;H^{2}(\mathcal{E}))\cap H^{1}(0,T;L^{2}(\mathcal{E}))$, and satisfies the energy estimate

[TABLE]

with $C$ a positive constant depending on $T$ and the coefficients $a,\,b$ and $p$.

Remark 2.2**.**

Theorem˜2.1* establishes a maximal-regularity result for parabolic PDEs on graphs. To our knowledge, this maximal regularity has not been established previously. The proof proceeds in two steps. First, one shows that the homogeneous (non–self-adjoint) operator $A=-\partial_{x}(a\partial_{x}\,\cdot)+b\partial_{x}+p$ generates an analytic semigroup. Second, one identifies $D(A^{1/2})$ with a suitable interpolation space and then proves that this interpolation space coincides with $H_{0}^{1}(\mathcal{E})$.*

The following theorem addresses the well-posedness of the random system (7).

Theorem 2.2** (Well-posedness of (7)).**

Let $y^{0}\in H^{1}(\mathcal{E})$, $f\in L^{2}(0,T;L^{2}(\mathcal{E}))$, the coefficients $(a,b,p)\in W^{1,\infty}_{pw}(\mathcal{E})\times L^{\infty}(\mathcal{E})\times L^{\infty}(\mathcal{E})$ satisfying (8)-(9), and $g\in H^{1}(0,T;l^{2}(\mathcal{V}_{b}))$ satisfy the compatibility condition $y^{0}_{e}(v)=g_{v}(0)$ for each $v\in\mathcal{V}_{b}$. Then, there exists a unique solution $z$ of (7) belonging to the class $H^{1}(0,T;L^{2}(\mathcal{E}))\cap C([0,T];H^{1}(\mathcal{E}))$, and the following energy estimate holds:

[TABLE]

where $C>0$ is a constant depending on the coefficients and $T$, but independent of $h$.

Remark 2.3**.**

The proof of Theorem˜2.2 hinges on two facts. (i) By Theorem˜2.1, on each interval $I_{k}$, the dynamics restricted to the active subgraph admits a unique solution whose terminal value lies in $H^{1}$ on that subgraph. (ii) The vertex–continuity condition in (7) propagates this $H^{1}$ regularity to the full graph. Starting with $z_{0}\in H^{1}(\mathcal{E})$, we therefore obtain $z(t_{1})\in H^{1}(\mathcal{E})$ after the first step, and the argument can be re-applied on $I_{2}$. Induction over $k\in[K]$ yields global well-posedness.

To quantify the variance of the RBM, we introduce $\Lambda:(0,T)\rightarrow\mathbb{R}$ by

[TABLE]

where $y$ is the solution of (1), and $f$ and $(a,b,p)$ are the respective source term and coefficients. Note that when $N=1$, systems (1) and (7) coincide and $\Lambda\equiv 0$.

Now, we present our main convergence theorem.

Theorem 2.3** (Convergence).**

Let $y^{0}\in H^{1}(\mathcal{E})$, $f\in L^{2}(0,T;L^{2}(\mathcal{E}))$, the coefficients $(a,b,p)\in W^{1,\infty}_{pw}(\mathcal{E})\times L^{\infty}(\mathcal{E})\times L^{\infty}(\mathcal{E})$ satisfying (8)-(9), and $g\in H^{1}(0,T;l^{2}(\mathcal{V}_{b}))$ satisfy the compatibility condition $y^{0}_{e}(v)=g_{v}(0)$ for each $v\in\mathcal{V}_{b}$. Suppose that Assumption A1 holds. Let $y$ be the solution of (1) and $z$ the solution of (7). Then, $\Lambda\in L^{1}(0,T)$ is uniformly bounded with respect to $h$, and

[TABLE]

where $C>0$ depends on the problem data and $T$, but is independent of $h$.

Remark 2.4**.**

Note that the regularity assumptions in Theorem˜2.3 are the same as those in Theorem˜2.1 and Theorem˜2.2 when the maximal regularity is established. This is not merely a coincidence; the proof of Theorem˜2.3 relies on energy estimates, for which $H^{1}(0,T;L^{2}(\mathcal{E}))$ and $L^{2}(0,T;H^{2}_{pw}(\mathcal{E}))$ regularity of $z$ and $y$, respectively, is required. Then, after several energy estimates, and ensuring constants independent of $h$, the proof is concluded by an appropriate application of Grönwall’s lemma.

Finally, the following corollary is a direct consequence of Theorem 2.3.

Corollary 2.1**.**

Let $y^{0}\in H^{1}(\mathcal{E})$, $f\in L^{2}(0,T;L^{2}(\mathcal{E}))$, the coefficients $(a,b,p)\in W^{1,\infty}_{pw}(\mathcal{E})\times L^{\infty}(\mathcal{E})\times L^{\infty}(\mathcal{E})$ satisfying (8)-(9), and $g\in H^{1}(0,T;l^{2}(\mathcal{V}_{b}))$ satisfy the compatibility condition $y^{0}_{e}(v)=g_{v}(0)$ for each $v\in\mathcal{V}_{b}$. Let $y$ be the solution of (1) and $z$ the solution of (7). Then, the following statements hold:

• (i)

There exists a constant $C>0$, independent of $h>0$, such that

[TABLE]

• (ii)

There exists a constant $C>0$, independent of $h>0$, such that

[TABLE]

• (iii)

For all $\varepsilon>0$, there exists a constant $C>0$, independent of $h>0$, such that

[TABLE]

Proof.

Part (i) follows from Theorem 2.3 and Jensen’s inequality. Part (ii) is derived using the relation $\operatorname{Var}[X]=\mathbb{E}[X^{2}]-(\mathbb{E}[X])^{2}$ along with part (i). Finally, part (iii) is a consequence of Markov’s inequality applied to Theorem 2.3. ∎

Corollary˜2.1 ensures that the expected trajectory of the continuous RBM approaches the actual trajectory as $h\rightarrow 0$. Additionally, as $h$ decreases, the variance of the RBM solution diminishes. This is particularly significant because it implies that for sufficiently small $h$, a single realization of the random system serves as a representative trajectory of the expected trajectory, thereby closely approximating the solution of the parabolic equation.

3 A decomposition example

We now illustrate a non-overlapping decomposition of $\mathcal{G}$ following Section˜2.1.

1. Feasible decomposition: We will consider the graph and decomposition introduced in Figure˜1. By inspection, we note that the decomposition is pairwise edge-disjoint. Moreover, for each interior vertex $v\in\mathcal{V}_{s,0}^{i}$ one checks that $\mathcal{E}^{i}_{s}(v)=\mathcal{E}(v)$, for every $i\in[4]$. For example, at $v_{3}\in\mathcal{V}_{s,0}^{1}$ we have $\mathcal{E}^{1}_{s}(v_{3})=\{e_{1},e_{2},e_{3}\}=\mathcal{E}(v_{3})$. Then, within this decomposition, we introduce the following family of batches

• •

Option 1: $\mathcal{B}_{1}=\{1\}$, $\mathcal{B}_{2}=\{2\}$, $\mathcal{B}_{3}=\{3\}$, $\mathcal{B}_{4}=\{4\}$, $\mathcal{B}_{5}=\{1,2,3\}$, and $\mathcal{B}_{6}=\{2,3,4\}$; or

• •

Option 2: $\mathcal{B}_{1}=\{1\}$, $\mathcal{B}_{2}=\{2\}$, $\mathcal{B}_{3}=\{3\}$, $\mathcal{B}_{4}=\{4\}$, $\mathcal{B}_{5}=[4]$.

Note that neither $v_{4}$ nor $v_{7}$ lies in any $\mathcal{V}_{s,0}^{i}$. However, both option 1 and option 2, guarantee ˜A1. In fact, for option 1 we have that $v_{4}\in\mathcal{V}^{5}$ and $v_{7}\in\mathcal{V}^{6}$ (see table below), while for option 2 we have that $v_{4},\,v_{7}\in\mathcal{V}^{5}$. Note that, as long as there exists $\mathcal{B}_{*}=[M]$, ˜A1 automatically holds.

2. Localization function: Let us first analyze the sets of vertices defined in (2), associated with the batches from option 1. We have the following table:

Now, we can analyze the function $\zeta^{j}_{\psi}$ defined in (4) for $j\in[6]$. By definition:

[TABLE]

Figure 2 illustrates $\zeta^{1}_{\psi}$, $\zeta^{5}_{\psi}$, and $\zeta^{6}_{\psi}$ when $\psi\equiv 1$, and $p_{j}=1/6$ for $j\in[6]$. In this setting, we have that

[TABLE]

Note that $\zeta^{j}_{1}$ can be discontinuous at $\mathcal{V}_{0}$. Thus, even if $\psi$ is continuous across the vertices, $\zeta^{j}_{1}$ might not be. In general, if $\psi\in W^{1,p}_{pw}(\mathcal{E})$, then $\zeta^{j}_{\psi}\in W^{1,p}_{pw}(\mathcal{E})$ for every $j\in[N]$ and $p\in[1,\infty]$.

4 Derivation of the random system

Let us derive the random system (7). Initially, the system is defined separately for each time interval. In the following, for simplicity, for every $k\in[K]$ we denote by $\zeta^{k}=\zeta^{\omega_{k}}$ and $\mathcal{G}^{k}=(\mathcal{E}^{k},\mathcal{V}^{k})=(\mathcal{E}^{\omega_{k}},\mathcal{V}^{\omega_{k}})=\mathcal{G}^{\omega_{k}}$. Now, for each $k\in[K]$, we solve the ODE

[TABLE]

for every $t\in I_{k}$. Here, the function $z_{e}^{k}$ is defined by

[TABLE]

for every $k\in[K]$, with $z^{0}_{e}(x,0)=y_{e}^{0}(x)$, and $\eta$ is the solution of the parabolic equation

[TABLE]

for every $t\in I_{k}$ and $k\in[K]$. Here $\mathcal{V}^{k}_{b}\cap\mathcal{V}_{0}$ are the vertices on the interface between the chosen subgraph and its complement (see Figure˜3). The second condition in (16) corresponds to a continuity condition inside the chosen subgraph $\mathcal{G}^{k}$. The third and fourth conditions in (16) correspond to the boundary data of the problem. The fifth condition is the flux condition inside $\mathcal{V}_{0}^{k}$.

Remark 4.1**.**

On each interval $I_{k}$, equation (14) simply freezes the solution from the previous time interval on the inactive edges by setting $\xi^{k}_{e}\equiv z^{k-1}_{e}$ for $e\in(\mathcal{E}^{k})^{c},$ while (16) evolves the solution on the active subgraph $\mathcal{G}^{k}$, imposing continuity across the interface vertices and the prescribed Dirichlet data on $\mathcal{V}_{b}^{k}\cap\mathcal{V}_{b}$. Since no flux condition is enforced at the interface (to avoid overdetermination of the system), the glued solution $z^{k}$ generally does not belong to $H^{2}(\mathcal{E})$, although the continuity condition guarantees $z^{k}\in H^{1}(\mathcal{E})$ at each step.

One can verify that, on each $I_{k}$, (7) on $(\mathcal{G}^{k})^{c}$ reduces to (14) and (7) on $\mathcal{G}^{k}$ reduces to (16). System (7) will be primarily used to establish the convergence of the method, while systems (14)-(16) are more suitable (yet equivalent) for numerical implementation. In particular, system (14)-(16) can be solved by following Algorithm˜1.

5 Well-posedness of the parabolic equation on a graph

The proof of Theorem˜2.1 proceeds in two main steps. First, we show that the corresponding homogeneous system is well-posed by proving that its associated operator is m-accretive. By the Hille–Yosida theorem, this guarantees a unique solution in $C([0,T];L^{2}(\mathcal{E}))$. Second, we prove the maximal regularity by establishing that the operator associated with the homogeneous system is analytic. Hence, for any initial condition in the real interpolation space $(D(A),L^{2}(\mathcal{E}))_{1/2}$, there exists a unique solution in $H^{1}(0,T;L^{2}(\mathcal{E}))\cap L^{2}(0,T;D(A))$, where $D(A)$ is the operator domain. We then show that $(D(A),L^{2}(\mathcal{E}))_{1/2}\cong H_{0}^{1}(\mathcal{E})$. These technical results can be found in Appendix˜A. With these results established, the existence and uniqueness of the non-homogeneous system are followed by lifting the boundary data to the right-hand side of the PDE.

5.1 Well-posedness of the homogeneous system

We proceed by showing the well-posedness of the homogeneous Dirichlet boundary problem:

[TABLE]

for every $t\in(0,T)$, where the coefficients $(a,b,p)\in W^{1,\infty}_{pw}(\mathcal{E})\times L^{\infty}(\mathcal{E})\times L^{\infty}(\mathcal{E})$ satisfy (8) and (9). Then, the following theorem holds.

Lemma 5.1**.**

Let $u^{0}\in L^{2}(\mathcal{E})$ and $\vartheta\in L^{2}(0,T;L^{2}(\mathcal{E}))$. Then, (17) has a unique solution $u\in C([0,T];L^{2}(\mathcal{E}))$. Moreover, if $u^{0}\in H_{0}^{1}(\mathcal{E})$, then (17) has a unique solution $u\in C([0,T];H_{0}^{1}(\mathcal{E}))\cap H^{1}(0,T;L^{2}(\mathcal{E}))\cap L^{2}(0,T;H^{2}(\mathcal{E})),$ and there exists a constant $C>0$ such that the following energy estimate holds:

[TABLE]

Proof.

Define the linear operator $A:D(A)\subset L^{2}(\mathcal{E})\rightarrow L^{2}(\mathcal{E})$ by

[TABLE]

By Lemma˜B.1, $A$ is the infinitesimal generator of a strongly continuous semigroup on $L^{2}(\mathcal{E})$. Given $u^{0}\in L^{2}(\mathcal{E})$, $(a,b,p)\in W^{1,\infty}_{pw}(\mathcal{E})\times L^{\infty}(\mathcal{E})\times L^{\infty}(\mathcal{E})$ satisfying (8), and $\vartheta\in L^{2}(0,T;L^{2}(\mathcal{E}))$, [4, Proposition 3.1] ensures the existence of a unique solution $u$ of (17) in $C([0,T];L^{2}(\mathcal{E}))$.

Furthermore, due to Theorem˜B.1, we have that if $u^{0}\in H_{0}^{1}(\mathcal{E})$ and the coefficient $b$ satisfies (9), then the unique solution $u$ of (17) fulfills

[TABLE]

and the energy estimate (18) holds, concluding the result. ∎

5.2 Well-posedness of the non-homogeneous system

Proof of Theorem˜2.1.

Consider a quadratic function $\mu$ defined by $\mu_{e}(x,t)=\alpha_{e}(t)x^{2}+\beta_{e}(t)x+\gamma_{e}(t),$ satisfying the boundary conditions $\mu_{e}(v,t)=\partial_{x}\mu_{e}(v,t)=0$ for every $v\in\mathcal{V}_{0}$ and $\mu_{e}(v,t)=g_{v}(t)$ for every $v\in\mathcal{V}_{b}$. As $g_{v}\in H^{1}(0,T)$, choosing $\alpha,\beta,\gamma\in H^{1}(0,T)$ yields $\mu\in H^{1}(0,T;H^{2}(\mathcal{E}))$. Define $u$ to be the solution to the homogeneous problem (17) with the modified source term $\vartheta_{e}:=\partial_{x}(a_{e}\partial_{x}\mu_{e})-b_{e}\partial_{x}\mu_{e}-p_{e}\mu_{e}+f_{e},$ and the initial condition $u^{0}_{e}(x)=y^{0}_{e}(x)-\mu_{e}(x,0)$. Since $(a,b,p)\in W^{1,\infty}_{pw}(\mathcal{E})\times L^{\infty}(\mathcal{E})\times L^{\infty}(\mathcal{E})$, we have that $\vartheta_{e}\in L^{2}(0,T;L^{2}(\mathcal{E}))$. By Lemma˜5.1, there exists a unique solution $u\in C([0,T];L^{2}(\mathcal{E}))$ to (17). Moreover, if $y^{0}\in H^{1}(\mathcal{E})$ and the compatibility condition $u^{0}_{e}(v)=g_{v}(0)$ holds for every $v\in\mathcal{V}_{b}$, then $u$ satisfies

[TABLE]

and the energy estimate

[TABLE]

holds, where $C>0$ is a constant independent of $h$. Define $y=u+\mu$, then $y$ satisfies (1) with the initial condition $y_{e}(x,0)=u_{e}(x,0)+\mu_{e}(x,0)=y^{0}_{e}(x)$ and boundary conditions given by $g_{v}(t)$. For the uniqueness of $y$, let $y_{1}$ and $y_{2}$ be two solutions of (1), then $\tilde{u}:=y_{1}-y_{2}$ is solution of (17) with $\vartheta\equiv u_{0}\equiv 0$. By (18), we deduce the uniqueness. Moreover, using the regularity of $u$ and $\mu$, we obtain

[TABLE]

The energy estimate for $y$ follows from the estimates on $u$ and the properties of $\mu$. Specifically, applying (18) and utilizing the continuity of the right-inverse of the trace operator in $H_{0}^{1}(\mathcal{E})$ (see [26, Chapter I, Theorem 8.3]), we obtain

[TABLE]

where $C>0$ is a constant independent of $h$. ∎

5.3 Well-posedness of the continuous RBM

Let us prove Theorem˜2.2.

Proof of Theorem˜2.2.

We prove the well-posedness by induction on $k\in[K]$.

Base Case ($k=1$): Fix a subgraph $\mathcal{G}^{1}=(\mathcal{E}^{1},\mathcal{V}^{1})$ selected at random. Since $y^{0}\in H^{1}(\mathcal{E})$, it follows that z^{0}\in H^{1}\bigl{(}(\mathcal{E}^{1})^{c}\bigr{)}. Hence, $\xi^{1}(x,t)=z^{0}(x)$ is the unique solution of (14) on $I_{1}\times(\mathcal{E}^{1})^{c}$. Next, fixing $\xi^{1}(v,t)$ as boundary data in $(\mathcal{V}_{b}^{1}\cap\mathcal{V}_{0})$ and $g_{v}(t)$ on $(\mathcal{V}_{b}^{1}\cap\mathcal{V}_{b})$, Theorem˜2.1 ensures the existence of a unique solution

[TABLE]

to (16). Define

[TABLE]

Then z^{1}\in C\bigl{(}\overline{I_{1}};H^{1}(\mathcal{E})\bigr{)}\cap H^{1}\bigl{(}I_{1};L^{2}(\mathcal{E})\bigr{)} is the unique solution of (7) on $I_{1}\times\mathcal{E}$.

Inductive step: For some $n\in[K-1]$, suppose that there is a unique solution

[TABLE]

of (7) on $I_{n}\times\mathcal{E}$. We now prove the existence and uniqueness for $n+1$.

Pick a random subgraph $\mathcal{G}^{n+1}=(\mathcal{E}^{n+1},\mathcal{V}^{n+1})$. Let $z^{n}(x,t_{n})\in H^{1}(\mathcal{E})$ be the initial condition for (16) and (14) in their respective subgraphs. Then, $\xi^{n+1}(x,t)=z^{n}(x,t_{n})$ is the unique solution of (14) in $I_{n+1}\times(\mathcal{E}^{n+1})^{c}$. Considering $\xi^{n+1}(v,t)$ as boundary data in $\mathcal{V}_{b}^{n+1}\cap\mathcal{V}_{0}$ and $g_{v}(t)$ on $\mathcal{V}_{b}^{n+1}\cap\mathcal{V}_{b}$, Theorem˜2.1 guarantees the existence of unique solution of (16) in the class

[TABLE]

Define

[TABLE]

Then $z^{n+1}$ is the unique solution of (7) on $I_{n+1}\times\mathcal{E}$. By Theorem˜2.1,

[TABLE]

This completes the induction step. Thus, by induction, for each $k\in[K]$, there is a unique solution of (7) on $I_{k}\times\mathcal{E}$. Consequently, $z(x,t)=z^{k_{t}}(x,t),$ corresponds to the unique solution of (7) in $(0,T)\times\mathcal{E}$. Furthermore, z\in C\bigl{(}[0,T];H^{1}(\mathcal{E})\bigr{)}\cap H^{1}\bigl{(}0,T;L^{2}(\mathcal{E})\bigr{)}, and satisfies z\in L^{2}\bigl{(}I_{k};H^{2}(\mathcal{E}^{k})\bigr{)} for each $k\in[K]$. ∎

6 Convergence of the random scheme

In this section, we prove the convergence in expectation of the continuous RBM. We first show the following lemma.

Lemma 6.1**.**

Let $\Lambda$ be defined as in (12). Then, under the assumptions of Theorem˜2.1, we have $\Lambda\in L^{1}(0,T)$ and

[TABLE]

where $C>0$ depends on $(a,b,p)$, but it is independent of $h$.

Proof of Lemma˜6.1.

We observe that for each $k\in[K]$ and $t\in I_{k}$ we have that

[TABLE]

Denote by $\rho_{k}:=1-\sum_{i\in\mathcal{B}_{\omega_{k}}}\mathbf{1}_{S^{i}_{\omega_{k}}}/\pi_{i}$ the piece-wise function on space. Then, applying Hölder’s inequality we have

[TABLE]

for $t\in I_{k}$. On the other hand, let $\pi_{\ast}:=\min_{i}\pi_{i}>0$, then for each $e\in\mathcal{E}^{k}$ we have that $\mathbb{E}[|\rho_{k}|^{2}]\leq\frac{1}{\pi_{\ast}}-1=:C_{\pi_{\ast}}$ for every $k$. Therefore, the same arguments lead us to

[TABLE]

Since Theorem˜2.1 holds, we conclude by applying (10). ∎

Now, we continue with the convergence proof.

Proof of Theorem˜2.3.

Let $y$ and $z$ denote the solutions of (1) and (7), respectively. Fix an arbitrary index $k\in[K]$. For every $t\in I_{k}$, we have

[TABLE]

In the following, we split the proof into three steps. First, we estimate $P_{1}$, and then,the terms $P_{2},P_{3},$ and $P_{4}$. Finally, we will conclude the proof by using the regularity of $z$ on time and applying Grönwall’s inequality.

**Step 1. $P_{1}$ estimation. ** For simplicity, we denote by $\zeta_{\psi}=\zeta_{\psi}^{k}=\zeta^{\omega_{k}}_{\psi}$ for every $t\in I_{k}$. Now, for $P_{1}$, integrating by parts yields

[TABLE]

Let us analyze the boundary term. Since $\operatorname{supp}(\zeta_{a})\subset\mathcal{G}^{k}$ on $[t_{k-1},t_{k}]$, we have that

[TABLE]

where \Psi_{e}(v)=\zeta_{a}(v)\,\bigl{(}y_{e}-z_{e}\bigr{)}(v)\partial_{x}(y_{e}-z_{e})(v)\,n_{e}(v), and $\mathcal{V}^{k}$ is the set of vertices of $\mathcal{G}^{k}$. Observe that the set of vertices $\mathcal{V}^{k}$ can be decomposed as

[TABLE]

see Figure˜3. Therefore, since (22) and (23), we deduce that

[TABLE]

Now, by definition of $\zeta_{a}$, we have that

[TABLE]

Moreover, by construction

[TABLE]

Finally, it follows from (24), (25) and (26) that

[TABLE]

where we have used that $y,z\in H^{1}(\mathcal{E}^{k})$ and both satisfy Kirchhoff’s condition on $\mathcal{V}_{0}^{k}$. Consequently, $P_{1}$ becomes

[TABLE]

**Step 2. $P_{2}$, $P_{3}$ and $P_{4}$ estimation. ** Since for every $t\in(t_{k-1},t_{k})$ we have $\operatorname{supp}(\zeta_{b})\subset\mathcal{E}^{k}$, for $P_{2}$ we obtain that

[TABLE]

Since $a$ satisfies (8), $\zeta_{a}^{k}>0$ in $\mathcal{E}^{k}$. Thus, using Hölder’s inequality in (6) yields

[TABLE]

Now, applying Young’s inequality to (6) we get

[TABLE]

Finally, we estimate $P_{3}+P_{4}$; by applying Hölder’s inequality, we get

[TABLE]

Before proceeding, note that $z(t_{k-1})$ depends on the realization of the random variable $\omega_{k-1}$. By the independence of the sequence $\{\omega_{k}\}_{k=1}^{K}$, the value $z(t_{k-1})$ is independent of $\zeta_{\psi}^{k}$, which is defined through $\omega_{k}$. Thus, for any functions $\psi,\rho,\phi:\mathcal{G}\rightarrow\mathbb{R}$, we have

[TABLE]

Now, taking expectation in (6) and substituting (6), (6), (6), and (6), some terms vanish due to (6), we obtain

[TABLE]

where C_{\zeta_{a},\zeta_{b}}=\bigl{\|}\tfrac{\zeta_{b}}{\sqrt{\zeta_{a}}}\bigr{\|}_{L^{\infty}(\mathcal{E}^{k})}^{2} and $C_{p}=\|p\|_{L^{\infty}(\mathcal{E})}^{2}/\pi_{\ast}$ with $\pi_{\ast}:=\min_{i}\pi_{i}>0$. Then, due to the Cauchy-Schwarz inequality and the definition of $\zeta_{\psi}$, we have that

[TABLE]

being $a_{0}>0$ the ellipticity constant from (8). Now, using Hölder’s inequality in (6) for the $L^{2}(\mathcal{E})$ product and for the expectation, we get

[TABLE]

where $C=2(C_{a_{0},b}+C_{p})$, and $\Lambda$ is defined in (12).

Step 3. Error composition. Integrating (34) over $(t_{k-1},t)$ gives

[TABLE]

Now, composing the estimate (6) for every $k\in[k_{t}]$ yields

[TABLE]

Since $y(0)=z(0)$, applying the integral Grönwall lemma (see [15, Appendix B.k]) gives

[TABLE]

Applying the Cauchy-Schwarz inequality to the sum yields

[TABLE]

On the other hand, since $z\in H^{1}(0,T;L^{2}(\mathcal{E}))$, we have

[TABLE]

Substituting (6) into (6) we get

[TABLE]

Now, applying the energy estimate from Theorem 2.2, we obtain

[TABLE]

with $\hat{C}=C_{T}(\|g\|_{H^{1}(0,T;l^{2}(\mathcal{V}_{b}))}+\|f\|_{L^{2}(0,T;L^{2}(\mathcal{E}))}+\|y^{0}\|_{H^{1}(\mathcal{E})}).$ Completing the proof. ∎

7 Numerical simulations

In this section, numerical examples are presented to validate Theorem˜2.3. Moreover, we compare the memory usage and computational time required to solve (7) with those needed to solve (1) on the full graph. We consider the graph depicted in Figure˜1. Two examples are considered: the first uses the batches from option 2 in Section˜3, and the second uses those from option 1.

7.1 Example 1

We begin by introducing the finite element framework for systems (1) and (7).

Finite Element Setting. For the discretization of the parabolic equation (1) on the graph, we employ a finite element discretization using standard piecewise linear ($P1$) basis functions (see [17] for an illustrative introduction). Let $n\geq 1$, and denote by $\mathfrak{V}$ the finite element space spanned by the $P1$ basis functions $\{\phi_{i}\}_{i=1}^{n}$. The semi-discrete formulation of (1) consists of finding $y_{n}\in C^{1}(0,T;\mathfrak{V})$ that satisfies

[TABLE]

for every $t\in(0,T)$ and $i\in[n]$, with the condition $y_{n}(t,v)=g(t)$ for all $v\in\mathcal{V}_{b}$.

Next, consider the decomposition illustrated in Figure˜1. In this decomposition, the graph is partitioned into four subgraphs $\{\mathcal{G}_{i}\}_{i=1}^{4}$. To introduce the RBM framework, we consider the batches from option 2 in Section˜3:

[TABLE]

As is observed in Section˜3, this decomposition and the associated batches satisfy ˜A1. This batch selection imposes that either one subgraph is chosen or the full graph is chosen. For $K\in\mathbb{N}$, we introduce the random variables $\{\omega_{k}\}_{k=1}^{K}$, with associated probabilities $p_{j}=1/5$ for each $j\in[5]$. Note that $\pi_{i}=2/5$ for $i\in[4]$.

The implementation of (7) is carried out using finite elements together with Algorithm˜1. For each $k\in[K]$, on the randomly selected graph $\mathcal{E}^{k}$, the semi-discrete problem reads as:

[TABLE]

for every $t\in I_{k}$ and $i\in[n]$. The function $\eta_{n}$ is subject to the boundary conditions $\eta_{n}(v,t)=\xi_{n}(v,t)$ for every $v\in\mathcal{V}_{0}\cap\mathcal{V}_{b}^{k}$ and $\eta_{n}(v,t)=g_{v}(t)$ for every $v\in\mathcal{V}_{b}\cap\mathcal{V}_{b}^{k}$, with $\xi_{n}(x,t)=z_{n}(x,t_{k-1})$ for all $(x,t)\in e\times I_{k}$ for every $e\in(\mathcal{E}^{k})^{c}$, and $z_{n}(x,0)=y^{0}(x)$ for all $x\in e$ for every $e\in\mathcal{E}$.

From now on, each edge $e\in\mathcal{E}$ is identified with the interval $(0,1)$ (after fixing an orientation). Also, the PDE coefficients will be taken as

[TABLE]

To construct an explicit solution, the method of manufactured solutions (see [3, Chapter 12]) is used. Specifically, we prescribe a function $y$ on $\mathcal{G}\times(0,T)$ that satisfies the continuity and Kirchhoff conditions. Then, the corresponding source term $f$, boundary data $g$, and initial condition $y_{0}$ are computed as

[TABLE]

for every $e\in\mathcal{E}$ and $v\in\mathcal{V}_{b}$. We assume that

[TABLE]

and $v(t)=\sin(2\pi t)$, where the coefficients $\alpha_{e}$ and $\beta_{e}$ are given in Table˜1. These coefficients are chosen such that the exact solution satisfies both the continuity and Kirchhoff conditions.

Numerical Implementation. For the time discretization of (38) and (40), we compare the first-order methods: implicit Euler (IE), the $\theta$-method, and the semi-implicit Euler method (SIEM), as well as the second-order method: Crank–Nicolson (CN).

For each edge $e\in\mathcal{E}$, $100$ basis functions are introduced per edge; therefore, we use $1000$ for all edges. Moreover, to include the vertex (interior and boundary) conditions, we consider a basis function for each vertex. Since the graph in (1) has $10$ vertices, the finite element space $\mathfrak{V}$ is spanned by $n=1010$ basis functions. We set $T=1$ and use a time step of $\Delta t=0.002$. Using the previous schemes, the computational results obtained by solving (38) are summarized in Table˜2. Since the CN scheme is second-order, its error is smaller than that of the other methods, albeit at the expense of higher computational cost. In Table˜2, we define

[TABLE]

For the RBM, equation (7) is solved using $100$ basis functions per edge on the randomly selected graph, with a time step $\Delta t=0.002$. For $h\in\{0.002,0.004,0.006\}$, the results obtained by using the same time-schemes are shown in Tables˜3 and 4. To approximate the expected value, $20$ realizations are employed. The results for the case $h=0.002$ are shown in Figure˜4. In Table˜3, we define

[TABLE]

Further, in Table˜3, the average execution time and memory usage refer to the mean computational time and memory consumption, respectively. Note that the average execution time is approximately half that required to solve the equation on the full graph. Although RBM uses less memory, the difference is not substantial because, while the equation is solved on a small subgraph when one of the batches $\mathcal{B}_{1}$, $\mathcal{B}_{2}$, $\mathcal{B}_{3}$, or $\mathcal{B}_{4}$ is selected, the full graph is used when $\mathcal{B}_{5}$ is chosen. Furthermore, the best approximation is achieved for $h=0.002$, in accordance with Theorem˜2.3. Note that employing the CN scheme does not enhance accuracy in this context, which suggests that RBM converges more slowly than a second-order method. This observation does not contradict Theorem˜2.3, which guarantees an order of at least one.

We now analyze the limit $h\rightarrow 0$ for IE. Table 5 presents various values of $h$ alongside the corresponding $\text{Error}_{1}$ and $\text{Error}_{2}$ defined in (41). Clearly, the error decreases as $h$ diminishes. Furthermore, Figure˜5 compares these errors with the linear rate of convergence predicted by Theorem˜2.3 (and Corollary˜2.1).

7.2 Example 2

Now, for the graph decomposition described in Figure˜1, we consider the batches from option 1 introduced in Section˜3; that is,

[TABLE]

By setting the probabilities $p_{j}=1/6$ for $j\in\{1,\dots,6\}$, it follows that $\pi_{1}=\pi_{4}=1/3$ and $\pi_{2}=\pi_{3}=1/2$. With this configuration, (40) is solved using $100$ basis functions per edge, employing the previous time-schemes with a time step $\Delta t=0.002$, and for various values of $h\in\{0.002,0.004,0.006\}$. The corresponding results are presented in Tables˜6 and 7. Note that the accuracy obtained with the batches (42) is similar to that in Section˜7.1. However, in this instance, the equation is never solved on the full graph, resulting in lower average memory usage compared to that in Tables˜3 and 4. Moreover, RBM reduces memory consumption by approximately one-half relative to the full graph approach (see Table˜2).

Computational Setup. The previous numerical experiments were executed on a device featuring an AMD Ryzen 9 5900HS processor with Radeon Graphics, clocked at $3.30$ GHz, and equipped with $16.0$ GB of RAM (with $15.4$ GB usable). The experiments were implemented in Python, utilizing the NumPy and SciPy libraries.

8 Conclusion and open problems

In this work, we have introduced a fully continuous Random Batch Method for linear parabolic equations on graphs and proved in Theorem 2.3 that the mean square error between the exact solution and the RBM approximation decays like $O(h)$. Numerical experiments confirm this convergence rate and show that the method delivers savings in both CPU time and memory, and these gains are independent of the chosen time discretization scheme.

We conclude by highlighting several directions for future work:

$\bullet$ Optimal control. Extending RBM to optimal control entails introducing a randomized version of the usual cost functional, penalizing both state and control in the RBM-driven dynamics. In the distributed-control setting with an $\alpha$-convex (LQR) cost, one can adapt the energy-estimate proofs from [33]. For general Lipschitz (nonconvex) functionals, the Pontryagin principle and associated adjoint states must be invoked, and one might only expect convergence up to a constant reflecting the fact that small perturbations of the dynamics may drive to different local minimizers of the cost. Boundary-control problems pose additional difficulties, since Theorem 2.3 requires $H^{1}$-in-time boundary data; here, transposition methods or refined trace-space frameworks may be necessary.

$\bullet$ Abstract operator framework. One can ask whether RBM admits a formulation purely in operator-theoretic terms: let $(A,D(A))$ be a densely defined operator on a Hilbert space $H$ generating an analytic semigroup, and assume initial data in the real interpolation space $(D(A),H)_{1/2}$ to secure maximal regularity. In principle, the continuous RBM can then be recast by randomizing the action of $A$ (or its resolvent) on suitable subspaces, mirroring our edge-wise decomposition. However, proving convergence in this abstract setting will likely hinge on delicate issues such as the boundedness of the fractional power $A^{1/2}$ and the admissibility of any boundary (or control) operators in $H$, which may fail without additional structure.

$\bullet$ Optimal decomposition. As the graph grows in size and complexity, it becomes essential to choose both the subgraph partition and the batch-selection probabilities in a way that reflects the underlying PDE. One may borrow from spectral graph theory—using eigenvectors of the weighted Laplacian to suggest cuts—but standard spectral weights ignore PDE data (coefficients, sources, initial conditions). To make these partitions truly “physics-informed,” one can incorporate problem data into the Laplacian weights before computing the spectrum. An alternative is to train a graph-neural network on the dual graph, with problem data as input features, to predict which edges to retain in each batch. Finally, another possibility is that any candidate decomposition can be evaluated by its induced variance operator $\Lambda$, and therefore the goal is to minimize $\|\Lambda\|_{L^{1}(0,T)}$. This approach minimizes the multiplicative error in Theorem 2.3 and speeds up convergence. It remains open which of these strategies—spectral weighting, learned pruning, or variance-based optimization—offers the best trade-off between computational cost and accuracy.

Appendix A Functional spaces on graphs

Let $\mathcal{G}:=(\mathcal{E},\mathcal{V},l)$ be a metric graph with vertices $v\in\mathcal{V}$ and edges $e\in\mathcal{E}$, each having length $l_{e}\in(0,\infty)$. Define $\mathbb{N}_{0}:=\mathbb{N}\cup\{0\}$. Every edge $e\in\mathcal{E}$ is identified with the interval $[0,l_{e}]$, so that a function $y:\mathcal{G}\rightarrow\mathbb{R}$ corresponds to a collection $\{y_{e}\}_{e\in\mathcal{E}}$, where each $y_{e}:[0,l_{e}]\rightarrow\mathbb{R}$.

A.1 Sobolev spaces

We summarize the function spaces on graphs used in this work; see [27, Chapter 3] for additional details.

Piecewise spaces. For $m\in\mathbb{N}_{0}$ and $p\in[1,\infty]$, the piecewise Sobolev space is

[TABLE]

We write $L^{p}(\mathcal{E}):=W_{pw}^{0,p}(\mathcal{E})$ and $H_{pw}^{m}(\mathcal{E}):=W_{pw}^{m,2}(\mathcal{E})$. We endow $W_{pw}^{m,p}$ with the norm

[TABLE]

For $k\in\mathbb{N}_{0}\cup\{\infty\}$, define the piecewise continuous space

[TABLE]

Vertex-continuous spaces. By [15, Chap. 5.6, Thm. 6], we have that $W^{m,p}(0,l_{e})\subset C^{m-1,1-1/p}(0,l_{e})$ if $m,\,p>1$. Thus, $W_{pw}^{m,p}(\mathcal{E})\subset C^{m-1}_{pw}(\mathcal{E})$ if $m,\,p>1$. Thus, functions in $W_{pw}^{m,p}(\mathcal{E})$ are continuous on each edge $e$, but may be discontinuous at vertices. On the other hand, every function in $W^{m,p}(0,l_{e})$ with $m\geq 1$ and $p\in[1,\infty)$ has a well-defined trace ([15, Chapter 5.5, Theorem 1]). This last observation motivates the Sobolev spaces of functions continuous at the vertices:

[TABLE]

where $\mathcal{E}(v)$ is the set of all edges incident at $v$. For $p=2$, we write $H^{m}(\mathcal{E}):=W^{m,2}(\mathcal{E})$. If the graph has boundary vertices $\mathcal{V}_{b}\subset\mathcal{V}$, then

[TABLE]

It is noteworthy that the spaces $H^{m}(\mathcal{E})$ and $H_{0}^{m}(\mathcal{E})$ are Hilbert spaces since they are closed subspaces of the Hilbert space $\prod_{e\in\mathcal{E}}H^{m}(e)$ with (finite) constraints defined by continuous linear maps. Finally, for $k\in\mathbb{N}_{0}\cup\{\infty\}$, we define the space of $k$-times continuously differentiable functions on the graph by

[TABLE]

Here $\partial^{k}y$ denotes the $k$-th derivative of $y$. If $k=0$, we write $C(\mathcal{E}):=C^{0}(\mathcal{E})$. For simplicity, we retain the notation $C(\mathcal{E})$, even though functions in $C(\mathcal{E})$ belong to $C(\mathcal{G})$ and admit well-defined values at the vertices. For $k=\infty$, we let $C^{\infty}(\mathcal{E})$ be the set of infinitely differentiable functions satisfying vertex continuity of all derivatives. Moreover, $C^{\infty}_{c}(\mathcal{E})\subset C^{\infty}(\mathcal{E})$ denotes those functions with compact support in $\mathcal{E}$, and $C^{\infty}_{0}(\mathcal{E})\subset C^{\infty}(\mathcal{E})$ denotes those functions vanishing at every boundary vertex in $\mathcal{V}_{b}$.

Spaces on vertices. For $p\in[1,\infty)$, we consider $l^{p}(\mathcal{V})$, the space of real-valued functions defined on the set of vertices $\mathcal{V}$, given by

[TABLE]

In the remainder, we focus primarily on the space $l^{2}(\mathcal{V}_{b})$.

Bochner spaces. Finally, we also consider the standard definitions of Bochner spaces $L^{p}(0,T;X)$, continuous spaces $C([0,T];X)$, and Sobolev spaces $W^{m,p}(0,T;X)$ (with $H^{m}(0,T;X)=W^{m,2}(0,T;X)$), where $X$ will be taken as a functional space defined on a graph. See [15, Section 5.9.2] for a comprehensive introduction of Bochner spaces.

A.2 Density results

Regarding the relation between the previous spaces, we have the following lemmas.

Lemma A.1**.**

The space $C^{\infty}_{0}(\mathcal{E})$ is dense in $H_{0}^{1}(\mathcal{E})$.

Proof.

Take any $u\in H_{0}^{1}(\mathcal{E})$. By definition: $u_{e}\in H^{1}(e)$ for each $e\in\mathcal{E}$, $u_{e}(v)=0$ for every $v\in\mathcal{V}_{b}$, and $u_{e_{1}}(v)=u_{e_{2}}(v)$ for every $v\in\mathcal{V}_{0}$, and $e_{1},\,e_{2}\in\mathcal{E}(v)$. We will approximate $u$ in $H_{0}^{1}(\mathcal{E})$ by functions $\rho$ belonging to $C^{\infty}_{0}(\mathcal{E})$. We proceed in two steps:

Step 1. Approximation on each edge. Let $e\in\mathcal{E}$ be fixed. We have two cases:

-Case (a). $e$ meets a boundary vertex $v\in\mathcal{V}_{b}$. In this case $u_{e}\in H^{1}_{v}(e)$ where $H^{1}_{v}(e)$ denotes the subspace of $H^{1}(e)$ whose functions also vanish at $v$. Then it follows from [5, Section 8.3] that $C^{\infty}_{0(v)}(\overline{e})$ is dense in $H^{1}_{v}(e)$. This implies the existence of a sequence $\{\rho_{e}^{n}\}_{n\in\mathbb{N}}\subset C^{\infty}_{0(v)}(\overline{e})$ with $\rho_{e}^{n}\rightarrow u_{e}$ in $H^{1}(e)$. Here $C^{\infty}_{0(v)}(\overline{e})$ denotes smooth functions on $\overline{e}$ that vanish at $v$.

-Case (b). $e$ does not meet any boundary vertex. Then $u_{e}\in H^{1}(e)$, and by [15, Section 5.3.3], $C^{\infty}(\overline{e})$ is dense in $H^{1}(e)$. Hence there is a sequence $\{\rho_{e}^{n}\}_{n\in\mathbb{N}}\subset C^{\infty}(\overline{e})$ such that $\rho_{e}^{n}\rightarrow u_{e}$ in $H^{1}(e)$.

Step 2. Global smooth function. Let $v\in\mathcal{V}_{0}$ be an interior vertex, and suppose it is incident to $|E|\geq 2$ edges, labeled $e_{1},\dots,e_{k}\in\mathcal{E}(v)$. Since $u_{e_{i}}(v)$ is the same for all $i\in[|E|]$, each sequence $\{\rho_{e_{i}}^{n}\}_{n\in\mathbb{N}}$ converges to $u_{e_{i}}$ in $H^{1}(e_{i})$, and by Sobolev embedding $H^{1}(e_{i})\hookrightarrow C(\overline{e_{i}})$, we have $\rho_{e_{i}}^{n}(v)\rightarrow u(v)$ in the supremum norm for each $i$. Hence for large $n$, the differences $|\rho_{e_{i}}^{n}(v)-\rho_{e_{j}}^{n}(v)|$ become arbitrarily small, but may not vanish. To obtain a family $\{\tilde{\rho}_{e_{i}}^{n}\}_{n\in\mathbb{N}}$, whose values exactly coincide at $v$, for each $j\in\{2,\dots,|E|\}$, and every $n\in\mathbb{N}$ define

[TABLE]

where $\eta_{v}^{n}$ is a smooth cutoff supported in a small neighborhood of $v$ (on edge $e_{j}$) such that $\eta_{v}^{n}(v)=1$ and $\|\eta_{v}^{n}\|_{H^{1}(e_{j})}$ can be made arbitrarily small. This adjustment forces $\tilde{\rho}_{e_{1}}^{n}(v)=\tilde{\rho}_{e_{2}}^{n}(v)=\cdots=\tilde{\rho}_{e_{|E|}}^{n}(v),$ for every $n\in\mathbb{N}$. Finally, by defining $\psi^{n}\colon\mathcal{E}\rightarrow\mathbb{R}$ as $\psi^{n}|_{e}=\tilde{\rho}_{e}^{n}$, we have that $\{\psi^{n}\}_{n\in\mathbb{N}}\subset C^{\infty}_{0}(\mathcal{E})$. Moreover, by construction $\psi^{n}\rightarrow u$ in $H_{0}^{1}(\mathcal{E})$ as $n\rightarrow\infty$. Concluding the proof.

∎

Lemma A.2**.**

The space $H_{0}^{1}(\mathcal{E})$ is dense in $L^{2}(\mathcal{E})$.

Proof.

Since $C_{0}^{\infty}(e)$ is dense in $L^{2}(e)$ for each $e\in\mathcal{E}$, it follows that $C_{0}^{\infty}(\mathcal{E})$ is dense in $L^{2}(\mathcal{E})$. Because $C_{0}^{\infty}(\mathcal{E})\subset H_{0}^{1}(\mathcal{E})\subset L^{2}(\mathcal{E})$, we can take closure in the $L^{2}(\mathcal{E})$ topology, and deduce the result. ∎

Lemma A.3**.**

Let $H^{-1}(\mathcal{E})$ be the dual of $H_{0}^{1}(\mathcal{E})$. Then $L^{2}(\mathcal{E})$ is dense in $H^{-1}(\mathcal{E})$.

Proof.

For every $f\in L^{2}(\mathcal{E})$ consider $f\mapsto F_{f},$ where $F_{f}(u)=(f,u)_{L^{2}(\mathcal{E})}$ for all $u\in H_{0}^{1}(\mathcal{E}).$ Assume for contradiction that $L^{2}(\mathcal{E})$ is not dense in $H^{-1}(\mathcal{E})$. By the Hahn-Banach theorem, there exists a non-zero linear functional $\Phi\in(H^{-1}(\mathcal{E}))^{*}$ such that $\Phi(F_{f})=0,$ for all $f\in L^{2}(\mathcal{E}).$ By reflexivity, $(H^{-1}(\mathcal{E}))^{*}$ is isometrically isomorphic to $H_{0}^{1}(\mathcal{E})$; hence there exists a non-zero $u\in H_{0}^{1}(\mathcal{E})$ satisfying $\Phi(F)=F(u)$ for all $F\in H^{-1}(\mathcal{E}).$ In particular, for all $f\in L^{2}(\mathcal{E})$ we have

[TABLE]

This implies $u=0$ in $L^{2}(\mathcal{E})$ and consequently $u=0$ in $H_{0}^{1}(\mathcal{E})$, contradicting $u\neq 0$. Hence, $L^{2}(\mathcal{E})$ must be dense in $H^{-1}(\mathcal{E})$. ∎

Appendix B Properties of the main operator on the graph

B.1 M-accretivity and analyticity

Lemma B.1**.**

Consider $A:D(A)\subset L^{2}(\mathcal{E})\rightarrow L^{2}(\mathcal{E})$ defined in (19). Then $(A+\lambda I)$ is $M$-accretive for $\lambda>0$, and $A$ generates an analytic semigroup on $L^{2}(\mathcal{E})$.

Proof.

We will start first by showing the M-accretivity of the operator $A$.

1. M-accretivity. We split the argument into two parts:

(i) Proof of accretivity. Let $y\in D(A)$. By integrating by parts,

[TABLE]

Because $y(v)=0$ for every $v\in\mathcal{V}_{b}$ and

[TABLE]

the boundary term vanishes, giving

[TABLE]

Using Hölder’s and Young’s inequalities and the uniform ellipticity of $a$, we deduce

[TABLE]

Therefore, by choosing $\varepsilon=a_{0}/\|b\|_{L^{\infty}(\mathcal{E})}^{2}$, it follows that

[TABLE]

with $\lambda=\|b\|_{L^{\infty}(\mathcal{E})}^{2}/2\varepsilon+\|p\|_{L^{\infty}(\mathcal{E})}$. Hence, $A+\lambda I$ is accretive.

(ii) Surjectivity of $A+\lambda I$. To conclude the $M$-accretivity, it is enough to show that $\text{Range}(\lambda I-A)=L^{2}(\mathcal{E})$. Let us consider the bilinear form $\mathfrak{a}:H_{0}^{1}(\mathcal{E})\times H_{0}^{1}(\mathcal{E})\rightarrow\mathbb{R}$

[TABLE]

To verify that $\mathfrak{a}$ is continuous, let us apply Hölder’s and Young’s inequalities to get

[TABLE]

where $C_{p}>0$ is the Poincaré constant. Moreover, analogous to (45) the coercivity of $\mathfrak{a}_{\lambda}$ follows. Therefore, Lax–Milgram theorem guarantees that for every $f\in L^{2}(\mathcal{E})$ there exists a unique $y\in H_{0}^{1}(\mathcal{E})$ such that

[TABLE]

In particular, we deduce that

[TABLE]

Therefore, the distributional derivative of $a\partial_{x}y$ is given by

[TABLE]

Since $a\in W^{1,\infty}_{pw}(\mathcal{E})$, we deduce that $y\in H^{2}(\mathcal{E})\cap H_{0}^{1}(\mathcal{E})$. On the other hand, since $C^{\infty}_{0}(\mathcal{E})$ is densely embedded into $H^{1}_{0}(\mathcal{E})$, integrating by parts (47) we get

[TABLE]

for every $z\in H^{1}_{0}(\mathcal{E})$. Since $z(v)=0$ for every $v\in\mathcal{V}_{b}$, and $z_{e}(v)=z_{v}$ for every $v\in\mathcal{V}_{0}$, using (48), (49) is reduced to

[TABLE]

for all $z\in H^{1}_{0}(\mathcal{E})$. Consequently, $y\in D(A)$, and therefore, $A+\lambda I$ is M-accretive. In consequence of Lumer-Phillips theorem (see [4, Theorem 2.8]), $A+\lambda I$ generates a semigroup of contractions on $L^{2}(\mathcal{E})$, implying that $A$ generates a semigroup in $L^{2}(\mathcal{E})$.

2. Analytic semigroup. Now let us prove that this semigroup is analytic. For that let us consider the bilinear form $\mathfrak{b}:H_{0}^{1}(\mathcal{E})\times H_{0}^{1}(\mathcal{E})\rightarrow\mathbb{R}$ defined by

[TABLE]

Observe that the bilinear form $\mathfrak{b}$ is continuous. Moreover, in a similar way to (45), the bilinear form $\mathfrak{b}$ satisfies the Gårding inequality (or $H-V$ coercivity), i.e.,

[TABLE]

for some $\lambda,\,\beta>0$ ($\lambda$ can be taken as in (45) and $\beta=a_{0}/2$). Now, consider the operator $\mathfrak{\mathcal{B}}$ associated to $\mathfrak{b}$ defined by

[TABLE]

We will show that, in fact, $(\mathfrak{\mathcal{B}},D(\mathfrak{\mathcal{B}}))=(A,D(A))$. Let us characterize $D(\mathfrak{\mathcal{B}})$, for that, integrating by parts the first product of $\mathfrak{b}$, we observe

[TABLE]

Then, using the fact that $z_{e}(v)=0$ for every $v\in\mathcal{V}_{b}$, and applying Hölder’s inequality

[TABLE]

Therefore, for $z\mapsto\mathfrak{b}(y,z)$ to be continuous into $L^{2}(\mathcal{E})$, we require $y\in H^{2}_{pw}(\mathcal{E})$ and $\sum_{e\in\mathcal{E}(v)}a_{e}(v)\partial_{x}y_{e}(v)n_{e}(v)=0.$ Consequently, $y\in D(A)$ and we have that $D(\mathfrak{\mathcal{B}})\subset D(A)$, a completely analogous argument show that $D(A)\subset D(\mathfrak{\mathcal{B}})$. Thus, from (B.1), it follows that $\mathfrak{\mathcal{B}}=A$. Finally, due toLemma˜A.2 and Lemma˜A.3, $(H_{0}^{1}(\mathcal{E}),L^{2}(\mathcal{E}),H^{-1}(\mathcal{E}))$ defines a Gelfant triplet. It follows from [4, Theorem 2.12] that, due to (51), we can conclude that $A$ generates an analytic semigroup. ∎

B.2 Maximal regularity

Lemma B.2**.**

Let $A:D(A)\subset L^{2}(\mathcal{E})\rightarrow L^{2}(\mathcal{E})$ be the operator defined in (19). Assume $(a,b,p)\in W^{1,\infty}_{pw}(\mathcal{E})\times L^{\infty}(\mathcal{E})\times L^{\infty}(\mathcal{E})$ satisfy (8) and (9). Consider

[TABLE]

where $W(2,0;D(A),L^{2}(\mathcal{E})):=H^{1}(0,T;L^{2}(\mathcal{E}))\cap L^{2}(0,T;D(A))$. Then, $D_{A}(1/2,2)$ coincides with the space $H^{1}_{0}(\mathcal{E})$.

Proof.

The argument is carried out in two steps.

Step 1. Identification via interpolation. By [4, Proposition 4.3], one obtains D_{A}(1/2,2)\simeq\bigl{(}D(A),\,L^{2}(\mathcal{E})\bigr{)}_{1/2}, where the space on the right-hand side denotes the real interpolation space. To characterize this interpolation space, first note that $(A+\lambda I)$ is $M$-accretive by Lemma˜B.1, and due to (45), $(A+\lambda I)$ is injective. It follows that $(A+\lambda I)^{-1}$ is well-defined.

We now verify that $(A+\lambda I)^{-1}$ is a bounded operator on $L^{2}(\mathcal{E})$. Let $f\in L^{2}(\mathcal{E})$ and define $y:=(A+\lambda I)^{-1}f$. If $y=0$, the claim is immediate. Otherwise, applying Poincaré’s inequality, Hölder’s inequality, and (45) yields,

[TABLE]

for every $f\in L^{2}(\mathcal{E})$. Dividing by $\|y\|_{L^{2}(\mathcal{E})}$ shows that $(A+\lambda I)^{-1}$ is indeed bounded on $L^{2}(\mathcal{E})$. Hence, by [4, Proposition 6.1], we have that

[TABLE]

Step 2. Characterization of the domain of $A^{1/2}$. Let us characterize $A^{*}$, the adjoint operator of $A$. Consider $y\in D(A)$. Integrating by parts gives

[TABLE]

For the products to be well-defined in $L^{2}(\mathcal{E})$, one needs $z\in H^{2}_{pw}(\mathcal{E})$. We impose vertex conditions for $z$ such that the boundary term vanishes. Since $y\in H_{0}^{1}(\mathcal{E})$, setting

[TABLE]

the first boundary term vanishes. Further, because $y\in D(A)$, requiring $z_{e}\in H_{0}^{1}(\mathcal{E})$ makes the second boundary term vanishing as well. Lastly, the remaining boundary term vanishes since $z\in H_{0}^{1}(\mathcal{E})$ and (9) holds. Consequently, for all $y,z\in D(A)$,

[TABLE]

which shows $D(A^{*})=D(A)$. Applying [28, Theorem 8.2] (see also [25]) we deduce that D\bigl{(}A^{1/2}\bigr{)}=D\bigl{(}A^{*1/2}\bigr{)}=D(\mathfrak{b})=H_{0}^{1}(\mathcal{E}), where $\mathfrak{b}$ is the bilinear form associated with $A$ from (50). Finally, (53) completes the argument. ∎

Theorem B.1**.**

Let $A:D(A)\subset L^{2}(\mathcal{E})\rightarrow L^{2}(\mathcal{E})$ be the operator defined in (19). Then, if $u^{0}\in H^{1}_{0}(\mathcal{E})$ and $f\in L^{2}(0,T;L^{2}(\mathcal{E}))$, the equation

[TABLE]

has a unique solution $u$ in $H^{1}(0,T;L^{2}(\mathcal{E}))\cap L^{2}(0,T;D(A))\cap C([0,T];H_{0}^{1}(\mathcal{E})).$ Moreover, $u$ satisfies the following inequality

[TABLE]

Proof.

Due to Lemma˜B.1, $A$ generates an analytic semigroup. Then, we know from [4, Theorem 3.1] that the map

[TABLE]

is an isomorphism. On the other hand, due to lemma˜B.2, we have that $D_{A}(1/2,2)$ coincides with the space $H^{1}_{0}(\mathcal{E})$. Therefore, under the assumptions of Theorem˜B.1, we can guarantee that (54) has a unique solution $u$ in $H^{1}(0,T;L^{2}(\mathcal{E}))\cap L^{2}(0,T;D(A))$. Finally, using the analyticity of $A$, [4, Remark 4.2] ensures that $H^{1}(0,T;L^{2}(\mathcal{E}))\cap L^{2}(0,T;D(A))\subset C([0,T];H_{0}^{1}(\mathcal{E})),$ concluding the proof.

∎

Acknowledgments

M. Hernández has been funded by the Transregio 154 Project, Mathematical Modelling, Simulation, and Optimization Using the Example of Gas Networks of the DFG, project C07, the fellowship ‘ANID-DAAD bilateral agreement", and DAAD/CAPES grant 57703041 “Control and numerical analysis of complex systems", and DFG and NRF. Südkorea-NRF-DFG-2023 programme, grant 530756074.

The author wishes to express sincere gratitude to H. Meinlschmidt for valuable discussions on interpolation spaces, and to S. Zamorano for a careful and critical reading of the manuscript. The author is also grateful to E. Zuazua for his unconditional support and encouragement.