MeGA-MP: Metric Graph Advection Message Passing A Physics-Informed Message Passing Operator for Advection-Dominated Metric Graphs
Janine Strotherm∗
Bielefeld University
Luca Hermes∗
Bielefeld University
André Artelt
arXiv:2607.05167v1 [cs.LG] 6 Jul 2026
Bielefeld University
Barbara Hammer
Bielefeld University
Abstract Many real-world systems are organized as networks where spatio-temporal dynamics unfold along connections and not discretely between nodes. Examples include utility networks such as water distribution systems or gas networks, electrical grids, and traffic flow networks. Such systems are naturally modeled as metric graphs, where edges correspond to one-dimensional Euclidean subspaces connected at vertices. Metric graphs are independent of an underlying global Euclidean space, limiting direct application of typical PINNs and operator-learning methods. Especially transport dynamics like advection require a methodology able to capture antisymmetric and long-range dependencies on graphs, which is itself a challenge. We propose a novel physics-informed message passing operator that encodes linear advection on metric graphs as an inductive bias. In the purely advective setting, the operator provably recovers the exact dynamics up to a theoretically derived discretization error without any training. Combined with trainable components like MLPs, our message passing operator extends to realistic advection-reaction dynamics in water distribution systems, where we achieve superior performance compared to baselines and zero-shot generalization across different graph topologies.
1
Introduction
In the last years, physics-informed graph neural networks (GNNs) have shown impressive zero-shot generalization capabilities when addressing dynamical systems on conventional graphs (Thangamuthu et al., 2022). However, many real-world systems do not exhibit discrete node-to-node dynamics. Instead, dynamics unfold along connections between nodes. Examples include quantum dynamics, water distribution systems, gas networks, electrical grids, and traffic flow networks (Böttcher & Porter, 2024; Rossman et al., 2020). While a common approach is to discretize the domain and recover a conventional graph, metric graphs naturally capture the geometry of these systems. In a metric graph, each edge is associated with a one-dimensional domain – an interval of edge-dependent length – governed by edge dynamics, expressed as an ordinary differential equation (ODE) or, more commonly, a partial differential equation (PDE). In combination with suitable coupling conditions at nodes, such as continuity or conservation of mass or flux, this defines a coupled system of PDEs (Böttcher & Porter, 2024; Blechschmidt et al., 2025). We identify advection, i.e., the transport of a substance through the presence of a flow field, as the main driver for common edge dynamics on a variety of application domains, such as utility networks like water distribution systems (Rossman et al., 2020) and gas networks (Gugat & Herty, 2022), traffic flow networks ∗ Authors contributed equally.
1
(Gugat et al., 2005), telecommunication networks, and blood flow through vascular systems (Bressan et al., 2014). From the machine learning (ML) perspective, an intuitive approach is to model such dynamics with GNNs. However, spatio-temporal advection-dominated dynamics on metric graphs pose well-known challenges in the context of GNNs. First, the directional nature of transport requires antisymmetric message passing. Second, spatial dependencies can become global over time, requiring long-range spatial integration. Furthermore, sharp gradients and shocks are preserved and propagated rather than smoothed out as in diffusion-dominated settings. These high-frequency signals cannot be captured by GNNs with inherent over-smoothing characteristics. In this work, we make use of a suitable inductive bias for linear advection to overcome these challenges and model the dynamics accurately on node level. More precisely, we avoid computationally expensive sampling of function evaluations along each edge space by solving the linear advection PDE using the Method of Characteristics (MoC) and integrating the results in an iterative message passing architecture. This enables an efficient and exact realization of the linear advection dynamics while maintaining end-to-end adaptability if observational data is available. Contributions We propose metric graph advection message passing (MeGA-MP), a message passing neural network (MPNN) framework that provably models linear advection on a metric graph iteratively without requiring learning. More precisely, we rigorously derive MeGA-MP directly from the underlying physical priors and provide a theoretical bound on the discretization error (Section 3). We show how MeGA-MP can be used as a numerical solver for the spatio-temporal node forecasting task on a large-scale and realistic metric graph with advection edge dynamics and demonstrate its overall power by extending it with learning components to solve advection-reaction edge dynamics. Moreover, we investigate the effect of both the purely physics-driven and the learnable components of the model in two ablation studies (Section 4). To benchmark MeGA-MP against classical solvers, we use the 1-dimensional Euclidean domain as a special case of a metric graph (Section 5). A browser-version of our model with live control over inflow conditions is available online1 . Related Work MeGA-MP can serve as a standalone numerical solver for linear advection, expressed as a physics-informed MPNN-architecture. This formulation also allows a straight-forward extension with learnable components into a trainable MPNN to solve more complex, advection-dominated problems on metric graphs. This positions the work at the intersection of numerical PDE solvers for metric graphs, physicsinformed ML, and ML for metric graphs. In this order, we will now discuss relevant existing approaches. Numerical methods for metric graphs largely are specifically designed for quantum graphs, i.e., metric graphs with edge dynamics and coupling conditions associated with the self-adjoint Laplace operator (Arioli & Benzi, 2018; Brio et al., 2022; Böttcher & Porter, 2024). Such works benefit from the rich existing theory on the Laplace operator (Berkolaiko, 2017; Berkolaiko & Kuchment, 2013; Kuchment, 2003; 2005; 2008). In particular, the Laplace operator on a metric graph with continuity and Kirchhoff-Neumann coupling conditions has an orthogonal eigenbasis, which can be utilized in different ways, e.g., by eigenfunction-based spectral solutions (Arioli & Benzi, 2018; Brio et al., 2022; Böttcher & Porter, 2024). Besides quantum graphs, numeric approaches also exist for specific application domains, such as water distribution systems (Rossman et al., 2020; Shang et al., 2023), gas networks (Gugat & Herty, 2022) or traffic networks (Bressan et al., 2014; Gugat et al., 2005). For a broader overview of application domains and related numeric approaches, we refer to Böttcher & Porter (2024). However, numerical approaches typically require discretization of the continuous edge domains, resulting in large numbers of degrees of freedom and increasing computational and memory costs as the spatial resolution grows. In addition, numerical solvers can become computationally demanding when high spatial or temporal resolution is required. These challenges have motivated the development of ML-based approaches for solving dynamical systems. To this end, Raissi et al. (2019) introduced physics-informed neural networks (PINNs), a framework for neural networks (NNs) that learns the unknown solution of a given PDE on a continuous time and space domain. Since PDE solutions depend on configurations like initial conditions, boundary conditions, and control signals, PINNs typically learn the solution to one such configuration. Neural operators (Kovachki 1 https://anonymous.4open.science/w/tmp-preprint-BB4F/
2
et al., 2023) address this limitation by learning a map from functions that configure a PDE to its solution. Prominent examples include DeepONet by Lu et al. (2019; 2021) and a variation of neural operators by Li et al. (2020a;b;c; 2024). The majority of these methods are developed for Euclidean domains and do not directly transfer to metric graphs. This issue was recently addressed in different works that apply domain decomposition techniques. This transforms the learning problem into solving dynamics on multiple Euclidean domains – the edge spaces – while constraining correct coupling at nodes. Blechschmidt et al. (2022) use edge-specific PINNs in such a setup to solve drift-diffusion equations on metric graphs. Since PDE parametrizations are – in line with typical PINN setups – not a model input, the approach cannot generalize to changing conditions. Laczkó et al. (2025) apply PINNs in a Richardson iterative scheme to solve the time-independent Schrödinger equation on metric graphs. This allows the graph structure to be excluded from model training, making the approach graph agnostic. However, generalization to changing PDE parameters is limited. Blechschmidt et al. (2025) address this by utilizing DeepONets and building a lego-like graph-agnostic architecture to solve drift-diffusion equations on metric graphs. However, this framework requires a node-level optimization procedure to determine interface boundary values for each new problem instance, introducing additional computational overhead at graph interfaces. Moreover, all edge dynamics considered in these works are diffusion-dominated and thus accompanied by the aforementioned Kirchhoff-Neumann conditions – a coupling condition fundamentally different to the ones required for advection-dominated edge dynamics. An orthogonal line of work has shown that the discretization of PDEs in time and space yield powerful physics-motivated MPNNs suitable for conventional tasks on graphs, such as node classification (Chamberlain et al., 2021; Eliasof et al., 2021; Rusch et al., 2022; Choi et al., 2023; Eliasof et al., 2024), but they do not directly address the continuous edge domains characteristic of metric graphs. The reason is the difference in geometry between a conventional graph with discrete node-to-node relations vs. a metric graph with continuous edge spaces. See Appendix A for a detailed discussion, which we recommend reading after the next section. The existing ML-based methods discussed above either learn solutions for individual PDE instances (PINNs), rely on interface optimization procedures (DeepONet-based approaches), or are designed for conventional graphs without continuous edge domains. In contrast, MeGA-MP directly embeds the numerical structure of linear advection on metric graphs into an MPNN architecture, providing a graph-native framework that can be extended with learnable components while preserving the underlying physical coupling mechanisms.
2
Background
The common use case for graphs in dynamical system modeling is as a means of discretizing a continuous Euclidean domain Ω ⊂ Rn , with nodes representing spatial locations and edges encoding proximity in Ω. For example, a grid is a graph that approximates a multidimensional Euclidean space, preserves the structure of the continuum, and is commonly used when solving PDEs with numerical schemes (Brandstetter et al., 2022). The same concept has also been adopted to formulate physics-motivated MPNN architectures such as ADR-GNN by Eliasof et al. (2024). Such a setup is not practical in cases where dynamical systems govern the evolution of signals in networks such as water distribution systems, gas networks, electrical grids or traffic flow networks. Instead, metric graphs can be employed, as they naturally capture the topology of the original domain. In particular, the edges euv ∈ E of a metric graph G = (V, E) are parametrized by lengths luv > 0 and, potentially, capacities αuv > 0, and are associated with one-dimensional Euclidean subdomains Ωuv = (0, luv ) aligned with the edge direction. Here, the graph G describes the system Ω := ⊔euv ∈E Ωuv of such subdomains coupled at nodes V , which correspond to the boundaries ∂Ωuv = {0, luv } of the edge spaces Ωuv . We visualize the difference between metric graphs and graphs that discretize an underlying continuous space using the example of advection in Figure 10 in Appendix A. Dynamical Systems on Metric Graphs In this work, we focus on metric graphs with associated dynamics that govern the spatio-temporal evolution of signals along the graph. We model this system as a finite, connected, and symmetric directed graph G = (V, E) together with the following properties which link 3
signals cv : R → R≥0 , t 7→ cv (t) and cuv : R × Ωuv → R≥0 , (t, z) 7→ cuv (t, z) over time R, at each node v ∈ V and along the edge space Ωuv = (0, luv ) of each edge euv ∈ E, respectively. 1. The edge dynamics are usually governed by a PDE that describes the change of the signal cuv over time R and along the edge space Ωuv of an edge euv ∈ E, that is, ∂cuv (t, z) = f (t, z, cuv , νuv ) ∂t
(1)
for all t ∈ R, z ∈ Ωuv = (0, luv ) and known functions or coefficients νuv that parametrize the edge-wise PDE function f . 2. The coupling conditions define how signals propagate through nodes in between adjacent edge spaces. Specifically, a node v ∈ V can be understood as the boundary z = luv or z = 0 of multiple adjacent edge spaces Ωuv or Ωvu aligned with the direction of the corresponding incoming or outgoing edge euv ∈ E or evu ∈ E, respectively. Advection on Metric Graphs PDE
We model the edge dynamics, i.e., Equation (1), as the linear advection ∂cuv (t, z) ∂cuv (t, z) = −νuv (t) ∂t ∂z
(2)
which describes how the signal cuv is transported in the presence of an incompressible, i.e., divergence-free, flow field νuv : R → R over time R and along the edge space Ωuv . The flow field additionally satisfies νvu (t) = −νuv (t)
(3)
and we use the convention that the sign of a flow velocity νuv (t) is positive (negative) if the direction of the corresponding edge euv ∈ E aligns (does not align) with the direction of the physical flow, i.e., if the physical flow is from u ∈ N (v) to v ∈ V (from v ∈ V to u ∈ N (v)). The coupling conditions for advection-dominated dynamics describe how signals propagate between different edge spaces. On the one hand, the mixing at nodes P u∈N− (v,t) quv (t) · cuv (t, luv ) P cv (t) = (4) u∈N− (v,t) quv (t) describes how signals propagate from edges into nodes. On the other hand, the inflow continuity condition cuv (t, 0) = cu (t) for all u ∈ N− (v, t)
(5)
describes how signals propagate from nodes into edges and can be seen as the local inflow boundary condition of the edge dynamics (Eq. (2)). In Equation (4), quv (t) defines the flow rate quv = αuv νuv : R → R from (temporal) inflow neighbors u ∈ N− (v, t) = {u ∈ N (v) | νuv (t) > 0} of v ∈ V at time t ∈ R, that are neighbors u ∈ N (v) where at time t, the physical flow field is directed from u to v. For details, we refer to Appendix A. Remark 2.1 (Global vs. local boundary conditions). Note that Equation (4) and (5) define local boundary conditions of the edge dynamics, i.e., Equation (2). In contrast, global boundary conditions of the overall dynamical system of the metric graph G correspond to Dirichlet boundary nodes v ∈ Vb ⊊ V for which the signal cv (t) is known at all times t ∈ R, also see Definition 3.1.
3
Method: Metric Graph Advection Message Passing
In this section, we cast the advection edge dynamics (Eq. (2)) and corresponding coupling conditions (Eq. (4) and (5)) into a novel message-passing framework that only requires knowledge of the system states 4
c = (cv )v∈V at nodes and the flow field ν = (νe )e∈E along edges. This is based on the observation that Equation (4) is a special case of message-passing (cf. Appendix B) in the form of X cv (t) = wuv (t) ϕ(t, euv , c, ν) with u∈N (v)
(6)
ReLU(quv (t)) u∈N (v) ReLU(quv (t))
wuv (t) = P
and
ϕ(t, euv , c, ν) = cuv (t, luv ),
where wuv (t) are the normalized weights that select inflows quv (t) > 0 to v ∈ V only and scale the messages generated by the message function ϕ. In the current formulation, this requires knowledge of the signal cuv (t, luv ) at the outflow boundary of each edge space Ωuv = (0, luv ). However, instead of explicitly solving the dynamics along the edge spaces Ω = ⊔euv ∈E Ωuv , we derive a representation of ϕ that only relies on signals c at nodes V . This leads to the following general problem we are trying to solve. Definition 3.1 (Spatio-temporal node forecasting with global boundary conditions). Let Vb ⊊ V be a subset of global boundary nodes (cf. Remark 2.1) and T = {ti ∈ R | i ∈ I = {1, ..., n}} a set of discrete times. For a fixed time ti ∈ T , we define the history set Thi (ti ) := {t1 , ..., ti } with cardinality nhi (ti ) := |Thi (ti )| = i and the future set Tfu (ti ) := {ti+1 , ..., tn } with cardinality nfu (ti ) = |Tfu (ti )| = n − i. For simplicity, we omit the dependency of nhi , nfu , Thi and Tfu on ti if ti is clear from the context. The spatio-temporal node forecasting problem with global boundary conditions refers to the learning problem (cv (t))v∈V \Vb ,t∈Thi , (cv (t))v∈Vb ,t∈Thi ⊔Tfu ) 7−→ (cv (t))v∈V \Vb ,t∈Tfu , that is the problem of inferring the signal c at times in the future set from the signal c at times in the history set and at boundary nodes Vb . Throughout the rest of this section, we propose MeGA-MP, a MPNN that solves this problem for advectiondominated edge dynamics, where the signal c usually refers to a concentration. More specifically, in Section 3.1, we derive advection message passing by relating the edge dynamics to message generation and the coupling conditions to message aggregation and update. In Section 3.2, we present the overall MeGA-MP algorithm taking numerical aspects due to the discretization of time (Definition 3.1) into account. 3.1
Physical Derivation of MeGA-MP
Characteristic Curves We model advection of the concentration cuv : R × Ωuv → R≥0 over time R and along the edge space Ωuv = (0, luv ) in between two nodes v ∈ V and u ∈ N (v) as the transport process governed by Equation (2) with incompressible flow field νuv : R → R. In this case, by the Method of Characteristics (MoC), we can derive characteristic curves along which the function cuv remains constant over time R and space Ωuv . Lemma 3.2 (Constant concentration along characteristic curves). Let euv ∈ E hold. If the function cuv obeys Equation (2), for each t0 ∈ R and z0 ∈ R, the curve ´t γt0 ,z0 : Tt0 ,z0 −→ R≥0 , t 7−→ cuv t, z0 + t0 νuv (s) ds ´ t′ ′ is constant on Tt0 ,z0 = {t ∈ R | 0 < z0 + t0 νuv (s) ds < luv ∀t ∈ min{t0 , t}, max{t0 , t} }. An intuition on Lemma 3.2 can be found in Appendix C.1.1. Figure 1 (left) visualizes several characteristic curves for an exemplary flow field νuv : R → R along an edge euv ∈ E of length luv = 1[m]. Following Lemma 3.2, the concentration cuv is the same at all positions along such a curve. Consecutively, for the example given in Figure 1 (right), knowing the concentration cu (t − δtuv ) at the node u ∈ N (v) (z = 0) is sufficient to determine the concentration cv (t) at node v ∈ V (z = luv ). This can be observed for the heatmap in the background of the figure, which shows a high-concentration pulse traveling from one boundary at z = 0 to the other at z = luv in δtuv ≈ 5[s]. Transport Times To formulate this relationship as a message function, we need a way to compute the start and end coordinates of characteristic curves that correspond to the sender and receiver node, respectively, as well as the time shift δtuv , which we call transport time. We also use specific terminology to differentiate 5
t
v
tuv
t
euv
Edge Space [m]
Message function
u0
2
4
6 Time [s]
8
10
12 0
2
4
6 Time [s]
8
10
12
Figure 1: Characteristic curves along the edge euv ∈ E: If a curve originates at node u at time t − δtuv and terminates at node v at time t, we call it a pass-through (dark blue). If the direction is reversed, we call it inverse pass-through (light blue). If a curve starts and ends at the same node, we call it a self-loop (dark orange) if that node is u and inverse self-loop (light orange) if that node is v. The corresponding transport time δtuv > 0 (right) denotes the duration between the curve’s start and end. The shaded area (gray) in the background shows a high-concentration pulse entering the domain at node u, that is constant along characteristic curves. the following situations, which are also described visually in Figure 1 (left). To get a deeper understanding of the different kinds of transport times, details can also be found in Appendix C.1.1. Definition 3.3 (Transport times). Let euv ∈ E and t ∈ R hold. If there exists a δtuv := δtuv (t) ∈ R, such that ˆ t νuv (s) ds = luv (−luv , 0, 0 ) and t−δtuv
ˆ t′
νuv (s) ds ∈ (0, luv )
((−luv , 0), (0, luv ), (−luv , 0))
t−δtuv ′
holds for all t ∈ (min{t − δtuv , t}, max{t − δtuv , t}), we call δtuv a transport time, or, more concretely, the pass-through (inverse pass-through, self-loop, inverse self-loop) time with respect to (w.r.t.) the edge euv ∈ E and time t ∈ R. Usually, euv ∈ E and t ∈ R are clear from the context and we just say pass-through (inverse pass-through, self-loop, inverse self-loop) time, just as we omit the dependency of δtuv on t for better readability. Message Function If the flow νuv is constant over time, the transport time δtuv is also constant, and concentrations are shifted between the nodes u and v. If the flow νuv varies, so does the transport time δtuv , and instead of a clear shift, we observe more complex time warping. As this generalizes shift, we use the prior as our message function. To compute a unique positive transport time for a given time t ∈ R, we need to measure the distance traveled in an arbitrary time interval [t − δt, t]. For this purpose, we define the function ˆ t ý : R −→ R, δt 7−→ ý(δt) := νuv (s) ds, (7) t−δt
where we omit the dependency of ý on t ∈ R, euv ∈ E and the overall flow field ν for better readability. Using ý, we can compute a positive transport time as the solution to the following optimization problem (OP): n o δtuv := min δt > 0 ý(δt) ∈ {−luv , 0, luv } . (8) In Appendix C.1.1, we prove that if there exists a solution, this solution uniquely determines one of the transport times from Definition 3.3 (cf. Theorem C.12). This in turn allows to define our advectionmessage-generation functions as ( cu (t − δtuv ) if ý(δtuv ) = ±luv ϕ(t, euv , c, ν) := . (9) cv (t − δtuv ) if ý(δtuv ) = 0 6
Especially to mention, using ý again, ϕ automatically identifies whether information is carried along a(n) (inverse) pass-through or a(n) (inverse) self-loop and outputs the information from the neighbor u ∈ N (v) or the node v itself accordingly. A rigorous derivation of Equation (9) under consideration of physical priors and the intuition behind this can be found in Appendix C.1.1. Message Aggregation and Update We integrate Equation (9) into a complete message passing scheme that determines the concentration cv (t) at a node v ∈ V and time t ∈ R as X cv (t) = wuv (t) · ϕ(t, euv , c, ν) (10) u∈N (v)
with wuv (t) as defined in Equation (6) and ϕ(t, euv , c, ν) as defined in Equation (9). Theoretical Guarantees As elaborated in Section 2, advection on metric graphs is defined by multiple physical priors, the edge dynamics (Eq. (2)) and the coupling conditions (Eq. (4)) and Eq. (5)). While the former are defined on the open edge spaces Ωuv = (0, luv ) only, the latter rely on the values cuv (·, 0) and cuv (·, luv ) on the boundary ∂Ωuv = {0, luv } of Ωuv . Therefore, suitable smoothness assumptions on cuv and νuv are required in order to relate these principles. Assumption 3.4 (Smoothness assumptions). For each euv ∈ E, let cuv : R×(0, luv ) → R≥0 be differentiable and cuv : R × [0, luv ] → R≥0 and νuv : R → R continuous. Under these assumptions, in Appendix C.1, we extend the findings from Lemma 3.2 and investigate the different phenomena as introduced in Definition 3.3 more thoroughly. Based on that, we provide a rigorous theory of why our advection message-passing scheme models the described physics exactly, with the main results as follows: Theorem 3.5 (Physical alignment). Let v ∈ V and t ∈ R hold. We assume that Assumption 3.4 holds and that for each u ∈ N (v), the OP (8) has a solution. If the function cv obeys the mixing at nodes (Eq. (4)) and if for each u ∈ N (v), the function cuv obeys the advection edge dynamics (Eq. (2)) with inflow continuity condition (Eq. (5)), then Equation (10) with ϕ as defined in Equation (9) holds. 3.2
Final Algorithm: MeGA-MP
So far, we have introduced the advection message-passing operator (Eq. (10)). We will now explain how it is applied in practice. For this, two aspects remain to be addressed: First, Equation (10) is defined for continuous time t ∈ R, whereas in practice boundary values and flow fields are typically known as measurements at discrete points in time. We therefore adapt Equation (10) to handle discrete temporal information, described in paragraph Temporal Discretization. Second, the edge lengths l = (le )e∈E and the flow velocities ν = (νe )e∈E determine how far information travels over the metric graph. Under fast flow conditions, information can traverse multiple edges of the metric graph within the time interval [t1 , tn ] (cf. Definition 3.1 and Appendix D). This cannot be captured by a single application of Equation (10), as it depends only on the immediate 1-hop node neighborhood. We therefore adopt an iterative message-passing scheme, described in paragraph Model Iterations. We define the resulting iterative scheme described by Equation (11) as metric graph advection message passing (MeGA-MP). Temporal Discretization We model the concentration cv : R → R≥0 at a node v ∈ V as a discrete time series (cv (t))t∈T over equidistant discrete times t ∈ T = {ti ∈ R | i ∈ I = {1, ..., n}} = Thi ⊔Tfu with sampling interval dt = ti+1 −ti . In accordance to Definition 3.1, we aim to estimate the signal yv := (cv (t))t∈Tfu ∈ Rnfu at all non-boundary nodes v ∈ V \ Vb and times t ∈ Tfu . To accommodate known values at times t ∈ Thi and unknown values at times t ∈ Tfu , we pad the known history xv := (cv (t))t∈Thi ∈ Rnhi with zeros, yielding (0) ĉv := (xvT , 0T ) ∈ Rnhi +nfu . Here, the zero entries are placeholders for the predictions ŷv := (ĉv (t))t∈Tfu ∈ Rnfu of yv (cf. Fig. 2 and Eq. (11) below). Since the transport times δtuv > 0 are positive by definition (0) (cf. Eq. (8)), applying Equation (10) to ĉv for all v ∈ V \ Vb shifts/warps information forward in time, populating the padded entries with the solutions to the node forecasting problem. However, as both edge lengths luv and flows νuv are still continuous, time shifts t − δtuv of the message function ϕ (Eq. (9)) are 7
1. Semi-Lagrangian Backtracing 3. Message Aggregation Inflow Neighbors
2. Message Generation 4. Update
Figure 2: Initial iteration of MeGA-MP (Eq. (11)). Messages are passed from inflow neighbors u ∈ N− (v, t) to the target node v (for all v ∈ V ): First, flow velocities νuv are converted to the transport time δtuv (Eq. (8)). Second, messages are generated via time-warping (Eq. (9)). Third, results are aggregated and assigned to the target node v (Eq. (10)). The graph is directed according to the flow for visualization purposes. Grey entries correspond to padded entries and masked features according to Equation (11). also continuous and not constrained to the discrete time points T . Therefore, time warping is approximated by interpolating between the two values at times tj ∈ T and tj+1 ∈ T for which t − δtuv ∈ [tj , tj+1 ) holds. We denote the discrete version of ϕ by φ. In our experiments, we use linear interpolation as detailed in Appendix B.1, but depending on the data, different interpolation methods can be considered. (0)
Model Iterations As mentioned above, applying Equation (10) one time to ĉv for all v ∈ V \ Vb only fully solves the node forecasting problem if all evaluations of the kind cu (t − δtuv ) in Equation (9) index known concentration values. This is only the case if either u ∈ Vb is a boundary node, or if t − δtuv < ti lies temporally before the upper bound of the known history Thi . The latter corresponds to a constraint on the flow field ν. However, generalizing this to arbitrary flow fields can be achieved by iterating Equation (10) via ( cv (t) if v ∈ Vb or t ∈ Thi = {t1 , ..., ti } (0) ĉv (t) := , 0 if v ∈ V \ Vb and t ∈ Tfu = {ti+1 , ..., tn } ( 0 if v ∈ Vb or t ∈ Thi (k) (11) ĉv (t) := P for all k ≥ 1, (k−1) w (t) φ(t, e , ĉ , ν) if v ∈ V \ Vb and t ∈ Tfu uv u∈N− (v,t) uv ĉv (t) :=
∞ X
ĉ(k) v (t)
for v ∈ V \ Vb and t ∈ Tfu ,
k=1
where ĉ(k−1) ∈ R|V |×(nhi +nfu ) corresponds to the vectorized version of the function c = (cv )v∈V that contains (k−1) the values ĉv (t) for all v ∈ V and t ∈ T = Thi ∪ Tfu . To build an intuition for how these iterations relax the constraint on ν, note that in a single iteration k, a node v ∈ V \ Vb receives information at time t ∈ Tfu from its neighbors (or itself) at earlier times t − δtuv < t. This defines a spatio-temporal inflow tree of depth one, along which information is passed. With each iteration, this spatio-temporal inflow tree extends further through the graph. As a result, the node v receives information from increasingly distant nodes at progressively earlier times, eventually from either boundary nodes Vb or the known history Thi . (k) Importantly, only at these leaves, information is initially non-zero. Therefore, summing the values ĉv (t) over all k ∈ N correctly aggregates information from antisymmetrically distributed sending nodes along paths in the inflow tree. A rigorous formalization of this intuition can be found in Appendix C.2. The initial iteration (k = 0 7→ k = 1) is visualized in Figure 2. Subsequent iterations follow the same schema, but with node feature vectors that already contain information at future time steps Tfu . Theoretical Guarantees The following Theorem 3.6 ensures that the prediction ŷv = (ĉv (t))t∈Tfu defined by the last row of Equation (11) is computable in finite time for all v ∈ V \ Vb . 8
(k)
Theorem 3.6 (Convergence of model iterations). There exists a κ ∈ N, such that ĉv (t) as defined in Equation (11) is equal to zero for all k > κ, v ∈ V and t ∈ T = Thi ∪ Tfu . For a node v ∈ Vb or a time t ∈ Thi , this statement is trivial (κ(v, t) = 0). In the proof of Theorem 3.6, we (k) show that for the more interesting case of a node v ∈ V \ Vb and a time t ∈ Tfu , ĉv (t) is only non-zero for all iterations k ∈ N which correspond to the length lp of any path p ∈ P̂(v, t) in the before-mentioned inflow tree T̂ (v, t). This naturally defines κ(v, t) = maxp∈P̂(v,t) lp as the maximum path length over all paths in the inflow tree; and κ = maxv∈V maxt∈t∈Thi ⊔Tfu κ(v, t) as the maximum of the former over all nodes and time steps. In practice, we compute the inflow tree and the node- and time-dependent κ(v, t) on the fly by accumulating the transport times along paths of the inflow tree until either a boundary node v ∈ Vb is reached or the difference between the time t ∈ Tfu and the accumulated transport time lies temporally before the upper bound of the known history Thi . If this is the case, the end of a path in the inflow tree is reached. Finally, we show that MeGA-MP (Eq. (11)) reproduces the ground-truth (Eq. (10)) if the interpolation error is zero. An extended version of Theorem 3.7 with an error bound obtained through the temporal interpolation can be found in Appendix C.2. Theorem 3.7 (Interpolation error). Let v ∈ V and t ∈ T = Thi ∪ Tfu hold. Let Assumption C.21 and C.22 hold. If the temporal interpolation error is zero, i.e., if ϕ = φ holds, we obtain ĉv (t) = cv (t). Remark 3.8 (Initial value problems and efficiency). For the sake of efficiency in both computational and memory complexity, we represent the solution of advection on a metric graph as time series at nodes only (cf. Definition 3.1). The actual values along the 1-dimensional edge spaces are solved implicitly. This enables our model to scale linearly with the number of nodes and, consequently, to handle even large graphs. For this advantage, our method cannot directly receive classical initial value problems (IVPs), where dense values cuv (t, z) along the edge spaces Ωuv are known at t = 0, except when they are constant. However, initial values along edges can be transformed into time series at nodes and vice versa. This is visualized by the gray characteristic curves in Figure 1. We will demonstrate this in Appendix D.2 as an extension of Experiment 2 in Section 5. This can also be used to recover the explicit signal along the edge spaces, as used in the demo. Runtime Complexity The runtime of MeGA-MP depends on edge lengths luv , flow velocities νuv , number of history and predictions steps nhi and nfu , respectively, and those of edges |E| in the graph. We analyze runtime in two parts: Computation of transport times and the iterative message passing procedure (cf. Algorithms 1, 2). The number of iterations required to solve the dynamics depends on the number of edges that information traverses within the future horizon Tfu (cf. Appendix A). An upper bound on the number of iterations can luv −ti , where tmin = mineuv ∈E,t∈Tfu ,νuv (t)̸=0 |νuv be defined as kup = ttnmin (t)| . In line with other message-passing algorithms, MeGA-MP also applies a message function per edge that processes vectors of length n = nhi +nfu , resulting in a worst-case runtime of O(kup n|E|). The runtime for computing the transport times δtuv from the flow field νuv (cf. Algorithm 1) depends on edge lengths luv and velocities νuv . The worst case number of iterations of the inner loop is m = min(n, maxeuv ∈E,t∈Tfu ,νuv (t)̸=0 |νuvluv (t)|dt ), and since the outer loop iterates n times for every edge in E, the resulting complexity is O(|E|mn) with an upper bound of O(|E|n2 ). Therefore, the overall worst-case runtime amounts to O(|E|(kup n + n2 )). 3.3
Extended Algorithm: MeGA-MP+
We defined MeGA-MP as the algorithm introduced in Section 3.2, characterized by the message aggregation scheme (Eq. (10)), induced from the mixing at nodes (Eq. (4)), and the physics-informed advection message function ϕ (Eq. (9)) within, derived from the advection edge dynamics (Eq. (2)) with inflow continuity condition (Eq. (5)). Formulating MeGA-MP as a message-passing operator allows straight-forward adaptations to more flexible and also trainable models. This can be achieved by extending the (discrete) message-function φ (Eq. (11)/(17)). We refer to such extensions of φ as φ+ , and denote the overall extended method as MeGA-MP+ . Note that this name stands for a concept rather than a specific method, since depending on the specific formulation of φ+ , multiple variations are possible. On the one hand, added components can 9
depend on the local metric graph neighborhood N (v), just as φ does. The resulting MeGA-MP+ operator can adapt to spatial dynamics other than linear advection. On the other hand, added components that do not depend on the neighborhood yield a MeGA-MP+ operator that can adapt to additional non-spatial dynamics, while spatially still representing linear-advection. We showcase one such extension in Section 4.
4
Experiment 1: Zero-Shot Transfer Learning of Advection-Reaction Dynamics on Real-World Metric Graphs
In this section, we demonstrate the power of MeGA-MP using water quality data, a classical example of an advection-reaction metric graph. To demonstrate that MeGA-MP can be used as a standalone solver for linear advection, but also that it can be extended by a learnable component to account for additional dynamics, we consider two physical configurations: a pure linear advection system, and a coupled advectionreaction system. The objective of this experiment is to solve the node-level regression task, predicting the temporal evolution of signals at nodes, as outlined in Definition 3.1. We benchmark both MeGA-MP for advection and an example of MeGA-MP+ for advection-reaction dynamics against baseline models, and run an ablation study to show the effect of the different components of both models. Since advection is the only spatial dynamic in both configurations, MeGA-MP and MeGA-MP+ should theoretically generalize to unseen topologies without retraining. We validate this zero-shot capability by testing both configurations on an additional metric graph structure not contained in the training set. The datasets are described in detail in Section 4.1, the comparison to baselines is in Section 4.2, and the ablation study is in Section 4.3. 4.1
Domain and Dataset
Domain In this experiment, we consider water distribution systems (WDSs), in which nodes, such as reservoirs or consumer junctions, are connected through pipes. As the pipes are associated with lengths, this domain can be modeled as a metric graph. The WDSs Hanoi (Vrachimis et al., 2018) and L-Town (Vrachimis et al., 2020) are shown in Figures 3 and 4, respectively. Both constitute well-known benchmark networks, used by engineers in the water community. Pipe Junction Reservoir
Figure 3: The WDS Hanoi with 32 nodes (31 con- Figure 4: The WDS L-Town with 784 nodes (782 sumer nodes Vc and one reservoir node Vr ) and 34 consumer nodes Vc and two reservoir nodes Vr ) and 908 pipes used only for testing. pipes used for training and testing. A common task in WDSs is to estimate the water quality by predicting the distribution of a substance throughout the WDS. The two main drivers of such processes are advection and reactions (Rossman et al., 2020). Reactions describe purely temporal changes in state variables such as chemical and physical transformations, interactions with the environment, and decay. Typically, such dynamics are expressed as ODEs. When combined with advection, they extend Equation (2) as follows: ∂cuv (t, z) ∂cuv (t, z) = − νuv (t) + fr (t, z, euv ), ∂t ∂z 10
(12)
where fr is a substance-specific reaction term. For this experiment, we simulate the injection of chlorine into the WDS. Chlorine decays over time due to chemical and physical interactions. Details on the specific formalization can be found in Appendix D.1. Dataset We generate two training, validation and test datasets by simulating water quality dynamics on the Hanoi WDS under the two previously mentioned configurations. Specifically, the first dataset models pure linear advection (Eq. (2)), while the second dataset models coupled advection-reaction dynamics (Eq. (12)). To evaluate out-of-distribution topology generalization, we extend each test set by simulating the same dynamics on a different metric graph, the L-Town WDS. All ground truth simulations are generated using EPyt-Flow (Artelt et al., 2024) based on EPANET-MSX (Shang et al., 2023), a hydraulic and water quality simulator for WDSs that supports advection and user-defined reactions such as Equation (32). Details on the parametrization of these setups can be found in Appendix D.1. Each sample of the datasets consists of a time series T over n = 701 steps, where the first time step t1 corresponds to the history set Thi (t1 ) = {t1 } (nhi = 1) and the remaining nfu = 700 steps correspond to the future set Tfu (t1 ) (cf. Definition 3.1). 4.2
Experiment 1.1: Comparison to Baselines
Models We consider two variants of our method: MeGA-MP, the purely physics-informed message-passing operator derived in Section 3.2, and MeGA-MP+ , an extension with a non-spatial learnable component enabling adaptation to the advection-reaction dynamics (Eq. (12)) by learning the reaction part from data. In the spirit of Section 3.3, we modify MeGA-MP by extending its message function φ (Eq. (11)/(17)) with a multi layer perceptron (MLP). The new message function φ+ is designed to be suitable for first-order reaction kinetics, where the solution will be an exponential function. The resulting message function φ+ is given by φ+ (t, euv , ĉ(k) , ν) = φ(t, euv , ĉ(k) , ν) · exp(|δtuv | · MLP(ρ−1 uv )),, {z } | {z } | Advection
(13)
Reaction
where ρuv ∈ (0, 1] denotes normalized pipe diameters. Details on the architecture and the training of the MLP can be found in Appendix D.1. To show the ability of our model for transfer learning, we train it on the WDS Hanoi (cf. Figure 3) and test it on both Hanoi and the unseen, much larger WDS L-Town (cf. Figure 4). Baselines As this is the first approach for modeling advection on metric graphs, we select a few representative baselines from the general graph domain. We compare ourselves to the GNN by Gilmer et al. (2017) (“NNConv”) as a representative of conventional GNNs, the GNN specifically designed for WDSs by Hermes et al. (2025) (“PDE-GNN”), the transformer GNN by Rampášek et al. (2022) (“GPSConv”) as a representative for GNNs potentially able to handle global graph dependencies, and the dynamic-system based GNN by Gravina et al. (2023) (“A-DGN”) as a representative for GNNs potentially able to handle dynamical systems with long-range dependencies on graphs. For details regarding these baselines, we also refer to Appendix D.1. There, we also provide details on why the adaptation of the methods for metric graphs from Blechschmidt et al. (2022; 2025) and Laczkó et al. (2025) were not applicable to our task. Results and Analysis Table 1 reports the average mean absolute error (MAE)2 for both configurations and for both topologies together with the standard deviation over all samples, nodes and time steps. Compared to all baselines, MeGA-MP and MeGA-MP+ show superior performance in the dynamics they are designed for, while MeGA-MP can – as expected – not accurately model the additional reaction dynamics (cf. Section 4.3). Under zero-shot transfer from Hanoi to L-Town, they show no significant increase in error, whereas the performance of the baseline models degrades. An exception is the A-DGNDia model, however, this baseline underfits in the training setting, indicating that the improvement in performance may be due to factors other than actual generalization. Table 1 also reports the compute time required to run the models on individual samples (batch-size = 1). Since we use the transport times as input to each model, we report 2 We clip the predictions of all baselines to the upper and lower bound of chlorine injections as especially PDE-GNN explodes for some samples in the L-Town evaluation dataset.
11
Table 1: MAE between simulator ground truth (EPANET-MSX) and model prediction of the advection and advection-reaction dynamical systems averaged across samples, nodes and time steps of the test datasets. Hanoi is the graph seen in training, L-Town is the generalization network not seen in training. The computation time for both networks and for each method is reported below (batch size = 1). The preprocessing time is the time it takes to convert the flow field to transport times (Algorithm 1). For a fair comparison to the simulator, this time has to be added to the computation time of all ML-driven models.
Comp. Time (sec)
MAE
WDS
NNConv PDE-GNN GPSConv A-DGN A-DGNDia MeGA-MP MeGA-MP+ NNConv PDE-GNN GPSConv A-DGN A-DGNDia MeGA-MP MeGA-MP+ preprocessing Simulator
Advection Hanoi L-Town (train/test graph) (test graph)
Advection-Reaction Hanoi L-Town (train/test graph) (test graph)
|V | = 32, |E| = 34
|V | = 784, |E| = 908
|V | = 32, |E| = 34
|V | = 784, |E| = 908
0.1693±0.1610 0.0442±0.0770 0.1878±0.2911 0.0916±0.0885 0.4260±0.3247 0.0015±0.0076 0.0015±0.0076 0.0041±0.0278 1.8648±0.0168 0.0085±0.0004 0.0266±0.0004 0.0074±0.0001 0.0447±0.0055 0.0450±0.0050 0.0920±0.0295 0.6724±0.1567
0.2758±0.2345 0.1640±0.1909 0.7429±2.7174 0.2568±0.2225 0.2502±0.1612 0.0017±0.0061 0.0018±0.0061 0.0029±0.0003 1.9251±0.0162 0.0145±0.0005 0.0271±0.0004 0.0420±0.0003 0.3964±0.0189 0.4068±0.0194 1.5156±0.1350 31.6906±2.6501
0.1242±0.1132 0.0651±0.0808 0.1381±0.1137 0.1841±0.1310 0.4441±0.3340 0.1062±0.1240 0.0013±0.0037 0.0023±0.0001 1.8926±0.0231 0.0083±0.0001 0.0260±0.0004 0.0073±0.0001 0.0449±0.0049 0.0450±0.0050 0.0992±0.0299 0.7520±0.1642
0.2248±0.2113 0.4136±0.4206 0.7319±1.7254 0.2606±0.1280 0.1614±0.1060 0.2349±0.2236 0.0008±0.0018 0.0028±0.0001 1.8925±0.0170 0.0176±0.0085 0.0267±0.0001 0.0418±0.0004 0.4063±0.0202 0.4036±0.0173 1.5184±0.1270 34.8503±2.9193
the latter without the time it takes to convert the flow field to transport times and instead disclose it as a separate row. Even when adding the preprocessing time to the compute time, each of the models is faster on the large test network L-Town than the simulator that we use as the ground truth. The NNConv model is the fastest, as it is a relatively simple MPNN. The runtime complexity of MeGA-MP and MeGA-MP+ , on the other hand, additionally depend on the edge lengths and flow field, as detailed in Section 3.2. However, while MeGA-MP and MeGA-MP+ are slower than most baselines, it is faster than PDE-GNN, which has second-best performance on the in-distribution setting. We now compare the results qualitatively, first for the pure advection dataset, then for the advection-reaction dataset. Pure Advection Configuration: In Figure 5, we show the predictions of all models evaluated on the metric graph Hanoi. Note that in this case, MeGA-MP does not have to learn, while the other GNN-baselines do. The PDE-GNN models the time series the most accurate, the error mostly stems from misalignment, indicating a failure to properly model transport times. A-DGN can follow the dynamics but exhibits smoothing and fails to capture the high frequency components of the signal. NNConv and GPSConv clearly underfit the data. In contrast, MeGA-MP is able to accurately model the time series at all visualized nodes. In Figure 6, we show the predictions of all models evaluated on L-Town. It can be observed that none of the baselines can precisely capture the signal at all nodes. PDE-GNN and A-DGN can follow the dynamics on some nodes, like nodes 1, 2, 3, and 8. These are the nodes closest to the injection/boundary nodes (stars), which suggests a failure to capture the long-range spatio-temporal dependencies between more distant nodes. MeGA-MP does not show a clear deviation from the ground truth.
12
1
2
1
0 1
6
1
2
0 1
7
3
0 1
4
4
5 3
0 1
5
0 1
6
Pipe Junction Injection Nodes MeGA-MP NNConv PDE-GNN GPSConv A-DGN Simulator
0 1
8
7
0 1
10 2
8
10 1 Avg. Flow Velocity [ms ]
0
0
2
4
6 8 Time [hours]
10
Figure 5: Predictions of baselines and MeGA-MP for the advection dynamical system at a randomly chosen subset of nodes of the WDS Hanoi. Hanoi is the graph on which all models are trained. Edge color corresponds to the flow velocity averaged across all time steps of the sample. Pipe Junction Injection Nodes MeGA-MP NNConv PDE-GNN GPSConv A-DGN Simulator
5
0 1
2
0 1
3
0 1
7
4
4
3 1
8
1
1
2
0 1
5
0 1
6
0 1
6
7
0 1
10 1 10 2 Avg. Flow Velocity [ms ]
10 3
8
0
0
2
4
6 8 Time [hours]
10
Figure 6: Predictions of baselines and MeGA-MP for the advection dynamical system at a randomly chosen subset of nodes of the WDS L-Town. The L-Town metric graph is not part of the training set. Edge color corresponds to the flow velocity averaged across all time steps of the sample.
Advection-Reaction Configuration: In Figure 7, we show the predictions of all models evaluated on Hanoi. The observations are similar to the purely advective setting, the only significant difference is that A-DGN underfits here, and for most nodes resembles a highly-smoothed average signal, similar to the other two baselines NN-Conv and GPSConv. Again, MeGA-MP+ accurately models the signal. In Figure 8, 13
showing predictions of all models evaluated zero-shot (without any retraining) on L-Town, a similar degree of underfitting can be observed. Only MeGA-MP+ can model the signal with high precision. Overall, MeGA-MP and MeGA-MP+ are the only models that perform consistently accurate over multiple nodes and generalize well from the training WDS Hanoi to the test WDS L-Town on the advection and the advection-reaction dynamics, respectively. The different generalization behaviors of the different baselines emphasize the challenges in modeling advection-dominated dynamics on metric graphs and why a suitable physical prior is necessary.
Conc.
Conc.
Conc.
Conc.
MeGA-MP +
NNConv
1.0 1
PDE-GNN
GPSConv
A-DGN
6 Time [hours]
8
Simulator
2
0.5 0.0 1.0 3
4
0.5 0.0 1.0 5
6
0.5 0.0 1.0 7
8
0.5 0.0
0
2
4
6 Time [hours]
8
10
0
2
4
10
Figure 7: Predictions of baselines and MeGA-MP+ for the advection-reaction dynamical system on the WDS Hanoi. The selected nodes are the same as in Figure 5.
Conc.
Conc.
Conc.
Conc.
MeGA-MP + 1.0 1
NNConv
GPSConv
A-DGN
6 Time [hours]
8
Simulator
2
0.5 0.0 1.0 3
4
0.5 0.0 1.0 5
6
0.5 0.0 1.0 7
8
0.5 0.0
0
2
4
6 Time [hours]
8
10
0
2
4
10
Figure 8: Predictions of baselines and MeGA-MP+ for the advection-reaction dynamical system on the WDS L-Town. The selected nodes are the same as in Figure 6. In the shown data sample, PDE-GNN explodes and has therefore been excluded.
14
4.3
Experiment 1.2: Ablation Studies
In this section, we investigate how much the physics-informed components of MeGA-MP and MeGA-MP+ contribute to the superior performance observed in Section 4.2. Models In our first ablation study, we investigate the contribution of the two components in the message function φ+ of MeGA-MP+ , i.e., the non-learnable and physics-informed advection message function φ (Eq. (11)/(17)), derived from the advection edge dynamics (Eq. (2)) with inflow continuity condition (Eq. (5)), and the learnable reaction MLP (Eq. (13)). We do so by comparing the performance of MeGA-MP – designed for advection dynamics – and MeGA-MP+ – designed for advection-reaction dynamics – on both the dataset for advection and advection-reaction dynamics. In our second ablation study, we investigate the contribution of the message aggregation scheme (Eq. (10)), induced from the mixing at nodes (Eq. (4)). We do so by training a model that has the same message aggregation structure as given in Equation (10), but where we replace the message function φ not by φ+ , which yields MeGA-MP+ , but by a fully learnable MLP φ− . Specifically, we design φ− as a 2-layer MLP with SeLU activation (Klambauer et al., 2017) which has access to adjacent node information over the whole time horizon T = Thi ∪ Tfu and the corresponding pre-processed transport times: (k−1) , ĉ , δt φ− (t, euv , ĉ(k−1) , ν) = MLP ĉ(k−1) , uv v u t
(k−1)
∈ Rnhi +nfu corresponds to the vectorized where similar to ĉ(k−1) ∈ R|V |×(nhi +nfu ) (cf. Eq. (11)), ĉv (k−1) (t) for all t ∈ T = Thi ∪ Tfu , and δtuv version of the function c = (cv )v∈V containing the values ĉv corresponds to the pre-processed vector δtuv = (δtuv )t∈Thi ∪Tf u ∈ Rnhi +nfu of transport times. We denote this modification of MeGA-MP by MeGA-MP− . An overview of the different models in these ablations can be found in Table 2. The setup and datasets used in this ablation are equal to the setup of experiment 1 described in Section 4.2. Table 2: MAE between simulator ground truth (EPANET-MSX) and predictions of different modifications to MeGA-MP on the advection and advection-reaction dynamical systems averaged across samples, nodes and time steps of the test datasets. The differences between the models is highlighted by the first three columns.
MeGA-MP MeGA-MP+ MeGA-MP−
Msg. update via Eq. (10) ✓ ✓ ✓
Msg. function
Learn. comp.
φ φ+ φ−
× ✓ ✓
Advection Hanoi L-Town (train/test graph) (test graph)
Advection-Reaction Hanoi L-Town (train/test graph) (test graph)
|V | = 32, |E| = 34
|V | = 784, |E| = 908
|V | = 32, |E| = 34
|V | = 784, |E| = 908
0.0015±0.0076 0.0015±0.0076 0.0910±0.0875
0.0017±0.0061 0.0018±0.0061 0.2169±0.1722
0.1062±0.1240 0.0013±0.0037 0.0971±0.0852
0.2349±0.2236 0.0008±0.0018 0.1873±0.3045
Results and Analysis Table 2 also shows the performance of the different modifications to MeGA-MP. As an extension of Table 2, we provide corresponding figures in Appendix D.1. Firstly, it is evident that MeGA-MP− performs significantly worse than the other two architectures in the purely advective setting. This indicates that the full strength of our approach does not only come from the message aggregation scheme (Eq. (10)), induced from the mixing at nodes (Eq. (4), but also from the inductive bias in the advection message functions φ (Eq. (11)/(17)) and thus, in φ+ (Eq. (13)), derived from the advection edge dynamics (Eq. (2)) with inflow continuity condition (Eq. (5)). Removing the physical priors from φ and thus, φ+ by replacing it with a fully learnable function φ− does not suffice to model the complex antisymmetric and longrange dependencies induced by the advection dynamics. Secondly, the on-par-performance of MeGA-MP and MeGA-MP+ in the advection setting indicates that the additional MLP is not required when advection is the only dynamic. In the case of advection-reaction, only MeGA-MP+ is able to accurately capture the added complexity from the reaction term, emphasizing the utility of the simple modification in φ+ . Due to the inability to learn 15
dynamics, MeGA-MP fails to capture the reaction dynamic without the modification used in MeGA-MP+ . Furthermore, the fact that MeGA-MP− performs worse than MeGA-MP+ indicates that the other GNNbased baselines from Section 4.2 do not perform worse because of a different message aggregation function. In summary, our ablation studies show the utility of coding the advection characteristics directly into the model architecture. Without it, the model cannot reliably capture the challenging spatio-temporal dependencies. This once again emphasizes the strength of physics-informed ML and carefully constructed inductive bias over purely data-driven approaches in advection-dominated settings.
5
Experiment 2: Advection on a 1D Euclidean Domain
In this section, we empirically validate MeGA-MP by using it as a numerical solver for advection in a domain where well-known versatile numerical solvers exist: On a 1-dimensional Euclidean domain. The advantage of this domain is that it also resembles a metric path graph. The goal of this experiment is to compare the characteristics of MeGA-MP to those of classical solvers. To this end, we compare our method against two classical baselines: A semi-Lagrangian solver, as well as a method-of-lines formulation integrating over time with a Runge-Kutta (RK4) solver using a WENO5 scheme (Jiang & Shu, 1996) to approximate the spatial derivatives. For details regarding the baselines, we refer to Appendix D.2. Domain and Dataset The domain for this experiment is Ω = [0, l], describing a 1-dimensional space of length l = 100. We aim to predict the solution c(t, z) of a synthetic boundary value problem (BVP) at equidistant grid points t and z of the time and space interval [0, 100] and [0, 100] for different temporal and 100 100 100 100 100 , 316 , 100 , 32 , 10 }, respectively. Conceptually, this models an spatial discretizations dt and dz from { 1000 injection at the boundary z = 0, advecting along Ω over time t ∈ [0, 100]. To enable a fair comparison, we choose the BVP such that an analytical solution exists. The setup is described in detail in Appendix D.2. In light of our discussion on IVPs and efficiency (Remark 3.8), we discuss a similar setup posed as an IVP with boundary values in Appendix D.2, too. Results and Analysis Table 3 shows the MAE of the different solvers at different space-time resolutions. We find that MeGA-MP approximates the ground truth well and provides results competitive to other solvers. It can also be observed that MeGA-MP shows the best performance when the temporal resolution is high (dt = 0.1), while it is less affected by changes in spatial resolution. This is an expected result, as our method interpolates the signal at nodes temporally. Therefore, the interpolation error correlates with the temporal resolution, as derived in depth in Appendix C.2. This is in contrast to the baseline semi-Lagrangian solver, which interpolates spatially and thus benefits from higher spatial resolution. RK4 becomes unstable for certain combinations of spatial and temporal resolutions. The MeGA-MP solutions for dz = 1.0 are visualized in Figure 9, where the numerical smoothing from temporal interpolation is visible, as the solution degrades with increasing dt. This is amplified because of Table 3: MAE between analytical ground truth and different solver solutions over all temporal and spatial grid points. The entry “–” indicates divergence of the method. Find the full version of this table and the runtime analysis in Appendix D.2. dt 0.1 1.0 10.0
dz 0.10 1.00 10.00 0.10 1.00 10.00 0.10 1.00 10.00
semi- Lagr. 0.0016±0.0050 0.0108±0.0230 0.0562±0.0764 0.0080±0.0144 0.0097±0.0137 0.0545±0.0720 0.0763±0.0964 0.0728±0.0923 0.0671±0.0842
RK4 – 0.0035±0.0073 0.0288±0.0415 – 0.0088±0.0155 0.0291±0.0395 – – 0.0728±0.0936 16
MeGA-MP 0.0024±0.0081 0.0012±0.0068 0.0057±0.0182 0.0157±0.0298 0.0103±0.0208 0.0115±0.0222 0.0939±0.0981 0.0915±0.0952 0.0761±0.0763
dt = 0.100
dt = 0.316
dt = 1.000
dt = 3.125
50
100 0 50 100 0 Space z (dz = 1.0)
50
dt = 10.000
50 0
1.00 0.75 0.50 0.25 0.00
Solution c(z, t)
Time t
100
0
50
100 0
100 0
50
100
Figure 9: MeGA-MP solution of the BVP (33) for a fixed spatial discretization of dz = 1.0 for varying values of temporal discretizations dt. Each vertical pixel row corresponds to the time series at a node (i.e., a grid point), each horizontal row corresponds to one time step across all nodes. The sinusoidal flow field (Eq. (35)) is always positive, i.e. the injected signal at the left boundary (z = 0) is transported towards the right boundary (z = 100).
the sharp peak and therefore high-frequency content of the signal. One can also observe that smoothing increases in z-direction, i.e. in the direction of the flow field, as this coincides with the direction of spatial integration and error accumulation (cf. Appendix C.2). More detailed visualizations of all combinations of space-time discretization, an extended version of Table 3 and a runtime analysis can be found in Appendix D.2.
6
Discussion & Conclusion
In this work, we have introduced MeGA-MP, a message-passing formulation solving linear advection dynamical systems on metric graphs. We have shown theoretically and empirically that it is a valid solver with a provable error bound and that it can be extended to solve more challenging, advection-dominated dynamics. In this work, we used an MLP as the extension. However, following the structural archetype of MPNNs in principle also allows more complex integrations with existing GNNs. The comparison to other PDE solvers shows the saliency of our method, as our solution is comparable to that of established numerical solvers in scenarios where a comparison is possible. By experimenting on dynamics in WDSs we did not only show that our method can solve advection-dominated dynamics on metric graphs, but also that MeGA-MP(+ ) is directly applicable as a surrogate model for water quality, specifically the modeling of chlorine. It closely models the simulator EPANET-MSX, while being considerably faster and more flexible, as it integrates much better into the existing ML infrastructure. Successful zero-shot transfer to unseen metric graphs further strengthens the positioning of MeGA-MP(+ ) as a viable surrogate model. We attribute this transferability to novel geometries mainly to the inductive bias in our message-passing scheme, allowing a separation of spatial integration of graph signals by MeGA-MP and data-driven adaptation to the additional non-spatial reaction dynamic using MeGA-MP+ . It is unclear to what extent the transferability persists when learnable spatial components are added. We will discuss this in the following sections, but argue that in the presence of advection, MeGA-MP always contributes a suitable inductive bias at low cost. Limitations MeGA-MP(+ ) is specifically designed for linear advection-dominated dynamics on metric graphs and is naturally limited to PDEs where linear advection is the dominant spatial dynamic. The approximation error is controlled by two factors: the temporal discretization and the number of message passing iterations required to propagate information across the graph. The latter depends on the edge lengths and flow velocities; short lengths and fast flows require more iterations, leading to greater accumulation of temporal interpolation errors. Both effects are consistent with our theoretical error bound (Theorem 3.7) and empirically supported by Experiment 2 (Section 5). We therefore recommend choosing sufficiently fine temporal discretizations relative to the flow velocities in the system. Beyond numerical accuracy, MeGA-MP(+ ) currently assumes exact knowledge of the flow field, which could be unavailable in some practical settings. 17
To this end, MeGA-MP(+ ) does not provide explicit uncertainty quantification, but the low execution time allows feasible Monte-Carlo estimates of uncertainty. Future Work Several directions can build naturally on this work. First, extending MeGA-MP to more complex spatial dynamics, such as semi- and non-linear advection or advection-diffusion, and to more complex temporal dynamics, such as multi-species reactions, via additional learnable operators is a promising avenue. The latter, in particular, would bring this method closer to a full ML-driven surrogate model for the simulator EPANET-MSX for WDSs. Second, exploring inverse problems such as identifying flow fields or reaction coefficients from observational data would not only address the current limitation of requiring exact flow field knowledge, but would further emphasize the practical advantages of such surrogates over numerical simulators. Thirdly, a systematic investigation of robustness to noisy or incomplete parameters such as flow fields is an important prerequisite for deployment of such systems in the real-world. Finally, metric graphs governed by other PDEs beyond advection-dominated edge dynamics pose fundamentally different challenges. A systematic study of these different problem classes, including benchmark datasets and baseline evaluations, is an important direction for future work and the comparability of different methods across domains.
Impact Statement This work presents a physics-informed ML method for modeling spatio-temporal dynamics on advectiondominated metric graph domains. The data used in this work is self-generated and does not contain personal or sensitive information, ensuring no privacy concerns under GDPR or comparable regulations, or ethical concerns. However, real-world deployment of such methods may require sensor data or consumption records that can be linked to individual consumers. In such cases, data collection and processing must comply with applicable data protection regulations, and anonymization or aggregation techniques should be applied where possible. The primary intended application of this work is modeling in scientific and engineering context, with the goal of improving the efficiency and accuracy of simulations. We also acknowledge the potential broader impacts of physics-informed ML including misuse and emphasize that responsible development and independent validation are essential. For example, in safety-critical settings, incorrect model predictions could have negative consequences for the provision of essential services. Therefore, before deploying such systems in the real world, a thorough investigation of failure modes is essential, including error-proneness under distribution shift, sensitivity to misspecified physical parameters and the absence of uncertainty quantification in the current formulation.
LLM Statement This research was partially assisted by using large language models (LLMs) such as ChatGPT, Gemini or Claude. More precisely, LLMs were used for retrieval and discovery, such as finding related work. However, all related work cited in this paper was read and verified by the authors. Moreover, LLMs were used to generate boilerplate code that has been checked for correctness by the authors. Finally, LLMs were used to assist with editing and proofreading of selected parts of the manuscript with a focus on improving clarity and readability. The scientific content, the experimental design, the code implementing the novel methods presented in this work, the results, most parts of the text in this manuscript as well as the theory and equations were not generated by LLMs.
References Takuya Akiba, Shotaro Sano, Toshihiko Yanase, Takeru Ohta, and Masanori Koyama. Optuna: A nextgeneration hyperparameter optimization framework. In Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, 2019. Mario Arioli and Michele Benzi. A Finite Element Method for Quantum Graphs. IMA Journal of Numerical Analysis, 38(3):1119–1163, 2018. 18
André Artelt, Marios S. Kyriakou, Stelios G. Vrachimis, Demetrios G. Eliades, Barbara Hammer, and Marios M. Polycarpou. Epyt-flow: A toolkit for generating water distribution network data. Journal of Open Source Software, 9(103):7104, 2024. Gregory Berkolaiko. An Elementary Introduction to Quantum Graphs. Geometric and Computational Spectral Theory, 700:41–72, 2017. Gregory Berkolaiko and Peter Kuchment. Introduction to Quantum Graphs. Number 186. American Mathematical Soc., 2013. Jan Blechschmidt, Jan-Frederik Pietschman, Tom-Christian Riemer, Martin Stoll, and Max Winkler. A comparison of PINN approaches for drift-diffusion equations on metric graphs. arXiv preprint arXiv:2205.07195, 2022. Jan Blechschmidt, Tom-Christian Riemer, Max Winkler, Martin Stoll, and Jan-F Pietschmann. PhysicsInformed DeepONets for drift-diffusion on metric graphs: simulation and parameter identification. In Proceedings of the 42th International Conference on Machine Learning (ICML), volume 267 of Proceedings of Machine Learning Research, pp. 4482–4506. PMLR, 2025. Lucas Böttcher and Mason A Porter. arXiv:2401.00735, 2024.
Dynamical Processes on Metric Networks.
arXiv preprint
Johannes Brandstetter, Daniel Worrall, and Max Welling. Message Passing Neural PDE Solvers. In Proceedings of the 10th International Conference on Learning Representations (ICLR), 2022. Alberto Bressan, Sunčica Čanić, Mauro Garavello, Michael Herty, and Benedetto Piccoli. Flows on Networks: Recent Results and Perspectives. EMS Surveys in Mathematical Sciences, 1(1):47–111, 2014. Xavier Bresson and Thomas Laurent. Residual gated graph convnets. CoRR, abs/1711.07553, 2017. URL http://arxiv.org/abs/1711.07553. Moysey Brio, J-G Caputo, and Hannah Kravitz. Spectral Solutions of PDEs on Networks. Applied Numerical Mathematics, 172:99–117, 2022. Ben Chamberlain, James Rowbottom, Maria I Gorinova, Michael Bronstein, Stefan Webb, and Emanuele Rossi. GRAND: Graph Neural Diffusion. In Proceedings of the 38th International Conference on Machine Learning (ICML), volume 139 of Proceedings of Machine Learning Research, pp. 1407–1418. PMLR, 2021. Airlie Chapman. Advection on Graphs. Semi-Autonomous Networks: Effective Control of Networked Systems through Protocols, Design, and Modeling, pp. 3–16, 2015. Jeongwhan Choi, Seoyoung Hong, Noseong Park, and Sung-Bae Cho. GREAD: Graph Neural ReactionDiffusion Networks. In Proceedings of the 40th International Conference on Machine Learning (ICML), volume 202 of Proceedings of Machine Learning Research, pp. 5722–5747. PMLR, 2023. Moshe Eliasof, Eldad Haber, and Eran Treister. PDE-GCN: Novel Architectures for Graph Neural Networks Motivated by Partial Differential Equations. In Proceedings of the 35th International Conference on Neural Information Processing Systems (NeurIPS), 2021. Moshe Eliasof, Eldad Haber, and Eran Treister. Feature Transportation Improves Graph Neural Networks. In Proceedings of the 38th AAAI Conference on Artificial Intelligence, number 11, pp. 11874–11882, 2024. Justin Gilmer, Samuel S Schoenholz, Patrick F Riley, Oriol Vinyals, and George E Dahl. Neural Message Massing for Quantum Chemistry. In Proceedings of the 34th International Conference on Machine Learning (ICML), volume 70 of Proceedings of Machine Learning Research, pp. 1263–1272. PMLR, 2017. Alessio Gravina, Davide Bacciu, and Claudio Gallicchio. Anti-Symmetric DGN: A Stable Architecture for Deep Graph Networks. In Proceedings of the 11th International Conference on Learning Representations (ICLR), 2023. 19
Martin Gugat and Michael Herty. Modeling, Control, and Numerics of Gas Networks, volume 23, pp. 59–86. Elsevier, 2022. Martin Gugat, Michael Herty, Axel Klar, and Günter Leugering. Optimal Control for Traffic Flow Networks. Journal of Optimization Theory and Applications, 126(3):589–616, 2005. Richard Habermann. Applied Partial Differential Equations with Fourier Series and Boundary Value Problems, volume 5. Pearson, 2013. Luca Hermes, André Artelt, Stelios G Vrachimis, Marios M Polycarpou, and Barbara Hammer. A Benchmark for Physics-informed Machine Learning of Chlorine Concentration States in Water Distribution Networks. SN Computer Science, 6(5):522, 2025. Weihua Hu*, Bowen Liu*, Joseph Gomes, Marinka Zitnik, Percy Liang, Vijay Pande, and Jure Leskovec. Strategies for pre-training graph neural networks. In International Conference on Learning Representations, 2020. URL https://openreview.net/forum?id=HJlWWJSFDH. Guang-Shan Jiang and Chi-Wang Shu. Efficient implementation of weighted eno schemes. Journal of Computational Physics, 126(1):202–228, 1996. ISSN 0021-9991. doi: https://doi.org/10.1006/jcph.1996. 0130. URL https://www.sciencedirect.com/science/article/pii/S0021999196901308. Diederik P. Kingma and Jimmy Ba. Adam: A method for stochastic optimization. CoRR, abs/1412.6980, 2014. URL https://api.semanticscholar.org/CorpusID:6628106. Günter Klambauer, Thomas Unterthiner, Andreas Mayr, and Sepp Hochreiter. Self-normalizing neural networks. In I. Guyon, U. Von Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett (eds.), Advances in Neural Information Processing Systems, volume 30. Curran Associates, Inc., 2017. URL https://proceedings.neurips.cc/paper_files/paper/2017/file/ 5d44ee6f2c3f71b73125876103c8f6c4-Paper.pdf. Nikola Kovachki, Zongyi Li, Burigede Liu, Kamyar Azizzadenesheli, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Neural Operator: Learning Maps Between Function Spaces With Applications to PDEs. The Journal of Machine Learning Research, 24(1):1–97, 2023. Peter Kuchment. Quantum graphs: I. Some Basic Structures. Waves in Random media, 14(1):S107, 2003. Peter Kuchment. Quantum Graphs: II. Some Spectral Properties of Quantum and Combinatorial Graphs. Journal of Physics A: Mathematical and General, 38(22):4887–4900, 2005. Peter Kuchment. Quantum Graphs: An Introduction and a Brief Survey. arXiv preprint arXiv:0802.3442, 2008. Csongor L Laczkó, Mihály A Vághy, and Mihály Kovács. A transferable PINN-based method for quantum graphs with unseen structure. IFAC-PapersOnLine, 59(1):67–72, 2025. Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Fourier Neural Operator for Parametric Partial Differential Equations. In Proceedings of the 9th International Conference on Learning Representations (ICLR), 2020a. Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Neural Operator: Graph Kernel Network for Partial Differential Equations. arXiv preprint arXiv:2003.03485, 2020b. Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Andrew Stuart, Kaushik Bhattacharya, and Anima Anandkumar. Multipole Graph Neural Operator for Parametric Partial Differential Equations. In Proceedings of the 34th International Conference on Neural Information Processing Systems (NeurIPS), 2020c. 20
Zongyi Li, Hongkai Zheng, Nikola Kovachki, David Jin, Haoxuan Chen, Burigede Liu, Kamyar Azizzadenesheli, and Anima Anandkumar. Physics-Informed Neural Operator for Learning Partial Differential Equations. ACM/IMS Journal of Data Science, 1(3):1–27, 2024. Lu Lu, Pengzhan Jin, and George Em Karniadakis. Deeponet: Learning nonlinear operators for identifying differential equations based on the universal approximation theorem of operators. arXiv preprint arXiv:1910.03193, 2019. Lu Lu, Pengzhan Jin, Guofei Pang, Zhongqiang Zhang, and George Em Karniadakis. Learning nonlinear operators via DeepONet based on the universal approximation theorem of operators. Nature Machine Intelligence, 3(3):218–229, 2021. Maziar Raissi, Paris Perdikaris, and George E Karniadakis. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics, 378:686–707, 2019. Ladislav Rampášek, Mikhail Galkin, Vijay Dwivedi, Anh Luu, Guy Wolf, and Dominique Beaini. Recipe for a General, Powerful, Scalable Graph Transformer. In Proceedings of the 36th International Conference on Neural Information Processing Systems (NeurIPS), 2022. Lewis A Rossman, H Woo, M Tryby, F Shang, R Janke, and T Haxton. EPANET 2.2 Users Manual. United States Environmental Protection Agency (EPA)., 2020. T Konstantin Rusch, Ben Chamberlain, James Rowbottom, Siddhartha Mishra, and Michael Bronstein. Graph-Coupled Oscillator Networks. In Proceedings of the 39th International Conference on Machine Learning (ICML), volume 162 of Proceedings of Machine Learning Research, pp. 18888–18909. PMLR, 2022. Feng Shang, Lewis A Rossman, and James G Uber. EPANET-MSX 2.0 User Manual. United States Environmental Protection Agency (EPA), 2023. Abishek Thangamuthu, Gunjan Kumar, Suresh Bishnoi, Ravinder Bhattoo, NM Krishnan, and Sayan Ranu. Unravelling the Performance of Physics-informed Graph Neural Networks for Dynamical Systems. In Proceedings of the 36th International Conference on Neural Information Processing Systems (NeurIPS), 2022. Stelios G. Vrachimis, Marios S. Kyriakou, Demetrios G. Eliades, and Marios M. Polycarpou. LeakDB: A Benchmark Dataset for Leakage Diagnosis in Water Distribution Networks. In Proceedings of the 1st International Joint Conference on Water Distribution Systems Analysis & Computing and Control for the Water Industry (WDSA/CCWI), 2018. Stelios G. Vrachimis, Demetrios G. Eliades, Riccardo Taormina, Avi Ostfeld, Zoran Kapelan, Shuming Liu, Marios Kyriakou, Pavlos Pavlou, Mengning Qiu, and Marios M. Polycarpou. BattLeDIM: Battle of the Leakage Detection and Isolation Methods. In Proceedings of the 2nd International Joint Conference on Water Distribution Systems Analysis & Computing and Control for the Water Industry (WDSA/CCWI), pp. 1–6, 2020. Shuhei Watanabe. Tree-structured parzen estimator: Understanding its algorithm components and their roles for better empirical performance, 2025. URL https://arxiv.org/abs/2304.11127.
21
Organization of the Supplementary Material A Details on the Background
23
B Details on the Method
25
B.1 Final Algorithm: MeGA-MP . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C Details on Theoretical Guarantees
25 27
C.1 Physical Derivation of MeGA-MP . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
27
C.1.1 Message Generation: Advective Transport Along Edges . . . . . . . . . . . . . . . . .
28
C.1.2 Message Aggregation and Update: Mixing at Nodes . . . . . . . . . . . . . . . . . . .
34
C.2 Error Bounds on MeGA-MP . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
35
C.2.1 Preliminaries: Trees with Children-Aggregation Scheme . . . . . . . . . . . . . . . . .
35
C.2.2 Construction of Inflow Trees
36
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
D Details on the Experiments
39
D.1 Experiment 1: Zero-Shot Transfer Learning of Advection-Reaction Dynamics on Real-World Metric Graphs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
39
D.2 Experiment 2: Advection on a 1D Euclidean Domain . . . . . . . . . . . . . . . . . . . . . . .
44
E Proofs
48
E.1 Proof of Lemma 3.2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
48
E.2 Proof of Lemma C.3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
49
E.3 Proof of Theorem C.4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
51
E.4 Proof of Lemma C.6 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
52
E.5 Proof of Lemma C.7 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
54
E.6 Proof of Lemma C.8 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
56
E.7 Proof of Lemma C.9 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
57
E.8 Proof of Theorem C.10 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
57
E.9 Proof of Theorem C.11 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
60
E.10 Proof of Theorem C.12 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
63
E.11 Proof of Lemma C.14
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
66
E.12 Proof of Theorem 3.5 and C.15 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
67
E.13 Proof of Theorem C.16 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
69
E.14 Proof of Theorem C.19 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
71
E.15 Proof of Theorem C.20 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
72
E.16 Proof of Theorem 3.6 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
74
E.17 Proof of Theorem 3.7 and Theorem C.24 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
75
22
A
Details on the Background
Advection on Metric Graphs Generally speaking, advection describes transport through the presence of a flow field νuv : R × Ωuv → R. For each time t ∈ R, the flow νuv (t, ·) : Ωuv → R defines a one-dimensional (t,·) vector field along the edge space Ωuv = (0, luv ). By definition of incompressibility, div(νuv (t, ·)) = ∂νuv = ∂z 0 holds, or equivalently, the vector field νuv (t, ·) is constant on Ωuv = (0, luv ), modeling the edge euv ∈ E. Analogously, the vector field νvu (t, ·) is constant on Ωvu = (0, lvu ) = (0, luv ), modeling the contrasting edge evu ∈ E. Since both constant vector fields νuv (t, ·) and νvu (t, ·) model the same physical phenomenon, namely the constant, actually observable physical flow field between the nodes v ∈ V and u ∈ N (v), but along contrasting directions (the beginning of the edge euv is the end of the edge evu and vice versa), the orientation of νvu (t, ·) must be the opposite of the one of νuv (t, ·). Since they both correspond to the same, constant physical flow field, νvu (t, ·) = −νuv (t, ·) (Eq. (3)) must hold, and since they do not depend on space, we omit the spatial dependency. Together with the convention of signs (cf. Section 2), this allows to formally define the inflow neighbors from v ∈ V at time t ∈ R, i.e., neighbors where there is a flow from u ∈ N (v) into v ∈ V at time t ∈ R, as N− (v, t) := {u ∈ N (v) | sgn(νvu (t)) = −1} = {u ∈ N (v) | sgn(νuv (t)) = 1}.
(14)
While incompressibility of a single flow field νuv (t, ·) implies the independence of the spatial variable z ∈ Ωuv , the incompressibility of the overall flow field ν = (νe )e∈E translates to the fact that the sum of flow rates quv (t) = auv νuv (t) into a node v ∈ V is equal to the sum of the flow rates qvu (t) = avu νvu (t) out of v: X X quv (t) = qvu (t). u∈N− (v,t)
u∈N+ (v,t)
The Challenges of Advection-dominated Dynamics on Metric Graphs Advection-dominated dynamical systems on metric graphs entail the challenge of complex, antisymmetric, and spatio-temporal dependency structures. We want visualize this using the simple example of an advection metric path graph, that is, a sequence of nodes (v0 , v1 , ..., vn ) with edges ej(j+1) in-between two consecutive nodes vj and vj+1 . We assume that the corresponding edge spaces Ωj(j+1) are of equal size lj(j+1) = 1m and have a constant flow field of νj(j+1) = 2m/s for all ej(j+1) ∈ E. If we access information at each of the nodes vj at discrete time points t0 , ..., ti , ti+1 , ... with temporal discretization dt = ti+1 − ti = 60s, information travels dt·νj(j+1) = 120 hops in the path dt · νj(j+1) = 60s · 2m/s = 120m per time step dt, which corresponds to lj(j+1) graph. By this, the signal cv (ti ) at a node v at time ti depends on the signal cu (tj ) at time tj and a node u which lies 120 hops upstream of v. If such path does not exist, the path will be shorter and will originate at boundary nodes. The above example provides a sketch how the spatial dependency structure depends on the lengths, the flow field and the temporal discretization. For more complex structured graphs and more complex configurations of PDE-parameters, the spatial dependency structure will be a tree that extends from the target (root) node v at time t to upstream nodes (leafs) as the node u in the example above (also see Appendix C.2). The dependency tree will extend further in directions with high flow rate and less in direction where the flow rate is lower. This results in two major challenges for standard MPNNs: Firstly, the extent of the receptive field can vary drastically across the space of the metric graph and over time due to varying PDE-parameters, making it impossible to find a single suitable number of message passing layers. Secondly, the receptive field can be highly asymmetric. These two properties require highly selective and anisotropic message passing that depends on the structure of the metric graph, flow field and the temporal discretization. While there already exist work that discusses graphs and GNNs specifically taylored to advection problems (Chapman, 2015; Eliasof et al., 2024), the underlying graph structure corresponds to a discretization of the Euclidean space; a paradigm different from the one on metric graphs (cf. Section 2). We highlight this difference in Figure 10.
23
The continuous advection operator:
The discrete advection operator:
Figure 10: Advection on a metric graph (top) vs. on a non-metric graph (bottom): A metric graph is a system of one-dimensional subspaces Ωuv = (0, luv ), corresponding to edges euv ∈ E, which are coupled at nodes V . Each edge euv is associated with a one-dimensional PDE that defines the dynamics of the signal cuv (t, z) over time t ∈ R and along the edge space z ∈ Ωuv = (0, luv ), parametrized by a known function νuv . In case of advection, the known functions ν = (νuv )euv ∈E correspond to the time-dependent flow field of the metric graph affecting the spatio-temporal evolution of the signals cuv (dark blue, light blue, dark orange and light orange lines). The coupling condition that defines the signal cv (t) at a node v ∈ V and over time t ∈ R can be expressed as a weighted aggregation of features ϕ(t, euv , c, ν) = cuv (t, luv ) at the boundary z = luv (orange-shaded ball around v) of the edge spaces Ωuv associated with incoming edges euv (cf. Eq. (6)). Leveraging the edge dynamics, their value can be traced back to either features cu (t − δtuv ) of neighboring nodes u ∈ N− (v, t) (dark blue lines) or features cv (t − δtuv ) of the node itself (light orange lines) and at earlier, over the time and neighboring nodes varying times t − δtuv (cf. Eq. (9)). In contrast, a (conventional, non-metric) graph can be considered as a discretization of the spatial Euclidean space Ω = Rn (gray area). Here, edges euv ∈ E preserve proximity in between nodes V , corresponding to coordinates in the Euclidean domain. The overall space Ω = Rn is associated with a multi-dimensional PDE that defines the dynamics of the signal c(t, v) := cv (t) over time t ∈ R and space v ∈ V ⊂ Ω, parametrized by a function ν. In case of advection, the discretized advection-operator of the underlying multi-dimensional v (t) and continuous advection-PDE ∂c∂t = −ν(t) · ∇cv (t, z) defines a node update cv (t + dt) that aggregates the net inflow from neighboring nodes u ∈ Nin (v) minus the outflow to neighboring nodes u ∈ Nout (v) and at a fixed time t ∈ R (black and gray lines and orange-shaded area aorund v). In practice, the discretized version ν = (νuv )euv ∈E of the continuous flow field is a learned parameter of a GNN.
24
B
Details on the Method
Message Passing has become a common principle in GNNs (Gilmer et al., 2017). It describes how local information at the nodes of a graph is passed to neighboring nodes. Message passing as defined in Gilmer et al. (2017) consists of three consecutive steps: Message generation, message aggregation, and node update. Message generation describes how a feature vector at one node is transformed as it is sent to a neighbor. Message aggregation and node update describe how multiple messages arriving at a node are accumulated into a single feature vector and update the features at the target node. A general formulation of this method is given by X xv (t + 1) = θ ϕ xu (t), xv (t), quv (t) , xv (t) , (15) u∈N (v)
P where ϕ is the message generation function, is a permutation-invariant aggregation function (e.g. sum or mean), θ is a function that updates the features xv of node v ∈ V , and quv are optional edge features of the edge euv ∈ E from node u ∈ N (v) to node v ∈ V . Moreover, t denotes the message-passing iteration step. B.1
Final Algorithm: MeGA-MP
Interpolation Different interpolation schemes can be applied for time warping, depending on the use case. A version of linear interpolation of a value cu (t − δtuv ) can be modeled as cu (t − δtuv ) ≈ a(tj ) · cu (tj ) + a(tj+1 ) · cu (tj+1 ) with a(tj+1 ) :=
t−δtuv −tj tj+1 −tj
Model Iterations
and a(tj ) := 1 −
t−δtuv −tj tj+1 −tj
(16)
for the j ∈ I such that t − δtuv ∈ [tj , tj+1 ) holds.
Based on Equation (16), we formally define φ(t, euv , ĉ(k) , ν) from Equation (11) as
( (k)
φ(t, euv , ĉ
, ν) :=
(k)
(k)
a(tj ) · ĉu (tj ) + a(tj+1 ) · ĉu (tj+1 ) (k) (k) a(tj ) · ĉv (tj ) + a(tj+1 ) · ĉv (tj+1 )
if ý(δtuv ) ∈ {−luv , luv } if ý(δtuv ) = 0
(17)
for the j ∈ I for which t − δtuv ∈ [tj , tj+1 ) holds. In practice, we precompute the transport times δtuv (k) and implement φ as an interpolation operator I[cw ](t − δtuv ) at the correct node w ∈ {u, v} that does not depend on the flow field ν anymore (cf. Algorithm 2). Further details on the implementation can be found in algorithms 1 and 2.
25
Algorithm 1 Preprocessing: Converting flows νuv (t) to transport times δtuv (t) (Definition C.5). This algorithm is implemented in a parallelized/vectorized form for efficiency. 1: Inputs: Flow velocities νuv ∈ Rn for edge euv ∈ E for n time steps (n = nhi + nfu ); Edge length luv ∈ R
for edge euv ∈ E.
2: C ← cumsum(νuv ) ∈ Rn+1 3: ∆T ∈ Rn init. empty
▷ cumulative distance traveled over time, with C0 = 0 ▷ array to store transport times 4: M sl ∈ {0, 1}n init. zero ▷ array to store if self-loop occurs at time t 5: ▷ Compute (inverse) pass-through transport times 6: for t = 0 to n − 1 do ▷ Trace characteristic curve starting at time t 7: for x = t down to 0 do ▷ Trace backward along time steps 8: dx,t ← Ct − Cx ▷ distance traveled between x and t 9: if dx,t ≥ luv then ▷ true for pass-throughs 10: interpolate x̂ between Cx and Cx+1 if needed 11: ∆Tt ← t − x̂ ; break 12: else if dx,t ≤ −luv then ▷ true for inverse pass-throughs 13: interpolate x̂ between Cx and Cx+1 if needed 14: ∆Tt ← x̂ − t ; break 15: end if 16: end for 17: end for ▷ Compute self-loop transport times 18: for t = 1 to n − 1 do ▷ Trace characteristic curve starting at time t 19: for x = t down to 1 do 20: d1 ← Cx−1 − Ct ; d2 ← Cx − Ct ▷ distance differences between consecutive steps 21: if sign(d1 ) ̸= sign(d2 ) then ▷ zero-crossing indicates self-loop 22: interpolate zero crossing: x̂ ← x − d2 /(d2 − d1 ) 23: ∆Tt ← sign(d2 ) · (t − x̂) 24: Mtsl ← 1 ; break 25: end if 26: end for 27: end for 28: Output: ∆T ∈ Rn ; M sl ∈ {0, 1}n ▷ ∆Tt > 0 ⇔ transport in edge dir., ∆Tt < 0 ⇔ transport opposing edge dir. Algorithm 2 MeGA-MP: Metric Graph Advection Message Passing 1: Inputs: Graph structure G = (V, E); Initial node states / known history x ∈ R|V |×khi ; Transport times
δtuv (t) for each edge euv ∈ E and for each time step ti ∈ T , where |T | = n = nhi + nfu . ▷ pad initial condition by prediction length nfu 3: Apply Dirichlet boundary conditions to x(0) at boundary nodes 4: xout ← x(0) ▷ initialize output accumulator 5: for i = 1 to κ do ▷ Iterate until convergence, c.f. Theorem 3.6 6: for each node v ∈ V and time t do ▷ Message-Passing P (i) (i−1) xv,t ← u∈N (v) wuv (t)(I[xu ](t − δtuv (t))1Muv sl (t)=0 7: ▷ cf. Eq. 6 and Sec. 3.2 and B.1 (i−1) + I[xv ](t − δtuv (t))1Muv sl (t)=1 ) 8: end for (i) 9: xv,t ← 0 ∀v ∈ Vb , ∀t ∈ T ▷ Set to zero at boundary nodes. (i) 10: xout ← xout + x 11: end for 12: if boundary values exist then 13: Set xout at boundary nodes 14: end if 15: Output: xout ▷ final node states after advection message passing 2: x(0) ← pad(x, nfu )
26
C
Details on Theoretical Guarantees
C.1
Physical Derivation of MeGA-MP
While in Section 3.1, we provide an intuition on how MeGA-MP incorporates the physical priors from Section 2, in this section, we present a rigorous derivation of our metric graph advection message passing algorithm based on these priors. More precisely, we prove that Equation (10) is equal to the mixing at nodes (Eq. (4)) given the advection dynamics (Eq. (2)) with inflow continuity condition (Eq. (5)) and suitable smoothness assumptions (Assumption (3.4)). We repeat these priors for convenience while giving an overview of the rest of this section. In Section C.1.1, we derive the message-generation function ϕ as defined in Equation (9) based on the advection edge dynamics. For each euv ∈ E, this refers to the change of the concentration cuv : R×Ωuv → R≥0 over time R and along the edge space Ωuv = (0, luv ), determined by the incompressible flow field νuv : R → R over time R. Formally, this is defined by the PDE ∂cuv (t, z) ∂cuv (t, z) = −νuv (t) ∂t ∂z
(Equation (2)).
In Section C.1.2, we derive the message aggregation and node update scheme cv (t) as defined in Equation (10) based on the coupling conditions given by the mixing at nodes and the inflow continuity condition. For each v ∈ V , the former refers to the concentration cv : R → R≥0 over time R, determined by a weighted aggregation of the concentrations cuv (·, luv ) carried by the flow rates quv = αuv νuv : R → R from (temporal) inflow neighbors u ∈ N− (v, ·) over time R. Formally, this is defined as P u∈N− (v,t) quv (t) · cuv (t, luv ) P cv (t) = (Equation (4)). u∈N− (v,t) quv (t) For each edge euv ∈ E, the latter refers to the concentration cuv : R × Ωuv → R≥0 over time R and at the temporal inflow boundary z = 0 of the edge space Ωuv = (0, luv ), determined by the concentration cu : R → R of the temporal inflow neighbor u ∈ N− (v, ·) over time R. Formally, this is defined by cuv (t, 0) = cu (t) for all u ∈ N− (v, t)
(Equation (5)).
The mixing at nodes and the inflow continuity condition formalize the following intuition: On the one hand, if there is a (physical) flow from u ∈ N− (v, t) to v ∈ V entering the edge euv ∈ E at time t ∈ R, the concentration cuv (t, 0) at the beginning z = 0 of that edge euv is equal to the concentration cu (t) at the node u out of which the flow flows into that edge (local inflow boundary condition). On the other hand, if there is a (physical) flow from u ∈ N− (v, t) to v ∈ V leaving the edge euv ∈ E at time t ∈ R, the concentration cuv (t, luv ) at the end z = luv of that edge euv contributes to the weighted aggregation at the node v into which the flow flows out of that edge (local outflow boundary condition). While the edge dynamics (Eq. (1)) are defined on the open edge spaces Ωuv = (0, luv ), the coupling conditions (Eq. (4) and Eq. (5)) rely on the values cuv (·, 0) and cuv (·, luv ) on the boundary ∂Ωuv = {0, luv } of Ωuv . Therefore, suitable smoothness assumptions on cuv and νuv are required to be able to transfer information in-between the edge spaces and its coupling boundaries.3 Assumption (Assumption 3.4). For each euv ∈ E, let • cuv : R × (0, luv ) → R≥0 be differentiable, • cuv : R × [0, luv ] → R≥0 be continuous and • νuv : R → R be continuous. 3 Note that strictly speaking, we also need some part of the assumptions simply for Equation (2) or more generally, Equation
(1), to be well-defined.
27
Remark C.1 (Symmetry of signals). Note that next to the physical priors introduced above, simply because of the fact that we model a metric graph as symmetric directed graph G (cf. Section 2), cvu (t, z) = cuv (t, lvu − z)
(18)
holds for all v ∈ V , u ∈ N (v), t ∈ R and z ∈ Ωvu = [0, lvu ]. That is, given that we consider a symmetric directed graph G = (V, E) and the contrasting edges evu , euv ∈ E of length lvu = luv , we expect that the signal cvu (t, z) along the edge evu at time t and position z ∈ [0, lvu ] is equal to the signal along the edge euv at time t and position lvu − z = luv − z ∈ [0, luv ] = Ωuv . C.1.1
Message Generation: Advective Transport Along Edges
Throughout the rest of this section, we assume that for each euv ∈ E (that is, for each v ∈ V and u ∈ N (v)), Equation (2) with coupling conditions as defined by Equation (4) and (5) and smothness assumptions as defined by Assumption 3.4 holds. According to Equation (6), the message-generation function ϕ(t, euv , c, ν) = cuv (t, luv ) relies on the signal cuv (·, luv ) over time R and at the temporal outflow boundary z = luv of the edge space Ωuv = (0, luv ). Since we are only aware of the signal cuv (·, 0) over time R and at the temporal inflow boundary z = 0 of each edge space Ωuv = (0, luv ) through the inflow boundary condition (Eq. (5)), computing ϕ(t, euv , c, ν) would normally require to explicitly solve the edge dynamics along the whole edge space Ωuv = [0, luv ]. This is typically done by computationally expensive numerical solvers. However, our message-generation function ϕ builds upon the MoC to avoid the explicit computation of concentrations cuv (·, z) in between the boundaries, i.e., for all z ∈ (0, luv ). More specifically, we can derive characteristic curves along which the function cuv remains constant over time R and space Ωuv . Lemma (Constant concentration along characteristic curves (Lemma 3.2)). Let euv ∈ E hold. If the function cuv obeys Equation (2), for each t0 ∈ R and z0 ∈ R, the curve ´t γt0 ,z0 : Tt0 ,z0 −→ R≥0 , t 7−→ cuv t, z0 + t0 νuv (s) ds n o ´ t′ ′ is constant on Tt0 ,z0 = t ∈ R | 0 < z0 + t0 νuv (s) ds < luv ∀t ∈ min{t0 , t}, max{t0 , t} . Intuitively, in Lemma 3.2, we fix an arbitrary time t0 ∈ R and´a position z0 ∈ R and follow any particle that t is at this position z0 at this time t0 along the path z(t) := z0 + t0 νuv (s) ds it travels during the time t−t0 (or t0 − t) for some time t ∈ Tt0 ,z0 with possibly varying flow velocities νuv (s) for all s ∈ [min{t0 , t}, max{t0 , t}]. If we assign a concentration cuv (t0 , z0 ) to this particle, the concentration should stay constant along the path this particle travels. The choice of t ∈ Tt0 ,z0 guarantees that the before-mentioned path remains within the edge euv , or in other words, in the open sub-domain R × Ωuv of the function cuv on which Equation (2) is defined on, as we will investigate in detail in the next Remark C.2. Remark C.2 (On Tt0 ,z0 ). Since for euv ∈ E, the velocity νuv : R → R can change its sign, for each t0 ∈ R and z0 ∈ R, the function ˆ t z : R −→ R, t 7−→ z(t) := z0 + νuv (s) ds t0
is differentiable (and thus, continuous), but not monotone, meaning that even if z(t0 ) = z0 ∈ Ωuv ′ and z(t) ∈ Ωuv holds for a t ∈ R, this does not automatically mean that z(t ) ∈ Ωuv holds for all ′ t ∈ [min{t0 , t}, max{t0 , t}], too. However, since Equation (2) only holds on the open sub-domain R × Ωuv of ´t the function cuv : R × [0, luv ] → R≥0 , any path from z(t0 ) = z0 to z(t) = z0 + t0 νuv (s) ds parametrized by time needs to be fully included in the space interval Ωuv = (0, luv ). In other words, t ∈ Tt0 ,z0 needs to hold. Since the set Tt0 ,z0 can also be empty, we are interested in when this is not the case and what properties it satisfies in this case. The following Lemma C.3 answers these questions. Lemma C.3 (On Tt0 ,z0 ). Let euv ∈ E, t0 ∈ R, z0 ∈ R hold and let Tt0 ,z0 be defined as in Lemma 3.2. 28
1. If z0 ̸∈ Ωuv holds, then Tt0 ,z0 = ∅ holds. 2. If z0 ∈ Ωuv holds, then Tt0 ,z0 ̸= ∅ holds. Even more, t0 ∈ Tt0 ,z0 holds. 3. If z0 ∈ Ωuv holds, there exist tl ∈ [−∞, t0 ) and tr ∈ (t0 , ∞] such that ´t (a) if tl ̸= −∞ holds, z0 + t0l νuv (s) ds ∈ {0, luv } holds, ´t (b) if tr ̸= +∞ holds, z0 + t0r νuv (s) ds ∈ {0, luv } holds and (c) Tt0 ,z0 = (tl , tr ) = (t0 − (t0 − tl ), t0 + (tr − t0 )) holds. Equivalently, there exist δtl , δtr ∈ (0, ∞] such that Tt0 ,z0 = (t0 − δtl , t0 + δtr ) holds. Since according to the previous Lemmas 3.2 and C.3, for any t0 ∈ R and z0 ∈ Ωuv , the concentration cuv is constant along characteristic curves γt0 ,z0 on an non-empty time interval Tt0 ,z0 = (tl , tr ), it suffices to know the value of cuv at one position in this time interval to have knowledge of the concentration over all times in the interval: Theorem C.4 (Advective transport along edges). Let euv ∈ E hold. If the function cuv obeys Equation (2), for each t0 ∈ R and z0 ∈ Ωuv , ˆ t cuv t, z0 + νuv (s) ds = cuv (t, z0 + ν uv (t0 , t) (t − t0 )) = cuv (t0 , z0 ) t0
holds for all t ∈ Tt0 ,z0 with Tt0 ,z0 as defined in Lemma 3.2 and for the mean velocity (ffl t ( ´t 1 νuv (s) ds if t ̸= t0 ν (s) ds if t ̸= t0 t−t t t0 uv 0 0 ν uv (t0 , t) := = . 0 if t = t0 0 if t = t0 According to Lemma C.3.3, for any t0 ∈ R and z0 ∈ Ωuv , the set Tt0 ,z0 usually corresponds to a time interval Tt0 ,z0 = (tl , tr ) such that on its boundary ∂(tl , tr ) = {tl , tr }, the information reaches exactly one of the ends of the edge euv , i.e., the node v ∈ V or the node u ∈ N (v). In view of the coupling conditions (Eq. (4)), this motivates to leverage the advective transport along edges, i.e., Theorem C.4, and combine it with the inflow boundary conditions defined through Equation (5), in order to derive a representation of the advection coupling conditions (Eq. (4)) that does not depend on the concentrations cuv (t, luv ) at the end z = luv of the edges euv for each u ∈ N− (v, t), but on the concentrations cu (t − δtuv ) and cv (t − δtuv ) at suitable times t − δtuv ∈ R. This, in turn, will serve as the basis for our message-passing algorithm (cf. Eq. (6) and (9)). Intuitively, the time shift δtuv ∈ R will correspond to the time that any information requires to travel along the edge euv with the flow νuv : R → R from one of the nodes v ∈ V or u ∈ N (v) until it reaches the node v ∈ V again. From the sender perspective, two different scenarios can occur: Either, the information leaves the node u ∈ N (v) at time t − δtuv ∈ R and enters the node v ∈ V at time t ∈ R. We will call this scenario a pass-through. Or, the information leaves the node v ∈ V at time t − δtuv ∈ R and starts traveling to the node u ∈ N (v), but before it reaches that node, the flow νvu = −νuv : R → R changes its sign and the information enters the node v ∈ V at time t ∈ R again. We will call this scenario a self-loop. Of course, the sign of the flow can also change multiple times before the information enters the node v ∈ V . As a first step, we now formalize the aforementioned pass-through and self-loop times formally and collect some helpful properties to get a deeper understanding of such. In practice, also corresponding inverse phenomena will appear. In summary, we call these times transport times, a name which we will motivate later. Definition C.5 (Transport times). Let euv ∈ E and t ∈ R hold. If there exists a δtuv ∈ R, such that ˆ t νuv (s) ds = luv and
(19)
νuv (s) ds ∈ (0, luv )
(20)
t−δtuv
ˆ t′
t−δtuv
29
′
hold for all t ∈ (min{t − δtuv , t}, max{t − δtuv , t}), we call δtuv the pass-through time w.r.t. the edge euv ∈ E and time t ∈ R. If instead, there exists a δtuv ∈ R, such that ˆ t νuv (s) ds = − luv and (21) t−δtuv
ˆ t′
νuv (s) ds ∈ (−luv , 0)
(22)
t−δtuv ′
hold for all t ∈ (min{t − δtuv , t}, max{t − δtuv , t}), we call δtuv the inverse pass-through time w.r.t. the edge euv ∈ E and time t ∈ R. If instead, there exists a δtuv ∈ R̸=0 , such that ˆ t νuv (s) ds = 0 and (23) t−δtuv
ˆ t′
νuv (s) ds ∈ (0, luv ),
(24)
t−δtuv ′
hold for all t ∈ (min{t − δtuv , t}, max{t − δtuv , t}), we call δtuv the self-loop time w.r.t. the edge euv ∈ E and time t ∈ R. If instead, there exists a δtuv ∈ R̸=0 , such that ˆ t νuv (s) ds = 0 and (25) t−δtuv
ˆ t′
νuv (s) ds ∈ (−luv , 0),
(26)
t−δtuv ′
hold for all t ∈ (min{t − δtuv , t}, max{t − δtuv , t}), we call δtuv the inverse self-loop time w.r.t. the edge euv ∈ E and time t ∈ R. Usually, euv ∈ E and t ∈ R are clear from the context and we just say pass-through time, inverse pass-through time, self-loop time or inverse self-loop time. The following lemmas help to get a deeper understanding of the different kinds of transport times. Lemma C.6 (Pass-through time). Let euv ∈ E and t ∈ R hold. If there exists a pass-through time δtuv ∈ R w.r.t. euv ∈ E and t ∈ R, the following properties hold: 1. δtuv ∈ R̸=0 and thus, t − δtuv ̸= t holds. 2. There exists no other pass-through time of the same sign (i.e., there exists either exactly one or no positive pass-through time, and there exists exactly one or no negative pass-through time). 3. The mean velocity 1 ν uv (t) := ν uv (t − δtuv , t) = νuv (s) ds = δt uv t−δtuv t
ˆ t νuv (s) ds t−δtuv
during the time interval [min{t − δtuv , t}, max{t − δtuv , t}] as defined in Theorem C.4 is equal to ν uv (t) = δtluv and thus satisfies uv ν uv (t) > 0 ⇐⇒ δtuv > 0 and ν uv (t) < 0 ⇐⇒ δtuv < 0. q
4. The flows quv (t − δtuv ) and quv (t) linked to the velocity νuv (·) = αvu(·) satisfy vu ( quv (t − δtuv ), quv (t) ≥ 0 if δtuv > 0 . quv (t − δtuv ), quv (t) ≤ 0 if δtuv < 0 30
According to Lemma C.6, if a pass-through time δtuv > 0 is positive, so does ν uv (t) > 0 hold, and there is an average flow from the node u ∈ N (v) to the node v ∈ V during the time [t − δtuv , t].4 If in contrast, a pass-through time δtuv < 0 is negative, so does ν uv (t) < 0 hold, and there is an average flow from the node luv v ∈ V to the node u ∈ N (v) during the time [t, t − δtuv ].5 In both cases, the time δtuv = ν uv (t) is exactly the time that any particle transported by the average flow q uv (t) = αuv ν uv (t) needs to travel along the edge euv of length luv . Therefore, we call δtuv the pass-through time. The cases are intuitively connected by the equality6 ˆ t νuv (s) ds = t−δtuv
ˆ t−δtuv
−νuv (s) ds =
t
ˆ t−δtuv
νvu (s) ds,
t
where in dependence of the sign of the pass-through time δtuv ∈ R̸=0 , we use the oriented of the two integrals on the left- or right-hand side for interpretation. We can make similar observations for self-loop times. Lemma C.7 (Self-loop time). Let euv ∈ E and t ∈ R hold. If there exists a self-loop time δtuv ∈ R̸=0 w.r.t. euv ∈ E and t ∈ R, the following properties hold: 1. There exists no other self-loop time of the same sign (i.e., there exists either exactly one or no positive self-loop time, and there exists exactly one or no negative self-loop time). 2. The mean velocity t
ν uv (t) := ν uv (t − δtuv , t) =
νuv (s) ds = t−δtuv
1 δtuv
ˆ t νuv (s) ds t−δtuv
during the time interval [min{t − δtuv , t}, max{t − δtuv , t}] as defined in Theorem C.4 satisfies ν uv (t) = 0. 3. The flows quv (t − δtuv ) and quv (t) linked to the velocity νuv (·) = (
quv (t − δtuv ) ≥ 0, quv (t) ≤ 0 quv (t − δtuv ) ≤ 0, quv (t) ≥ 0
qvu(·) αvu
satisfy
if δtuv > 0 . if δtuv < 0
According to Lemma C.7, independent on whether a self-loop time δtuv ∈ R̸=0 is positive or negative, ν uv (t) = 0 holds. If a self-loop time δtuv > 0 is positive, there is an average flow from the node u ∈ N (v) to itself during the time [t − δtuv , t].7 If in contrast, a self-loop time δtuv < 0 is negative, there is an average flow from the node u ∈ N (v) to itself during the time [t, t − δtuv ].8 In both cases, the time δtuv is exactly the time that any particle transported by the average flow q uv (t) = αuv ν uv (t) = 0 needs to travel from u along some part of the edge euv until the flow changes its sign and the particle travels back to the node u it came from. Therefore, we call δtuv the self-loop time. As another result that will turn out to be useful in practise, by the symmetry of the underlying graph G and the conservation of flows (Eq. (3)), an inverse pass-through time induces a pass-through time and vice versa, while an inverse self-loop time induces a self-loop time and vice versa. 4 More precisely, since in this case, q (t − δt ), q (t) ≥ 0 holds, we know that at time t − δt uv uv uv uv ∈ R and t ∈ R, there is a flow from u ∈ N− (v, t − δtuv ) ∩ N− (v, t) to v ∈ V . Note that we tacitly ignored the case of zero flows here due to Remark E.1. 5 More precisely, since in this case, q (t − δt ), q (t) ≤ 0 holds, we know that at time t ∈ R and t − δt uv uv uv uv ∈ R, there is a flow from v ∈ V to u ∈ N+ (v, t) ∩ N+ (v, t − δtuv ). Note that we tacitly ignored the case of zero flows here due to Remark E.1. 6 The proof follows similar arguments as given in, e.g., the proof of Lemma C.8. 7 More precisely, since in this case, q (t − δt ) ≥ 0 and q (t) ≤ 0 holds, we know that at time t − δt uv uv uv uv ∈ R, there is a flow from u ∈ N− (v, t − δtuv ) to v ∈ V , while at time t ∈ R, there is a flow from v ∈ V to u ∈ N+ (v, t). Note that we tacitly ignored the case of zero flows here due to Remark E.2. 8 More precisely, since in this case, q (t) ≥ 0 and q (t − δt ) ≤ 0 holds, we know that at time t ∈ R, there is a flow from uv uv uv u ∈ N− (v, t) to v ∈ V , while at time t − δtuv ∈ R, there is a flow from v ∈ V to u ∈ N+ (v, t − δtuv ). Note that we tacitly ignored the case of zero flows here due to Remark E.2.
31
Lemma C.8 (Inverse pass-through time). Let euv ∈ E and t ∈ R hold. δtuv ∈ R is an inverse pass-through time w.r.t. euv ∈ E and t ∈ R if and only if (iff) δtvu := δtuv ∈ R is a pass-through time w.r.t. evu ∈ E and t∈R Lemma C.9 (Inverse self-loop time). Let euv ∈ E and t ∈ R hold. δtuv ∈ R is an inverse self-loop time w.r.t. euv ∈ E and t ∈ R iff δtvu := δtuv ∈ R is a self-loop time w.r.t. evu ∈ E and t ∈ R We now collected all relevant properties which help us to prove the main results of this section as motivated above: Given a pass-through time δtuv ∈ R w.r.t. euv ∈ E and t ∈ R, we can re-write the concentrations cuv (t, luv ) and cvu (t − δtuv , lvu ) at the end z = luv = lvu of the edges euv and evu using the concentrations cu : R → R≥0 and cv : R → R≥0 , respectively: Theorem C.10 (Pass-through time). Let euv ∈ E and t ∈ R hold. If the function cuv obeys Equation (2) and there exists a pass-through time δtuv ∈ R̸=0 w.r.t. euv ∈ E and t ∈ R, we obtain luv ,0 . cuv (t, luv ) = cuv (t − δtuv , 0) = cuv t − ν uv (t) Even more, the concentrations cuv (·, luv ) and cvu (·, lvu ) at the end z = luv = lvu of the edges euv ∈ E and evu ∈ E are connected to the concentrations cu and cv of the nodes of that edge, respectively, by ( cuv (t, luv ) = cu (t − δtuv ) if δtuv > 0 . cvu (t − δtuv , lvu ) = cv (t) if δtuv < 0 According to Theorem C.10, if a pass-through time δtuv > 0 is positive, the concentration cuv (t, luv ) at time t ∈ R and at the end of the edge euv ∈ E, along which there is an average flow from the node u ∈ N (v) to the node v ∈ V during the time [t − δtuv , t], is equal to the concentration cu (t − δtuv ) at the node u ∈ N (v) and −δtuv < 0 units of time earlier than t. If in contrast, a pass-through time δtuv < 0 is negative, the concentration cvu (t − δtuv , lvu ) at time t − δtuv ∈ R and at the end of the edge evu ∈ E, along which there is an average flow from the node v ∈ V to the node u ∈ N (v) during the time [t, t − δtuv ], is equal to the concentration cv (t) at the node v ∈ V and +δtuv < 0 units of time earlier than t − δtuv . Analogously, given a self-loop time δtuv ∈ R w.r.t. euv ∈ E and t ∈ R, we can re-write the concentration cvu (t, lvu ) at the end z = lvu = luv of the edge evu using the concentration cu : R → R≥0 : Theorem C.11 (Self-loop time). Let euv ∈ E and t ∈ R hold. If the function cuv obeys Equation (2) and there exists a self-loop time δtuv ∈ R̸=0 w.r.t. euv ∈ E and t ∈ R, we obtain cuv (t, 0) = cuv (t − δtuv , 0). Even more, the concentration cvu (·, lvu ) at the end z = lvu = luv of the edge evu ∈ E is connected to the concentration cu of the node of that edge by ( cvu (t, lvu ) = cu (t − δtuv ) if δtuv > 0 . cvu (t − δtuv , lvu ) = cu (t) if δtuv < 0 According to Theorem C.11, if a self-loop time δtuv > 0 is positive, the concentration cvu (t, lvu ) at time t ∈ R and at the end of the edge evu ∈ E, along which there is an average flow from the node u ∈ N (v) to itself during the time [t − δtuv , t], is equal to the concentration cu (t − δtuv ) at the node u ∈ N (v) and −δtuv < 0 units of time earlier than t. If in contrast, a self-loop time δtuv < 0 is negative, the concentration cvu (t − δtuv , lvu ) at time t − δtuv ∈ R and at the end of the edge evu ∈ E, along which there is an average flow from the node u ∈ N (v) to itself during the time [t, t − δtuv ], is equal to the concentration cu (t) at the node u ∈ N (v) and +δtuv < 0 units of time earlier than t − δtuv . In summary, Theorem C.10 and Theorem C.11 illustrate the time shift required to pass information from one node to another. Therefore we call these shifts transport times, and the underlying concepts serve 32
as a basis for our advection-message-generation function ϕ (cf. Eq. (6)). To be able to properly define it, we need to make sure that for each edge euv ∈ E and time t ∈ R, such a transport time exists and is unique. For the sake of brevity, for such a pair, we define the function ˆ t ý : R −→ R, δt 7−→ ý(δt) := νuv (s) ds (Equation (7)), t−δt
where we omit the dependency of ý on euv ∈ E, t ∈ R and the overall flow field ν for better readability. Then we obtain the following existence and uniqueness result: Theorem C.12 (Existence and uniqueness of positive transport times). Let euv ∈ E and t ∈ R hold. If νuv (t) ̸= 0 holds and the set ˆ t νuv (s) ds ∈ {−luv , 0, luv } Auv (t) := δt ∈ (0, ∞) t−δt = δt ∈ (0, ∞) ý(δt) ∈ {−luv , 0, luv } is not empty, the following properties hold: 1. δtuv := min Auv (t) = min{δt ∈ (0, ∞)
ý(δt) ∈ {−luv , 0, luv }} ∈ (0, ∞) exists.
2. δtuv is a positive pass-through time, inverse pass-through time, self-loop time or inverse self-loop time. More specifically, δtuv is (always w.r.t. euv and t) ... ... an inverse pass-through time
⇐⇒ ý(δtuv ) = −luv ,
... a self-loop time or an inverse self-loop time
⇐⇒ ý(δtuv ) = 0,
... a pass-through time
⇐⇒ ý(δtuv ) = luv .
3. There exists no other positive pass-through time, inverse pass-through time, self-loop time or inverse self-loop time than δtuv . Theorem C.12 guarantees that the mapping ù : R × E × F (R, Rne ) −→ (0, ∞), (t, euv , ν) 7−→ ù(t, euv , ν) := δtuv
(27)
is well-defined. Throughout most parts of this work, t ∈ R, euv ∈ E and ν = (νe )e∈E are known from the context and we use the notation δtuv as a shortcut for ù(t, euv , ν) to ease notation whenever it is reasonable. Finally, this allows to define the advection-message-generation function in accordance to the results of Theorems C.10 and C.11 as ϕ : R × E × F (R, Rnn ) × F (R, Rne ) −→ R≥0 , ( cu (t − δtuv ) (t, euv , c, ν) 7−→ ϕ(t, euv , c, ν) := cv (t − δtuv )
if ý(δtuv ) ∈ {−luv , luv } if ý(δtuv ) = 0
(Equation (9)).
Remark C.13 (Double assignment phenomenon). Note that through the relation9 ˆ t ˆ t νvu (s) ds = − νuv (s) ds, t−δtuv
t−δtuv
two neighboring nodes v ∈ V and u ∈ N (v) will be assigned to the same transport time δtuv , corresponding to either pass-through and inverse pass-through times (cf. Lemma C.8) or self-loop and inverse self-loop times (cf. Lemma C.9). We call this the double-assignment phenomenon. 9 See, e.g., the proof of Lemma C.8.
33
In the following Section C.1.2, we integrate our advection-message-generation function ϕ into the coupling conditions (Eq. (4)) and show that the resulting representation of such corresponds to our advectionmessage-passing algorithm as presented in Section 3 (Eq. (10)) and which automatically handles the doubleassignment phenomenon. C.1.2
Message Aggregation and Update: Mixing at Nodes
As investigated in the previous Section C.1.1, for a node v ∈ V and a time t ∈ R, we are interested in a representation of the advection coupling conditions (Eq. (4)) that does not depend on the concentrations cuv (t, luv ) at the end z = luv of the edges euv for each u ∈ N− (v, t), but on the concentrations cu (t − δtuv ) and cv (t − δtuv ) shifted by suitable transport times δtuv ∈ R. Since in practice, we can only assume to have knowledge of the past, the suitability refers to the requirement that δtuv > 0 should hold. Theorem C.12 guarantees the existence and uniqueness of such. In this subsection, we prove that the message-generation function as defined in Equation (9) indeed assigns to each time t ∈ R and each edge euv ∈ E a message ϕ(t, euv , c, ν) that determines the correct representation of cuv (t, luv ) from the correct source under consideration of the concentrations c and the flow field ν. The following lemma will help reaching this goal. Lemma C.14 (Transport times). Let euv ∈ E and t ∈ R hold. In the setting of Theorem C.12, the following property holds: 1. δtuv is (always w.r.t. euv and t) ... ... a pass-through time
⇐⇒ ý(δtuv ) ∈ {−luv , luv } and quv (t) > 0,
... an inverse pass-through time
⇐⇒ ý(δtuv ) ∈ {−luv , luv } and quv (t) < 0,
... a self-loop time
⇐⇒ ý(δtuv ) = 0 and quv (t) < 0,
... an inverse self-loop time
⇐⇒ ý(δtuv ) = 0 and quv (t) > 0.
Lemma C.14 justifies why in the definition of the advection-message-generation function ϕ (Eq. (9)), we do not need to distinguish between pass-through and inverse pass-through or self-loop and inverse self-loop times and though automatically handle the double assignment phenomenon (cf. Remark C.13): Taking into account the sign of the flow rate quv (t) at time t ∈ R, corresponding counterparts will vanish in the coupling condition due to the summation over inflow neighbors u ∈ N− (v, t) only (cf. Eq. (4)) and in our advection-message-passing algorithm due to the ReLU-function (cf. Eq. (6) and (10)). Moreover, the assumptions of Theorem C.12 and Lemma C.14 are naturally satisfied in the contexts where they are applied: On the one hand, summands where the velocity νuv (t) is equal to zero do not contribute by construction of the coupling conditions (4). On the other hand, whenever the set Auv (t) is empty (that is, no transport time exists), we assume to have knowledge of the values cv (t) for the respective times t: In cases of initial value problems, these values are known by definition (cf. Remark 3.8); for boundary value problems we assume an initial state of zero at non-boundary nodes, these values cv (t) for the respective times t are then also zero by definition. Consequently, we can combine the results from the advective transport along edges (Section C.1.1) with the results of this Section to obtain a message-passing representation of the coupling conditions (Eq. (4)) which indeed coincides with our final advection-message-passing scheme (10): Theorem C.15 (Physical alignment (full version)). Let v ∈ V and t ∈ R hold. We assume that Assumption 3.4 holds and that for each u ∈ N (v), Auv (t) is not empty (that is, the OP (8) has a solution). If the function cv obeys the mixing at nodes (Eq. (4)) and if for each u ∈ N (v), the function cuv obeys the advection edge dynamics (Eq. (2)) with inflow continuity condition (Eq. (5)), then cv (t) =
X u∈N (v) ý(δtuv )=±luv
wuv (t) · cu (t − δtuv ) +
X
wuv (t) · cv (t − δtuv ) =
u∈N (v) ý(δtuv )=0
34
X u∈N (v)
wuv (t) · φ(t, euv , c, ν)
holds, where δtuv = min Auv (t) is the transport time as defined in Theorem C.12. In other words, Equation (10) with ϕ as defined in Equation (9) holds. C.2
Error Bounds on MeGA-MP
In this section, we first provide some theory on a special graph structure, trees, and introduce further notation along the way (Section C.2.1). This does not only allow to formalize the full version of Theorem 3.7, but also to prove Theorem 3.6 and 3.7 (Section C.2.2). Nevertheless, the following sections can also help to get a deeper understanding of our iterative algorithm, Equation (11), whose intuition was already describe in Section 3.2. C.2.1
Preliminaries: Trees with Children-Aggregation Scheme
In this section, we first investigate on the intuition behind the message passing scheme (10) (and analogously, Eq. (11)), which – by the definition of the weights wuv (Eq. (6)) – can also be written as X wuv (t) · ϕ(t, euv , c, ν). cv (t) = u∈N− (v,t)
For simplicity, let us assume that the message function ϕ outputs ϕ(t, euv , c, ν) = cu (t − δtuv ) (cf. Eq. (9)), that is, that at time t, we observe a pass-through along the edge euv and not a an inverse self-loop (cf. Lemma C.14). In this case, Equation (10) can – with (1) a change in notation and (2) using the more formal notation from Equation (27) – be written as (1)
cv (t) =
X
(2)
wv1 v (t) · cv1 (t − δtv1 v ) =
v1 ∈N− (v,t)
X v1 ∈N− (v,t)
wv1 v (t) · cv1 (t − ù(t, ev1 v , ν)). {z } | =:t1
For each inflow neighbor v1 ∈ N− (v, t) of v at time t, we can again apply Equation (10), yielding X X cv (t) = wv1 v (t) · wv2 v1 (t1 ) · cv2 (t1 − ù(t1 , ev2 v1 , ν)). {z } | v1 ∈N− (v,t) v2 ∈N− (v1 ,t1 )
=:t2
Intuitively, we can repeat this process as often as we want until a neighborhood N− (vl , tl ) is empty, which by definition is the case if vl ∈ Vb is a source node. What we create along the way are paths p = ((v0 , t0 ) := (v, t), (v1 , t1 ), (v2 , t2 ), ...) of both locations vl ∈ V and times tl ∈ R for indices l = 1, 2, ... along which we trace back information that contributes to the signal cv (t).10 Knowledge of such paths, which we call inflow paths, can be helpful in case where the signal cv (t) at time t is unknown, expressing it through signals cvl (tl ) at earlier times tl < t. While the graph G is a symmetric directed graph (cf. Section 2), the set of all such inflow paths P(v, t) defines a tree structure11 , which we construct formally and iteratively from the root node (v0 , t0 ) = (v, t) to leaf nodes in the next Section C.2.2. Before, we provide a general result on trees with children-aggregation aggregation schemes: Theorem C.16 (Trees with children-aggregation scheme). Let T = (A, E) be a finite tree with root node a0 ∈ A, internal nodes Aint and leaf nodes Aleaf (that is, A = {a0 } ⊔ Aint ⊔ Aleaf holds). If each node a ∈ A is associated with a magnitude xa ∈ R that satisfies X xa = ωab · xb for all a ∈ {a0 } ⊔ Aint , b∈C(a)
where the set C(a) ⊂ A consists of the children of a ∈ A and wab is a weight of the edge eab ∈ E that depends on both the parent a ∈ A and its child b ∈ C(a), then the magnitude x0 := xa0 of the root node a0 ∈ A is 10 Note that the indices l = 1, 2, ... do not necessarily correspond to the discrete times indices I in Section 3.2. 11 Note that by definition, the transport time is positive (Theorem C.12). Therefore, t > t > t > ... holds. 1 2
35
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
...
Contribution of the path for any
to
:
Figure 11: A visualization of the concepts from Theorem C.16. determined by the magnitudes xlp := xalp of the leaf nodes alp ∈ Aleaf by X
x0 =
ω(a0 , ..., alp ) · xlp ,
p=(a0 ,...,alp )∈P(a0 ,Aleaf )
where P(a0 , Aleaf ) is the set of paths connecting the root node a0 to one of its leaf nodes Aleaf and lp −1
ω(a0 , ..., alp ) :=
Y
ωal al+1
l=0
is the product of weights along a(n) (arbitrary) path p = (a0 , ..., alp ) ∈ P with length lp ∈ N. A visualization of such a tree can be found in Figure 11. C.2.2
Construction of Inflow Trees
In this section, we first provide some intermediate results – Theorem C.19 and Theorem C.20 – that help proving Theorem 3.6 and 3.7. In order to prove the former, we construct suitable trees in the sense of Theorem C.16 as already motivated in Section C.2.1: 36
Definition C.17 (Inflow tree and interpolation inflow tree). Let G = (V, E) be an advection metric graph as defined in Section 2. Let v ∈ V and t ∈ R hold. We define the inflow tree T (v, t) of v at time t using the notation from Theorem C.16 by iteratively defining the nodes and children per layer l ∈ N0 of the tree until a stopping criteria is fulfilled: • l = 0-th layer: We choose the root node a0 = (v0 , t0 ) := (v, t). • l + 1-th layer: Given a parent node al = (ul , vl , tl )12 in the l-th layer, we choose its children al+1 = (ul+1 , vl+1 , tl+1 ) ∈ C(al ) as each combination of – an inflow neighbor ul+1 ∈ N− (vl , tl ) of vl at time tl and – the corresponding shifted time tl+1 := tl − ù(tl , eul+1 vl , ν) obtained through the time shift applied to tl using the (positive) transport time ù(tl , eul+1 vl , ν) > 0 (cf. Eq. (27)) along the edge eul+1 vl ∈ E and – the corresponding source node13 ( ul+1 if ý(ù(tl , eul+1 vl , ν)) ̸= 0, vl+1 := (28) vl if ý(ù(tl , eul+1 vl , ν)) = 0. • We stop the construction after an al+1 = (ul+1 , vl+1 , tl+1 ) if a stopping criteria is fulfilled, which for now can be any stopping criteria, or until naturally, the neighborhood N− (vl+1 , tl+1 ) is empty, which is by definition the case if vl+1 ∈ Vb is a source node. • We denote the set of paths connecting the root node a0 = (v0 , t0 ) = (v, v, t) to one of its leaf nodes Aleaf by P(v, t), which is an abbreviation for P(a0 , Aleaf ). We similarly define the interpolation inflow tree T̂ (v, t) of v at time t using the notation from Theorem C.16 by iteratively defining the nodes and children per layer l ∈ N0 of the tree until a stopping criteria is fulfilled: • l = 0-th layer: We choose the root node a0 = (v0 , t0 ) := (v, t). • l + 1-th layer: Given a parent node al = (ul , vl , tl ) in the l-th layer, we choose its children al+1 = (ul+1 , vl+1 , tl+1 ) ∈ C(al ) as each combination of – an inflow neighbor ul+1 ∈ N− (vl , tl ) of vl at time tl and – an corresponding interpolation time tl+1 ∈ {tj , tj+1 } for the j ∈ I for which tl −ù(tl , eul+1 vl , ν) ∈ [tj , tj+1 ) holds (cf. Eq. (17))14 and – the corresponding source node vl+1 analogously to Equation (28).15 • We stop the construction after an al+1 = (ul+1 , vl+1 , tl+1 ) for which vl+1 ∈ Vb or tl+1 ∈ Thi = {t1 , ..., ti } holds.16 • We denote the set of paths connecting the root node a0 = (v0 , t0 ) = (v, t) to one of its leaf nodes Aleaf by P̂(v, t), which is an abbreviation for P̂(a0 , Aleaf ). 12 Note that the definition of a l+1 does not depend on ul , which is why this step is well-defined for l = 0 and l + 1 = 1, too. 13 In the former case, we observe a pass-through along the edge e ul+1 vl , that is, information flows from vl+1 = ul+1 to vl .
In the latter case, we observe an inverse self-loop along the edge eul+1 vl , that is, information flows from vl+1 = vl to vl . For more details, see Section C.1. 14 Similar to the construction to the inflow tree T (v, t), t − ù(t , e l l ul+1 vl , ν) is the – possibly not in the observable, discrete subdomain of R lying – shifted time obtained through the time shift applied to tl using the (positive) transport time ù(tl , eul+1 vl , ν) > 0 (cf. Eq. (27)) along the edge eul+1 vl ∈ E. Note that for l ∈ N, tl already corresponds to one of the discrete time points that are known, that is, there exists an j ∈ I such that tl = tj holds. 15 Note that the definition of v l+1 does not depend on tl+1 , which is why the corresponding source node vl+1 of both tl+1 = tj and tl+1 = tj+1 is the same. 16 Note that similar to the construction to the inflow tree T (v, t), we can also chose any other suitable stopping criteria. We chose the one that we will need later. Also note that with this stopping criterium, if v ∈ Vb or t ∈ Thi holds, the creation of the tree T̂ (v, t) stops immediately and does only consist of the path p = ((v0 , t0 ) = (v, t)) of length zero.
37
Remark C.18 (Finite depths of interpolation inflow tree). We want to explicitly highlight that for each v ∈ V and t ∈ R, the stopping criteria for the interpolation inflow tree T̂ (v, t) is fulfilled at some point. In other words, the interpolation inflow tree T̂ (v, t) has a finite depth and each path p = ((v0 , t0 ) = (v, t), (u1 , v1 , t1 ), ..., (ulp , vlp , tlp )) ∈ P̂(v, t) has a finite length lp ∈ N. This is because by definition of the OP (8), we assume positive transport times (also see discussion above Theorem C.15). Consequently, t = t0 > t1 > t2 > ... holds and for arbitrary t ∈ R, at some point, tlp ∈ Thi must hold. We now leverage these trees together with Theorem C.16 to prove the following two intermediate results, (k) which allow to reconstruct the magnitudes cv (t) and ĉv (t) of interest along inflow paths: Theorem C.19 (Closed form of advection message passing). Let v ∈ V and t ∈ R hold. Let T (v, t) be the inflow tree of v at time t as defined in Definition C.17. If cv (t) is defined by Equation (10), then lp −1 X Y cv (t) = wul+1 vl (tl ) · cvlp (tlp ) (29) p=((v0 ,t0 )=(v,t),(u1 ,v1 ,t1 ),...,(ulp ,vlp ,tlp ))∈P(v,t)
l=0
with weights wul+1 vl (tl ) according to Equation (6) holds. Theorem C.20 (Closed form of advection message passing iterations). Let v ∈ V and t ∈ T = Thi ∪Tfu hold. Let T̂ (v, t) be the interpolation inflow tree of v at time t as defined in Definition C.17. Let k ≥ maxp∈P̂(v,t) lp (k)
hold. If ĉv (t) is defined by Equation (11), then
lp −1
ĉ(k) v (t) =
X
Y
p=((v0 ,t0 )=(v,t),(u1 ,v1 ,t1 ),...,(ulp ,vlp ,tlp ))∈P̂(v,t)
p) wul+1 vl (tl ) · a(tl+1 ) · ĉ(k−l (tlp ) vlp
(30)
l=0
with weights wul+1 vl (tl ) and a(tl+1 ) according to Equation (6) and (16), respectively, holds. Theorem C.20 already allows a proof of Theorem 3.6, which is given together with the other proofs in Appendix E. For the proof of Theorem 3.7, we need to make some further assumptions: While Equation (29) and (30) are already of similar structure, the number of paths p̂ ∈ P̂(v, t) grows exponentially by a factor of 2l , where l is the depth of the tree T (v, t) and T̂ (v, t), as compared to the number of paths p ∈ P (v, t). This is due the fact that per time that we observe in the inflow tree T (v, t), we observe two interpolation times in the interpolation inflow tree T̂ (v, t) (cf. Definition C.17). To make each summand in Equation (30) comparable to a summand in Equation (29), we need to make the following assumptions: Assumption C.21 (Identifiable paths). Let v ∈ V and t ∈ R hold. Let T (v, t) and T̂ (v, t) be the inflow tree and the interpolation inflow tree from Definition C.17, respectively. We assume that each path p̂ = ((v̂0 , t̂0 ) = (v, t), (û1 , v̂1 , t̂1 )..., (ûlp , v̂lp , t̂lp )) ∈ P̂(v, t) is uniquely identifiable with a path p = ((v0 , t0 ) = (v, t), (u1 , v1 , t1 ), ..., (ulp , vlp , tlp )) ∈ P(v, t) in the sense that ûl = ul and v̂l = vl holds for all l = 1, ..., lp . Assumption C.22 (Equal time for jumping into history). Let v ∈ V and t ∈ R hold. Let T̂ (v, t) be the interpolation inflow tree from Definition C.17. In the definition of each tl+1 ∈ {tj , tj+1 }, we assume that tj , tj+1 ∈ Tfu or tj , tj+1 ∈ Thi holds. Assumption C.21 is satisfied as long as the inflow neighborhoods N− (vl , tl ), from which we choose an inflow neighbor ul+1 when creating a path in T (v, t), are equal to the corresponding neighborhoods N− (v̂l , t̂l ), from which we choose an inflow neighbor ûl+1 when creating a path in T̂ (v, t). This is guaranteed if the underlying flow field ν that defines these neighborhoods does not change its orientation during the time intervals [tj , tj+1 ] for which tl ∈ [tj , tj+1 ) holds. If the assumption holds, it allows the following definition. Definition C.23 (Joint paths). Let v ∈ V and t ∈ R hold. Let T (v, t) and T̂ (v, t) be the inflow tree and the interpolation inflow tree from Definition C.17, respectively. Let assumption C.21 hold. If the path p̂ = ((v̂0 , t̂0 ) = (v, t), (û1 , v̂1 , t̂1 )..., (ûlp , v̂lp , t̂lp )) ∈ P̂(v, t) is uniquely identifiable with the path p = ((v0 , t0 ) = (v, t), (u1 , v1 , t1 ), ..., (ulp , vlp , tlp )) ∈ P(v, t), we define the joint path of p̂ by p̃ = ((v0 , t0 ) = (v, t), (u1 , v1 , t1 , t̂1 ), ..., (ulp , vlp , tlp , t̂lp )). 38
(31)
We can still consider p̃ as an element of T̂ (v, t), since we only added information that are unused, in, e.g., Theorem C.20. Moreover, Assumption C.22 guarantees that there exist two such paths that are exactly equal up to the last t̂lp ∈ {tj , tj+1 } for a j ∈ I in the sense of Definition C.17. Finally, using the property that per such two paths ending in such two left and right interpolation bounds, the factors a(tj ) and a(tj+1 ) add up to one, we can canonically extend the sum in Equation (29) to the sum in Equation (30) using the joint paths. This allows to prove Theorem 3.7, which in its full version is phrased as follows: Theorem C.24 (Interpolation error (full version)). Let v ∈ V and t ∈ T = Thi ∪ Tfu hold. Let T (v, t) and T̂ (v, t) be the inflow tree and the interpolation inflow tree of v at time t as defined in Definition C.17, respectively. Let Assumption C.21 and C.22 hold. Then the prediction error in-between cv (t) and ĉv (t) is given by |cv (t) − ĉv (t)|
≤
X p̃∈P̂(v,t)
a(t̂1 , ..., t̂lp ) · |cvlp (tlp ) − cvlp (t̂lp )| +|cvlp (tlp )| · |w(t0 , ..., tlp −1 ) − w(t̂0 , ..., t̂lp −1 )| | {z } {z } | Temporal interpolation error
Accumulated weighting error
with joint paths p̃ = ((v0 , t0 ) = (v, t), (u1 , v1 , t1 , t̂1 ), ..., (ulp , vlp , tlp , t̂lp )) ∈ P̂(v, t) according to Definition C.23, lp −1
lp −1
a(t̂1 , ..., t̂lp ) :=
Y
a(t̂l+1 )
w(t̂0 , ..., t̂lp −1 ) :=
and
Y
wul+1 vl (t̂l )
l=0
l=0
with weights wul+1 vl (tl ) and a(tl+1 ) according to Equation (6) and (16), respectively. If the temporal interpolation error is zero, we obtain |cv (t) − ĉv (t)| = 0. Theorem C.24 states that the error between the ground truth cv (t) and the prediction ĉv (t) solely stems from the time interpolation: On the one hand, the temporal interpolation error is the error from observing a value cvlp (t̂lp ) at a node vlp ∈ V and at a discrete time tlp ∈ Thi , when actually, we needed to observe the value cvlp (tlp ) at an eventually non-discrete time tlp ∈ R. On the other hand, for exactly the same reason, the weights wul+1 vl (t̂l ) that we observe at discrete times t̂l ∈ Thi ⊔ Tfu might differ from the weights wul+1 vl (tl ) that we needed to observe at non-discrete times tl ∈ R for all l = 1, ..., lp − 1, yielding the accumulated weighted error. Both errors decrease with decreasing content of high frequencies in concentrations c and flow fields ν.
D
Details on the Experiments
D.1
Experiment 1: Zero-Shot Transfer Learning of Advection-Reaction Dynamics on Real-World Metric Graphs
Domain In both networks Hanoi and L-Town, we model chlorine decay, which is part of the edge dynamics in Equation (12), as a first-order reaction over time and also include interactions with the pipe wall, which depend on the pipe diameter ρuv ∈ R≥0 . This is formally defined as 4kw fr (t, z, euv ) = − kb + cuv (t, z), ρuv
(32)
where the parameters kb ∈ R≥0 and kw ∈ R≥0 control the intensity of chlorine decay over time and due to wall reactions, respectively. 39
Dataset We generate datasets of 1000 training, 200 validation, and 200 testing samples for the Hanoi WDS, and 200 testing samples for the L-Town WDS. Each sample consists of n = 701 time steps with a time difference of dt = 60 (in s) and corresponds to a time series T that describes the distribution of chlorine c over discrete time points t ∈ T and through the WDS G, determined by advection or advection-reaction dynamics as described in section 4.1. The chlorine is injected over time t ∈ T at boundary nodes Vb . The samples differ in graph topology, in flow fields and injection pattern (boundary condition). We now describe this data generation process in more detail. As stated in Section 4, we use EPyT-Flow (Artelt et al., 2024) based on EPANET-MSX (Shang et al., 2023) as the simulator. • Fixed parameters: Boundary nodes Vb at which chlorine is injected correspond to reservoir nodes Vr . This corresponds to one and two boundary nodes in the case on Hanoi and L-Town, respectively. For the advection-reaction datasets, the values of the reaction coefficients in Equation (32) that model the decay of the concentration are set to kb = 0.04 and kw = 0.034. • Varying graph topology: In order to vary graph topology, we sample pipe lengths luv ∼ U (0.1, 80) (in m) and pipe diameters ρuv ∼ U (0.06, 0.025) (in m) uniformly per sample at random. To maintain the typical WDS topology of larger pipes at reservoirs and successively smaller pipes towards leaf nodes, diameters are ordered from low to high, as in the original WDSs Hanoi and L-Town. • Varying flow fields: In order to generate scenarios with different flow fields, we vary nodal demands d(t) = (dv (t))v∈V that directly determine the flow field ν(t) = (νe (t))e∈E through hydraulic equations (Rossman et al., 2020). Empirically, independent random demands yield spatially homogeneous flow fields. To introduce spatial variations of the flow fields, we apply spectral clustering s : V → {0, ..., nc } on the binary network adjacency matrix to partition the network into nc = 4 clusters for Hanoi and nc = 7 clusters for L-Town. Then the demand pattern for a node v ∈ V is sampled from a cluster-shared Gaussian process plus independent noise: (i − j)2 2 . dv (t) = fs(v) (t) + ϵv , ϵv ∼ N (0, 0.1 ), fs(v) ∼ GP(0, K), K(i, j) = exp − 2 · 0.12 · (n − 1)2 The resulting hydraulics show flow velocities ν that vary in space and time and can also change direction. Assigning demands dv (t) may result in physically implausible states that generate negative pressures, however, these states can be easily made plausible by offsetting the heads accordingly, which would have no effect on the flows. Therefore, the flows are physically plausible and since the quality simulation only depends on flows and not on pressures, the resulting evolution of chlorine in the system is still valid. • Varying injection patterns: The injection patterns cv (t) (Dirichlet boundary condition) in which chlorine is injected at boundary nodes v ∈ Vb are sampled from the Fourier function space (using the first M modes) as follows: cv (t) =
M X 1 t t a cos 2kπ + b sin 2kπ , k k kα 60 · 700 60 · 700
k=1
cv (t) =
cv (t) − mint∈T ĉv (t) ∈ [0, 1], maxt∈T cv (t) − mint∈T cv (t)
with M ∼ Uint (3, 30), α ∼ U(0.5, 1.5), and ak , bk ∼ N (0, 1), for all v ∈ Vb , t ∈ T . Here Uint corresponds to an integer-valued uniform distribution. Three samples from this function space are visualized in Figure 12. Model The MLP that is intended to learn the reactions consists of two linear layers with 8 neurons each, and a SELU activation (Klambauer et al., 2017) in between. We train the model with the Adam optimizer (Kingma & Ba, 2014) over 250 epochs with a batch size of 256 using a learning rate of 0.001 that is reduced by a factor of 0.2 if the training loss did not decrease for the last 3 epochs. The loss function is the MAE. Boundary nodes are masked. 40
Cl injection (mgl )
1.0 0.8 0.6 0.4 Injection Pattern 1 Injection Pattern 2 Injection Pattern 3
0.2 0.0
0
2
4
6 Time [hours]
8
10
12
Figure 12: Example patterns used for the chlorine injection at boundary nodes Vb . Baselines
The baselines we compare ourselves to are:
• NNConv: As a simple baseline, we use the kernel-based MPNN from Gilmer et al. (2017). The node feature matrix X ∈ R|V |×n with n = nhi + nfu is populated with the complete time series at boundary nodes Vb , as well as the first time step as the known history Thi at all other nodes V \ Vb . The edge feature matrix E ∈ R|E|×3n contains the time series of scaled edge flows, edge capacities, and transport times. Node and edge features are first mapped into a latent space via two separate linear layers. After that, GG-NN convolutional layers (Gilmer et al., 2017) (NNConv implementation from torch-geometric) are applied. The message function is implemented as a two-layer MLP with a ReLU activation between the two layers. The messages are aggregated via sum aggregation and decoded by a single linear layer calibrated manually. • PDE-GNN: Hermes et al. (2025) is a PDE-style GNN which solves initial and boundary value problems and has been proposed for chlorine estimation in WDSs. The model learns a PDE function that is integrated by Euler integration. Calibrated manually, we apply the same scheduler with 0.2 learning rate decay upon training loss plateau over 3 epochs and we set the initial learning rate to 0.001. The time step used for Euler integration is set to dt = 0.25. Some hyperparameters are tuned, other hyperparameters are left as specified in Hermes et al. (2025). • GPSConv: Rampášek et al. (2022) is a graph transformer that we include as a representative of an architecture capable of capturing global dependencies. As positional encodings, we utilize the graph Laplacian eigenvectors, since positional encodings based on random walks would degenerate to zero if the initial conditions are zero. The Laplacian positional encoding is truncated to 30 components, allowing the transfer-learning experiment from Hanoi (32 nodes) to the L-Town (784 nodes), and embedded to 8 dimensions via a linear layer. We tune the message-passing module, which corresponds to the parameter Φ in Table 4, candidates are NNConv, GINE, and GatedGCN (Gilmer et al., 2017; Hu* et al., 2020; Bresson & Laurent, 2017), which are architectures that can accommodate edge features. • A-DGN: Gravina et al. (2023) is a dynamical-systems based GNN, where the message-passing operation acts as a differential operator, similar to an ODE. As such, it computes updates ∂X ∂l of the node features X that are integrated over several layers l via Euler integration with a step size of ϵ. This model allows deep architectures that are stable and non-dissipative, mediated through antisymmetric weight matrices applied to the output of any message passing layer Φ. We use this model in the weight-sharing style and select the same options for Φ as for GPSConv. As generalization is a focus here, we evaluate two variations A-DGN and A-DGNDia , to account for the difference in size of the training and evaluation graph: The first is a hyperparameter tuned baseline, while for the second variant, we change the number of iterations of the trained model to the graph diameter. We have found that the latter works much better when the step size ϵ is adapted as well. num_iters , and num_iters′ = diameter(G). Specifically, we set ϵ′ = ϵ diameter(G) • PINN-based approaches on metric graphs (Blechschmidt et al., 2022; Laczkó et al., 2025) are effective for learning solutions of specific PDE settings on metric graphs and are suited for inverse problems. However, PINNs are not operator-learning frameworks: A change in PDE parameters or 41
boundary conditions typically requires retraining. Since our focus is on flexible models that generalize across varying flow conditions and topologies, PINN-based approaches are not an appropriate baseline in our setting. • PI DEEPONET is a LEGO-like neural operator for metric graphs (Blechschmidt et al., 2025) originally developed for drift–diffusion edge dynamics. While this approach is highly effective for diffusive dynamics, its direct application to purely hyperbolic problems such as linear advection is not straightforward. First, the method assumes flow-directed edges, effectively imposing a fixed transport direction that is incompatible with more general advection fields that can change direction. Second, at inner nodes the coupling is enforced through Kirchhoff–Neumann conditions, which conserve diffusive gradient fluxes ∂z c. For hyperbolic dynamics, the flux is defined as ν(t, z)c(t, z) which would require a different coupling law such as mixing. For these reasons we omit this work as a baseline for now. All models are trained using the Adam optimizer over a maximum of 1000 epochs. During training we validate the model every 10 epochs and training is stopped early if the validation loss did not improve over the last 50 epochs. For a fair comparison of MeGA-MP and MeGA-MP+ to the other methods, we run a small-scale hyperparameter tuning for every baseline method using Optuna (Akiba et al., 2019). The search space is specified in Table 4. Every configuration is trained for 60 epochs. Trials are terminated at an earlier epoch if the trial’s result is worse than the median of results of previous trials at the same epoch (median pruning). New configurations are sampled based on the tree-structured Parzen Estimator (TPE) (Watanabe, 2025) implemented in Optuna. Table 4: Optuna Search Space for Different Model Types Parameter
NNConv
PDE-GNN
GPSConv
A-DGN
{64, 96, 128, 160}
hidden_dim num_layers
{2, 3, 4, 5, 6}
learning_rate
[10−4 , 3 × 10−3 ] (log)
weight_decay
[10−6 , 10−1 ] (log) {16, 32, 64} (tests with larger batch-sizes were not successful)
batch_size clip_grad_norm
–
{10−4 , 0.5, 1.0, 2.0, −}
–
–
nfu (training)
700
dropout
–
{32, 100, 150}
700
700
[0.0, 0.3]
[0.0, 0.3]
–
num_heads
–
–
{2,4,8}
–
attn_dropout
–
–
[0.0, 0.3]
–
step_size / ϵ
–
{1/4, 1/3, 1/2, 1}
–
{ 1.0, 0.1, 0.01, 0.001 }
γ
–
–
–
{ 1.0, 0.1, 0.01, 0.001 }
Φ
–
–
{ NNConv, GINE, GatedGCN }
{ NNConv, GINE, GatedGCN }
num_iters
–
–
–
{ 1, 2, 5, 10, 20, 50 }
Fixed Hyperparameters nhi
1
nfu (evaluation)
700
loss_fn
L1
epochs
60 (hyperparameter tuning) / 1000 (training, with early stopping)
42
Results As an extension of our ablation studies in Section 4.3, in Figures 13 to 16, we compare the performance of MeGA-MP, MeGA-MP+ and MeGA-MP− on both the advection and the advection-reaction dynamical systems as well as on Hanoi and L-Town qualitatively.
Conc.
Conc.
Conc.
Conc.
MeGA-MP 1.0 1
MeGA-MP +
MeGA-MP
Simulator
2
0.5 0.0 1.0 3
4
0.5 0.0 1.0 5
6
0.5 0.0 1.0 7
8
0.5 0.0
0
2
4
6 Time [hours]
8
10
0
2
4
6 Time [hours]
8
10
Figure 13: Predictions of MeGA-MP, MeGA-MP+ and MeGA-MP− for the advection dynamical system on the WDS Hanoi. The selected nodes are the same as in Figure 5.
Conc.
Conc.
Conc.
Conc.
MeGA-MP 1.0 1
MeGA-MP +
MeGA-MP
Simulator
2
0.5 0.0 1.0 3
4
0.5 0.0 1.0 5
6
0.5 0.0 1.0 7
8
0.5 0.0
0
2
4
6 Time [hours]
8
10
0
2
4
6 Time [hours]
8
10
Figure 14: Predictions of MeGA-MP, MeGA-MP+ and MeGA-MP− for the advection dynamical system on the WDS L-Town. The selected nodes are the same as in Figure 6.
43
Conc.
Conc.
Conc.
Conc.
MeGA-MP 1.0 1
MeGA-MP +
MeGA-MP
Simulator
2
0.5 0.0 1.0 3
4
0.5 0.0 1.0 5
6
0.5 0.0 1.0 7
8
0.5 0.0
0
2
4
6 Time [hours]
8
10
0
2
4
6 Time [hours]
8
10
Figure 15: Predictions of MeGA-MP, MeGA-MP+ and MeGA-MP− for the advection-reaction dynamical system on the WDS Hanoi. The selected nodes are the same as in Figure 5.
Conc.
Conc.
Conc.
Conc.
MeGA-MP 1.0 1
MeGA-MP +
MeGA-MP
Simulator
2
0.5 0.0 1.0 3
4
0.5 0.0 1.0 5
6
0.5 0.0 1.0 7
8
0.5 0.0
0
2
4
6 Time [hours]
8
10
0
2
4
6 Time [hours]
8
10
Figure 16: Predictions of MeGA-MP, MeGA-MP+ and MeGA-MP− for the advection-reaction dynamical system on the WDS L-Town. The selected nodes are the same as in Figure 6. D.2
Experiment 2: Advection on a 1D Euclidean Domain The boundary value problem which we analyse in Section 5 is given by ∂c(t, z) ∂c(t, z) 2πt = −ν(t) with ν(t) = 0.3 · sin + 0.3 and ∂t ∂z 100 ! 2 1 t − µt c(t, 0) = exp 2 σt
Domain and Dataset
(33)
for all t ∈ [0, 100] and z ∈ [0, 100] and with µt = 38 and σt = 18. Its analytical solution is computed analogously to the solution of the initial value problem with boundary condition, as discussed next. 44
Additionally, in this Appendix, we analyze the initial value problem with boundary condition given by ∂c(t, z) ∂c(t, z) 2πt = −ν(t) with ν(t) = 0.3 · sin + 0.3 and ∂t ∂z 100 2 ! 2 ! (34) 1 z − µz 1 t − µt c(0, z) = exp , c(t, 0) = exp 2 σz 2 σt for all t ∈ [0, 100] and z ∈ [0, 100] and with µz = 8, σz = 2.5, µt = 38 and σt = 18. To compute an analytical solution, we can solve the initial value problem civp (t, z) and the boundary value problem cbvp (t, z) separately and sum the individual results. This decomposition is possible since for a positive and incompressible flow ν(t) > 0, characteristic lines do not cross and each point of the solution corresponds to a value of either the initial value function or the boundary value function. The initial value problem is solved by integrating the´ velocity field to compute the spatial shift of the t Gaussian c(0, z) for each t ∈ [0, 100], via t 7−→ ý0 (t) := 0 ν(s) ds. The antiderivative of ν(t) is given by 1 cos(4πt) + , ν̄(t) = 0.3 t − 4π 4π
(35)
and the solution of the initial value problem is civp (t, z) = c(0, z + ý0 (t)). To solve the boundary value problem, we find the transport time δt(z) from the boundary at z = 0 to every other point z in the domain. We leverage the same procedure as for our MeGA-MP (cf. Eq. 8) and compute ´ δt(z) δt(z) as the unique solution of 0 ν(s)ds = z. Since ν̄(t) is strictly increasing the solution is unique and well-defined. According to Theorem C.10, the coordinate (t + δt(z), z) lies on the characteristic curve that intersects the boundary at (t, 0). We can shift the boundary value Gaussian accordingly and obtain cbvp (t, z) = c(t + δt(z), 0). The solution to the initial value problem with boundary condition is then c(t, z) = civp (t, z) + cbvp (t, z) In practice, we compute the respective values numerically at each discrete spatial and temporal grid point. Next to the analytical solution, the numerical and analytical solvers we compare ourselves to are: • Semi-Lagrangian scheme: We use a custom semi-Lagrangian solver. For every grid point z ∈ Ω and iteration step n + 1, the method computes the departure point z ′ of a particle arriving at z by tracing backward along the velocity field ν(t, z). Since the departure point z ′ is a continuous coordinate, we interpolate spatially to compute the new solution c(tn+1 , z). If we assume a spatially constant and temporally piecewise constant flow field, the departure point is simply z ′ = z − dt · ν(tn , z). Otherwise, it is approximated by integrating the characteristic ODE dz dt = −ν(t, z) via the RK4 method. The result is then c(tn+1 , z) = I[c(tn , ·)](z ′ ), where I is the spatial interpolation operator like the linear interpolation in Eq. (16). • Runge-Kutta-4 is the classical order-four Runge-Kutta (RK4) numerical method for solving PDEs (and ODEs) of the kind ∂c(t, z) = f (t, z, c) ∂t with the update rule c(tn+1 , z) = c(tn , z) + 45
dt (k1 + 2k2 + 2k3 + k4 ) 6
where n ∈ N is the iteration step, dt is the step-size and k1 = f (tn , z, c), k2 = f (tn + dt 2 , c(tn , z) + dt dt dt 2 k1 ), k3 = f (tn + 2 , c(tn , z) + 2 k2 ) and k4 = f (tn + dt, c(tn , z) + dtk3 ). We apply the RK4 within a method-of-lines setup. The right-hand side function f (t, z, c) approximates the negative spatial derivative of the flux via a fifth-order Weighted Essentially Non-Oscillatory (WENO5) scheme (Jiang & Shu, 1996) using Lax-Friedrichs flux-splitting for numerical stability. Specifically, left- and right-moving flux components are reconstructed separately at cell interfaces, resulting in the spatial −F̂
F̂
derivative approximation at cell grid-points i: f (t, zi , c) ≈ − i+1/2dx i−1/2 where F̂i±1/2 represent the combined numerical interface fluxes. Baselines and Results In Figure 17, the analytical solution of the boundary value problem, which serves as ground truth, and the MeGA-MP solution are visualized at different discretizations and on the whole resulting grid. Moreover, in Figure 18, the analytical solution of the initial value problem with boundary values, which serves as ground truth, and the MeGA-MP solution are visualized at different discretizations and on the whole resulting grid. Additionally, in Figure 19, we visualize the solutions of the initial value problem with boundary conditions (34) for different solvers and different temporal and spatial discretizations at the last step t = 100. As an extension of Table 3, in Table 5, we compare our method to the numerical solvers at different discretizations. We also report the runtime of each method in Table 6. It can be observed that a larger number of edges (smaller dx) results in longer runtime, which is expected. Note that these values are only meant to convey a rough idea of time complexity since there may be several ways to improve the performance of the baselines. We ran all methods on the CPU.
0.1
Time
0.1
1.0
Time
1.0
0.32
Time
0.32
0.8
Time
0.8
Solution c(z, t)
dt1.0
Time
Time
dt1.0
0.6
Solution c(z, t)
0.6
3.12
Time
3.12
0.4
Time
0.4
10.0 Space
Space
Space
Space
Space
0.1
0.32
1.0
3.12
10.0
dz
0.0
(a) Analytical solution.
Time
10.0
0.2
Time
0.2
Space
Space
Space
Space
Space
0.1
0.32
1.0
3.12
10.0
dz
0.0
(b) MeGA-MP solution.
Figure 17: Solutions of the boundary value problem (33) for different temporal and spatial discretizations and on the whole resulting grid
46
0.1
Time
0.1
1.0
Time
1.0
0.32
Time
0.32
0.8
Time
0.8
Solution c(z, t)
dt1.0
Time
Time
dt1.0
0.6
Solution c(z, t)
0.6
3.12
Time
3.12
0.4
Time
0.4
10.0 Space
Space
Space
Space
Space
0.1
0.32
1.0
3.12
10.0
dz
Time
10.0
0.2
Time
0.2
0.0
Space
Space
Space
Space
Space
0.1
0.32
1.0
3.12
10.0
(a) Analytical solution.
dz
0.0
(b) MeGA-MP solution.
Figure 18: Solutions of initial value problem with boundary condition (34) for different temporal and spatial discretizations and on the whole resulting grid.
3.12
Time
10.0
Time
dt1.0
Analytical
Time
0.32
MeGA-MP
Time
0.1
RK4
Time
Semi-Lagrangian
Space 0.1
Space 0.32
Space 1.0
dz
Space 3.12
Space 10.0
Figure 19: Solutions of the initial value problem with boundary condition (33) for different solvers and different temporal and spatial discretizations at the last step t = 100.
47
Table 5: MAE between analytical ground truth and different solver solutions over all temporal and spatial grid points. The entry “–” indicates divergence of the method. dt 0.1
0.32
1.0
3.12
10.0
E
Proofs
E.1
Proof of Lemma 3.2
dx 0.10 0.32 1.00 3.12 10.00 0.10 0.32 1.00 3.12 10.00 0.10 0.32 1.00 3.12 10.00 0.10 0.32 1.00 3.12 10.00 0.10 0.32 1.00 3.12 10.00
semi- Lagr. 0.0016±0.0050 0.0046±0.0114 0.0108±0.0230 0.0239±0.0422 0.0562±0.0764 0.0027±0.0060 0.0038±0.0075 0.0095±0.0184 0.0227±0.0400 0.0557±0.0754 0.0080±0.0144 0.0080±0.0129 0.0097±0.0137 0.0206±0.0324 0.0545±0.0720 0.0250±0.0364 0.0243±0.0350 0.0231±0.0315 0.0254±0.0321 0.0504±0.0639 0.0763±0.0964 0.0752±0.0954 0.0728±0.0923 0.0675±0.0852 0.0671±0.0842
RK4 – 0.0013±0.0043 0.0035±0.0073 0.0105±0.0152 0.0288±0.0415 – 0.0108±0.0325 0.0042±0.0082 0.0117±0.0160 0.0304±0.0418 – – 0.0088±0.0155 0.0100±0.0155 0.0291±0.0395 – – – 0.0231±0.0333 0.0299±0.0406 – – – – 0.0728±0.0936
MeGA-MP 0.0024±0.0081 0.0013±0.0058 0.0012±0.0068 0.0023±0.0113 0.0057±0.0182 0.0071±0.0162 0.0049±0.0123 0.0031±0.0092 0.0036±0.0126 0.0071±0.0189 0.0157±0.0298 0.0143±0.0277 0.0103±0.0208 0.0081±0.0165 0.0115±0.0222 0.0384±0.0535 0.0375±0.0524 0.0345±0.0491 0.0242±0.0370 0.0238±0.0314 0.0939±0.0981 0.0933±0.0977 0.0915±0.0952 0.0843±0.0904 0.0761±0.0763
Lemma (Lemma 3.2). Let euv ∈ E hold. If the function cuv obeys Equation (2), for each t0 ∈ R and z0 ∈ R, the curve ´t γt0 ,z0 : Tt0 ,z0 −→ R≥0 , t 7−→ cuv t, z0 + t0 νuv (s) ds n o ´ t′ ′ is constant on Tt0 ,z0 = t ∈ R | 0 < z0 + t0 νuv (s) ds < luv ∀t ∈ min{t0 , t}, max{t0 , t} . Proof. Let cuv : R × Ωuv → R≥0 obey Equation (2) (note that Equation (2) is only defined for z ∈ Ωuv = (0, luv )). We can consider its space coordinate z as a function z : R → Ωuv parametrized by time. Consequently, we can use the method of characteristics (cf. Habermann (2013)) to derive the following system of ODEs from Equation (2) for the composed, now only time-dependent function cuv : R → R≥0 , cuv (t) := cuv (t, z(t)): dcuv (t) = 0, dt
dz(t) = νuv (t), dt
dt = 1. dt
By the fundamental theorem of calculus, we observe that for each t0 ∈ R and z0 ∈ R, the function ˆ t z : R −→ R, t 7−→ z(t) := z0 + νuv (s) ds t0
48
Table 6: Time in seconds that it takes each method to solve the BVP. The time it takes for the preprocessing to convert flows into transport times (s. Algorithm 1) is included in these values. The entry “–” indicates divergence of the method. dt 0.1
0.32
1.0
3.12
10.0
dx 0.10 0.32 1.00 3.12 10.00 0.10 0.32 1.00 3.12 10.00 0.10 0.32 1.00 3.12 10.00 0.10 0.32 1.00 3.12 10.00 0.10 0.32 1.00 3.12 10.00
semi- Lagr. 58.1295 16.9899 5.5392 1.7895 0.6871 18.0842 5.4553 1.6860 0.6071 0.2222 6.0633 1.7784 0.5320 0.1810 0.0699 1.7991 0.5470 0.1691 0.0582 0.0223 0.5620 0.1711 0.0570 0.0189 0.0075
RK4 – 1.0589 0.8999 0.8881 0.8863 – 0.3194 0.2795 0.2794 0.2766 – – 0.0893 0.0883 0.0871 – – – 0.0283 0.0281 – – – – 0.0088
MeGA-MP 63.1235 5.9495 1.3102 0.5884 0.1444 13.9027 1.7582 0.5185 0.2051 0.3668 6.0126 1.1763 0.2651 0.0628 0.0218 3.6538 0.7568 0.1539 0.0470 0.0191 2.3683 0.4693 0.1159 0.0400 0.0172
Analytical 0.2867 0.0792 0.0416 0.0097 0.0039 0.0398 0.0130 0.0049 0.0019 0.0013 0.0109 0.0029 0.0012 0.0007 0.0005 0.0024 0.0012 0.0007 0.0006 0.0004 0.0012 0.0008 0.0005 0.0005 0.0005
satisfies the second ODE. Therefore, by´the first ODE, we can conclude that the function cuv is constant along t characteristic curves (t, z(t)) = (t, z0 + t0 νuv (s) ds) along the time R as well as along the space interval Ωuv , i.e., within the edge euv of length luv > 0. Therefore, taking into account the domain and co-domain of cuv , ′ ′ νuv and z as a function, t0 , t ∈ R and z(t ) ∈ Ωuv = (0, luv ) need to hold for all t ∈ [min{t0 , t}, max{t0 , t}].17 ´t In summary, the curve γz0 : Tt0 ,z0 → R≥0 , t 7→ cuv (t, z0 + t0 νuv (s) ds) is constant on ˆ t n o ′ Tt0 ,z0 = t ∈ R | 0 < z0 + νuv (s) ds < luv ∀t ∈ min{t0 , t}, max{t0 , t} . ′
t0
E.2
Proof of Lemma C.3
Lemma (Lemma C.3). Let euv ∈ E, t0 ∈ R, z0 ∈ R hold and let Tt0 ,z0 be defined as in Lemma 3.2. 1. If z0 ̸∈ Ωuv holds, then Tt0 ,z0 = ∅ holds. 2. If z0 ∈ Ωuv holds, then Tt0 ,z0 ̸= ∅ holds. Even more, t0 ∈ Tt0 ,z0 holds. 3. If z0 ∈ Ωuv holds, there exist tl ∈ [−∞, t0 ) and tr ∈ (t0 , ∞] such that ´ 17 Especially to mention, the requirement z + t ν (s) ds ∈ Ω uv does not suffices, cf. Remark C.2. 0 t0 uv
49
´t (a) if tl ̸= −∞ holds, z0 + t0l νuv (s) ds ∈ {0, luv } holds, ´t (b) if tr ̸= +∞ holds, z0 + t0r νuv (s) ds ∈ {0, luv } holds and (c) Tt0 ,z0 = (tl , tr ) = (t0 − (t0 − tl ), t0 + (tr − t0 )) holds. Equivalently, there exist δtl , δtr ∈ (0, ∞] such that Tt0 ,z0 = (t0 − δtl , t0 + δtr ) holds. Proof. Similarly to Remark C.2, we define the function ˆ t νuv (s) ds.
z : R −→ R, t 7−→ z(t) := z0 + t0
By the fundamental theorem of calculus, the function z is differentiable and therefore, continuous. Moreover, using the function z, we can re-write the set Tt0 ,z0 as n o ′ ′ Tt0 ,z0 = t ∈ R | z(t ) ∈ Ωuv = (0, luv ) ∀t ∈ min{t0 , t}, max{t0 , t} . 1.: If z0 ̸∈ Ωuv holds, let us assume that Tt0 ,z0 ̸= ∅ holds. By definition of Tt0 ,z0 , there exists a t ∈ R, such ′ ′ ′ that z(t ) ∈ Ωuv holds for all t ∈ [min{t0 , t}, max{t0 , t}]. However, for t = t0 ∈ [min{t0 , t}, max{t0 , t}], we ´ ′ t obtain z(t ) = z(t0 ) = z0 + t00 νuv (s) ds = z0 ̸∈ Ωuv – a contradiction. 2.: If z0 ∈ Ωuv holds, we obtain z(t0 ) = z0 + t0 ∈ Tt0 ,z0 and Tt0 ,z0 ̸= ∅ holds.
´ t0 t0
νuv (s) ds = z0 ∈ Ωuv and therefore, trivially,
3.: We define the sets ′
′
′
′
Al :={t ∈ (−∞, t0 ] | z(t ) ̸∈ Ωuv } and Ar :={t ∈ [t0 , +∞) | z(t ) ̸∈ Ωuv }. Based on these sets, we choose ( max Al tl := −∞
if Al ̸= ∅ if Al = ∅
and
( min Ar tr := +∞
if Ar ̸= ∅ . if Ar = ∅
If z0 ∈ Ωuv holds, as seen above, z(t0 ) = z0 ∈ Ωuv holds. Therefore, we observe that t0 ̸∈ Al and t0 ̸∈ Ar , or equivalently, Al ⊂ (−∞, t0 ) and Ar ⊂ (t0 , +∞) hold. Consequently, our choices of tl and tr indeed satisfy tl ∈ [−∞, t0 ) and tr ∈ (t0 , ∞] and even more, t0 ∈ (tl , tr ) holds. We now want to show that for these choices, indeed, 5(a) - 5(c) hold. 3(a), 3(b): If tl ̸= −∞ holds, z0 + ´ tr z0 + t0 νuv (s) ds ∈ {0, luv } holds.
´ tl t0
νuv (s) ds ∈ {0, luv } holds and if tr ̸= +∞ holds,
The claims are equivalent to show that if Al ̸= ∅ holds, z(tl ) ∈ {0, luv } holds and if Ar ̸= ∅ holds, z(tr ) ∈ {0, luv } holds. If Al ̸= ∅ holds, let us assume that z(tl ) ̸∈ {0, luv } holds. Since tl = max Al ∈ Al , z(tl ) > luv or z(tl ) < 0 must hold. Without loss of generality (W.l.o.g.), let us assume that z(tl ) < 0 holds. Since z is continuous and z(t0 ) = z0 ∈ Ωuv , i.e., 0 < z(t0 ) holds, according to the intermediate value theorem, there exists a ′ ′ ′ t ∈ (tl , t0 ) such that z(t ) = 0 holds. Therefore, we can conclude that max Al = tl < t ∈ Al holds – a contradiction to the choice of tl = max Al . If Ar ̸= ∅ holds, let us assume z(tr ) ̸∈ {0, luv } holds. Analogous arguments lead to a contradiction. 3(c), step 1: Tt0 ,z0 ⊂ (tl , tr ) holds. 50
Let t ∈ Tt0 ,z0 ⊂ R and let us assume that t ̸∈ (tl , tr ), or equivalently, t ≤ tl < t0 or t ≥ tr > t0 holds. Case 1: If t ≤ t0 and Al ̸= ∅ holds, we observe that t ≤ tl = max Al < t0 must hold. By definition of ′ ′ Tt0 ,z0 , z(t ) ∈ Ωuv holds for all t ∈ [min{t0 , t}, max{t0 , t}] = [t, t0 ]. However, by definition of Al , for ′ ′ t = tl ∈ [t, t0 ) ⊂ [t, t0 ], z(t ) = z(tl ) ̸∈ Ωuv holds – a contradiction to the first conclusion. Case 2: If t ≤ t0 and Al = ∅ holds, we observe that t ≤ tl = −∞ must hold – a contradiction. Case 3: If t ≥ t0 and Al ̸= ∅ holds, we observe that t ≥ tr = min Ar > t0 must hold. By definition of ′ ′ Tt0 ,z0 , z(t ) ∈ Ωuv holds for all t ∈ [min{t0 , t}, max{t0 , t}] = [t0 , t]. However, by definition of Ar , for ′ ′ t = tr ∈ (t0 , t] ⊂ [t, t0 ], z(t ) = z(tr ) ̸∈ Ωuv holds – a contradiction to the first conclusion. Case 4: If t ≥ t0 and Al = ∅ holds, we observe that t ≥ tr = +∞ must hold – a contradiction. 3(c), step 2: (tl , tr ) ⊂ Tt0 ,z0 holds. Let t ∈ (tl , tr ) ⊂ R and let us assume that t ̸∈ Tt0 ,z0 holds. ′ Case 1: If t ≤ t0 and Al ̸= ∅ holds, by definition of Tt0 ,z0 , there exists a t ∈ [min{t0 , t}, max{t0 , t}] = [t, t0 ] ⊂ ′ ′ (tl , t0 ] = (max Al , t0 ] such that z(t ) ̸∈ Ωuv holds. Therefore, we can conclude that max Al = tl < t ∈ Al holds – a contradiction to the choice of tl = max Al . ′ Case 2: If t ≤ t0 and Al ̸= ∅ holds, by definition of Tt0 ,z0 , there exists a t ∈ [min{t0 , t}, max{t0 , t}] = ′ ′ [t, t0 ] ⊂ (tl , t0 ] = (−∞, t0 ] such that z(t ) ̸∈ Ωuv holds. Therefore, we can conclude that t ∈ Al holds – a contradiction to Al ̸= ∅. ′ Case 3: If t ≥ t0 and Ar ̸= ∅ holds, by definition of Tt0 ,z0 , there exists a t ∈ [min{t0 , t}, max{t0 , t}] = [t0 , t] ⊂ ′ ′ [t0 , tr ) = [t0 , min Ar ) such that z(t ) ̸∈ Ωuv holds. Therefore, we can conclude that min Ar = tr > t ∈ Ar holds – a contradiction to the choice of tr = min Al . ′ Case 4: If t ≥ t0 and Ar ̸= ∅ holds, by definition of Tt0 ,z0 , there exists a t ∈ [min{t0 , t}, max{t0 , t}] = ′ ′ [t0 , t] ⊂ [t0 , tr ) = [t0 , +∞) such that z(t ) ̸∈ Ωuv holds. Therefore, we can conclude that t ∈ Ar holds – a contradiction to Ar ̸= ∅. E.3
Proof of Theorem C.4
Theorem (Theorem C.4). Let euv ∈ E hold. If the function cuv obeys Equation (2), for each t0 ∈ R and z0 ∈ Ωuv , ˆ t
cuv t, z0 +
νuv (s) ds = cuv (t, z0 + ν uv (t0 , t) (t − t0 )) = cuv (t0 , z0 )
t0
holds for all t ∈ Tt0 ,z0 with Tt0 ,z0 as defined in Lemma 3.2 and for the mean velocity (ffl t ν uv (t0 , t) :=
ν (s) ds t0 uv 0
if t ̸= t0 = if t = t0
(
1 t−t0
0
´t t0
νuv (s) ds
if t ̸= t0 . if t = t0
´t Proof. For each t0 ∈ R and z0 ∈ Ωuv ⊂ R, by Lemma 3.2, γt0 ,z0 (t) = cuv (t, z0 + t0 νuv (s) ds) = c holds for all t ∈ Tt0 ,z0 and a constant c ∈ R≥0 . Thus, it suffices to show that c = cuv (t0 , z0 ) holds. Let us assume that c ̸= cuv (t0 , z0 ) holds. Since z0 ∈ Ωuv , by Lemma C.3, t0 ∈ Tt0 ,z0 holds. However, for ´t t = t0 ∈ Tt0 ,z0 , we obtain c = γt0 ,z0 (t) = cuv (t, z0 + t00 νuv (s) ds) = cuv (t0 , z0 ) – a contradiction. Therefore, ´t cuv (t, z0 + t0 νuv (s) ds) = cuv (t0 , z0 ) holds. Finally, by definition of ν uv (t0 , t), for t ̸= t0 , we obtain ν uv (t0 , t) · (t − t0 ) =
t − t0 t − t0
ˆ t
ˆ t νuv (s) ds =
t0
νuv (s) ds, t0
´t allowing to conclude that cuv (t, z0 + t0 νuv (s) ds) = cuv (t, z0 + ν uv (t0 , t) (t − t0 )) holds. 51
E.4
Proof of Lemma C.6
Lemma (Lemma C.6). Let euv ∈ E and t ∈ R hold. If there exists a pass-through time δtuv ∈ R w.r.t. euv ∈ E and t ∈ R, the following properties hold: 1. δtuv ∈ R̸=0 and thus, t − δtuv ̸= t holds. 2. There exists no other pass-through time of the same sign (i.e., there exists either exactly one or no positive pass-through time, and there exists exactly one or no negative pass-through time). 3. The mean velocity t
ν uv (t) := ν uv (t − δtuv , t) =
νuv (s) ds = t−δtuv
1 δtuv
ˆ t νuv (s) ds t−δtuv
during the time interval [min{t − δtuv , t}, max{t − δtuv , t}] as defined in Theorem C.4 is equal to and thus satisfies ν uv (t) = δtluv uv ν uv (t) > 0 ⇐⇒ δtuv > 0 and ν uv (t) < 0 ⇐⇒ δtuv < 0. q
satisfy 4. The flows quv (t − δtuv ) and quv (t) linked to the velocity νuv (·) = αvu(·) vu ( quv (t − δtuv ), quv (t) ≥ 0 if δtuv > 0 . quv (t − δtuv ), quv (t) ≤ 0 if δtuv < 0 Proof. 1.: Let us assume that δtuv = 0 holds. Consequently, (1) by basic results from analysis, (2) condition (19) and (3) definition of the pipe length luv > 0, ˆ t ˆ t (3) (1) (1) (2) 0 = νuv (s) ds = νuv (s) ds = luv > 0 t−0
t−δtuv
holds – a contradiction. 2.: Let us assume that there exist two different pass-through times δt1uv , δt2uv ∈ R̸=0 with the same sign. If δt1uv , δt2uv > 0 holds, without loss of generality (w.l.o.g.), let 0 < δt1uv < δt2uv hold. Consequently, ′ for t = t − δt1uv ∈ (t − δt2uv , t) = (min{t − δt2uv , t}, max{t − δt2uv , t}), by (1) by basic results from analysis and (2) condition (19), we obtain ˆ t′
(1)
ˆ t
νuv (s) ds = t−δt2uv
ˆ t νuv (s) ds −
t−δt2uv
(2)
(1)
νuv (s) ds = luv − luv = 0 t−δt1uv
– a contradiction to condition (20). ′
If δt1uv , δt2uv < 0 holds, w.l.o.g., let δt2uv < δt1uv < 0 hold. Consequently, for t = t − δt1uv ∈ (t, t − δt2uv ) = (min{t − δt2uv , t}, max{t − δt2uv , t}), we can use the same calculation to again obtain a contradiction. 3.: By (1) definition of the mean velocity ν uv (t − δtuv , t) as defined in Theorem C.4 and (2) condition (19), we indeed obtain ˆ t t (1) 1 (2) luv (1) ν uv (t) := ν uv (t − δtuv , t) = νuv (s) ds = νuv (s) ds = . δtuv t−δtuv δtuv t−δtuv The second claim then follows by the fact that by definition of the pipe length, luv > 0 holds.
52
4.: If δtuv > 0 holds, let us assume that quv (t − δtuv ) < 0, and consequently, νuv (t − δtuv ) < 0 holds. Since ′ νuv ∈ C 0 (R, R) is continuous, there exists a t ∈ (t − δtuv , t) = (min{t − δtuv , t}, max{t − δtuv , t}) such that ′ νuv (s) < 0 holds for all s ∈ (t − δtuv , t ). Consequently, by (1) basic results from analysis, we obtain ˆ t′
(1)
νuv (s) ds < 0 t−δtuv | {z } <0
– a contradiction to condition (20). Similarly, if δtuv > 0 holds, let us assume that quv (t) < 0, and consequently, νuv (t) < 0 holds. Since ′ νuv ∈ C 0 (R, R) is continuous, there exists a t ∈ (t − δtuv , t) = (min{t − δtuv , t}, max{t − δtuv , t}) such that ′ νuv (s) < 0 holds for all s ∈ (t , t). Consequently, by (1) basic results from analysis and (2) condition (19), we obtain ˆ t′
(1)
ˆ t
ˆ t
νuv (s) ds = t−δtuv
t−δtuv
|
νuv (s) t′ | {z } <0
νuv (s) ds − {z }
(2)
|
= luv
{z
(1)
ds > luv }
<0
{z
|
>0
}
– again a contradiction to condition (20). In contrast, if δtuv < 0 holds, let us assume that quv (t − δtuv ) > 0, and consequently, νuv (t − δtuv ) > 0 holds. ′ Since νuv ∈ C 0 (R, R) is continuous, there exists a t ∈ (t, t − δtuv ) = (min{t − δtuv , t}, max{t − δtuv , t}) such ′ that νuv (s) > 0 holds for all s ∈ (t , t − δtuv ). Consequently, by (1) basic results from analysis, we obtain ˆ t′
(1)
νuv (s) ds = −
ˆ t−δtuv t′
t−δtuv
|
νuv (s) ds < 0 | {z } >0 {z }
(1)
>0
– again a contradiction to condition (20). Similarly, if δtuv < 0 holds, let us assume that quv (t) > 0, and consequently, νuv (t) > 0 holds. Since ′ νuv ∈ C 0 (R, R) is continuous, there exists a t ∈ (t, t − δtuv ) = (min{t − δtuv , t}, max{t − δtuv , t}) such that ′ νuv (s) > 0 holds for all s ∈ (t, t ). Consequently, by (1) basic results from analysis and (2) condition (19), we obtain ˆ t′
(1)
ˆ t
ˆ t
νuv (s) ds = t−δtuv
νuv (s) ds − t−δtuv
(1)
ˆ t
=
t−δtuv
|
t′
νuv (s) ds
ˆ t′ νuv (s) ds + {z }
(2)
t
|
= luv
νuv (s) ds > luv | {z } >0 {z } (1)
>0
|
{z
>0
}
– again a contradiction to condition (20). Remark E.1. Note that if one of the flow rates discussed above are zero, that they can only be zero in this point of time, but not on a whole sub-interval of the time interval [min{t − δtuv , t}, max{t − δtuv , t}]. More precisely, if δtuv > 0 holds, let us assume that there exists an ϵ > 0 such that for all s ∈ [t − δtuv , t − ′ δtuv + ϵ), quv (s) = 0, and consequently, νuv (s) = 0, holds. Thus, for t = min{t − δtuv + ϵ, t} ∈ (t − δtuv , t) = 53
(min{t − δtuv , t}, max{t − δtuv , t}), by (1) basic results from analysis, we obtain ˆ t′
(1)
νuv (s) ds = 0 t−δtuv | {z } =0
– a contradiction to condition (20). In combination with the findings above, therefore, quv (t − δtuv ) ≥ 0 needs to hold, and if the (unlikely) case quv (t − δtuv ) = 0 indeed applies, we observe that for all ϵ > 0 there exists an s ∈ [t − δtuv , t − δtuv + ϵ) such that quv (s) ̸= 0, and consequently, νuv (s) ̸= 0, holds. Along the same lines, one can also show that indeed, quv (s) > 0, and consequently, νuv (s) > 0, holds (otherwise, we would again obtain a contradiction to condition (24)). The other cases are analogous adaptations of this discussion and the findings above.
E.5
Proof of Lemma C.7
Lemma (Lemma C.7). Let euv ∈ E and t ∈ R hold. If there exists a self-loop time δtuv ∈ R̸=0 w.r.t. euv ∈ E and t ∈ R, the following properties hold: 1. There exists no other self-loop time of the same sign (i.e., there exists either exactly one or no positive self-loop time, and there exists exactly one or no negative self-loop time). 2. The mean velocity 1 ν uv (t) := ν uv (t − δtuv , t) = νuv (s) ds = δtuv t−δtuv t
ˆ t νuv (s) ds t−δtuv
during the time interval [min{t − δtuv , t}, max{t − δtuv , t}] as defined in Theorem C.4 satisfies ν uv (t) = 0. 3. The flows quv (t − δtuv ) and quv (t) linked to the velocity νuv (·) = (
quv (t − δtuv ) ≥ 0, quv (t) ≤ 0 quv (t − δtuv ) ≤ 0, quv (t) ≥ 0
qvu(·) αvu
satisfy
if δtuv > 0 . if δtuv < 0
The proof of Lemma C.7 is very similar to the proof of Lemma C.6 and only requires a few changes. The interested reader can compare both proofs and investigate why the changes are required, and how they induce the different results in Lemma C.7.3 as compared to Lemma C.6.4. Proof. 1.: Let us assume that there exist two different self-loop times δt1uv , δt2uv ∈ R̸=0 with the same sign. ′ If δt1uv , δt2uv > 0 holds, w.l.o.g., let 0 < δt1uv < δt2uv hold. Consequently, for t = t − δt1uv ∈ (t − δt2uv , t) = (min{t − δt2uv , t}, max{t − δt2uv , t}), by (1) by basic results from analysis and (2) condition (23), we obtain ˆ t′
(1)
ˆ t
νuv (s) ds = t−δt2uv
ˆ t νuv (s) ds −
t−δt2uv
(2)
(1)
νuv (s) ds = 0 − 0 = 0 t−δt1uv
– a contradiction to condition (20). ′
If δt1uv , δt2uv < 0 holds, w.l.o.g., let δt2uv < δt1uv < 0 hold. Consequently, for t = t − δt1uv ∈ (t, t − δt2uv ) = (min{t − δt2uv , t}, max{t − δt2uv , t}), we can use the same calculation to again obtain a contradiction.
54
2.: By (1) definition of the mean velocity ν uv (t − δtuv , t) as defined in Theorem C.4 and (2) condition (23), we indeed obtain ˆ t t (1) (1) 1 (2) ν uv (t) := ν uv (t − δtuv , t) = νuv (s) ds = νuv (s) ds = 0. δtuv t−δtuv t−δtuv 3.: If δtuv > 0 holds, let us assume that quv (t − δtuv ) < 0, and consequently, νuv (t − δtuv ) < 0 holds. Since ′ νuv ∈ C 0 (R, R) is continuous, there exists a t ∈ (t − δtuv , t) = (min{t − δtuv , t}, max{t − δtuv , t}) such that ′ νuv (s) < 0 holds for all s ∈ (t − δtuv , t ). Consequently, by (1) basic results from analysis, we obtain ˆ t′
(1)
νuv (s) ds < 0 t−δtuv | {z } <0
– a contradiction to condition (24). Similarly, if δtuv > 0 holds, let us assume that quv (t) > 0, and consequently, νuv (t) > 0 holds. Since ′ νuv ∈ C 0 (R, R) is continuous, there exists a t ∈ (t − δtuv , t) = (min{t − δtuv , t}, max{t − δtuv , t}) such that ′ νuv (s) > 0 holds for all s ∈ (t , t). Consequently, by (1) basic results from analysis and (2) condition (23), we obtain ˆ t′ ˆ t ˆ t (1) νuv (s) ds = νuv (s) ds − νuv (s) ds < 0 t−δtuv t−δtuv t′ | {z } | {z } >0 | {z } (2) =0
(1)
>0
{z
|
<0
}
– again a contradiction to condition (24). In contrast, if δtuv < 0 holds, let us assume that quv (t − δtuv ) > 0, and consequently, νuv (t − δtuv ) > 0 holds. ′ Since νuv ∈ C 0 (R, R) is continuous, there exists a t ∈ (t, t − δtuv ) = (min{t − δtuv , t}, max{t − δtuv , t}) such ′ that νuv (s) > 0 holds for all s ∈ (t , t − δtuv ). Consequently, by (1) basic results from analysis, we obtain ˆ t′
(1)
νuv (s) ds = −
ˆ t−δtuv t′
t−δtuv
|
νuv (s) ds < 0 | {z } >0 {z }
(1)
>0
– again a contradiction to condition (24). Similarly, if δtuv < 0 holds, let us assume that quv (t) < 0, and consequently, νuv (t) < 0 holds. Since ′ νuv ∈ C 0 (R, R) is continuous, there exists a t ∈ (t, t − δtuv ) = (min{t − δtuv , t}, max{t − δtuv , t}) such that ′ νuv (s) < 0 holds for all s ∈ (t, t ). Consequently, by (1) basic results from analysis and (2) condition (23), we obtain ˆ t′ ˆ t ˆ t (1) νuv (s) ds = νuv (s) ds − νuv (s) ds t−δtuv
t′
t−δtuv
(1)
ˆ t′
ˆ t
=
t−δtuv
|
νuv (s) ds + {z }
(2)
t
|
=0
νuv (s) ds < 0 | {z } <0 {z } (1)
<0
| – again a contradiction to condition (24). 55
{z
<0
}
Remark E.2. Note that if one of the flow rates discussed above are zero, that they can only be zero in this point of time, but not on a whole sub-interval of the time interval [min{t − δtuv , t}, max{t − δtuv , t}]. More precisely, if δtuv > 0 holds, let us assume that there exists an ϵ > 0 such that for all s ∈ [t − δtuv , t − ′ δtuv + ϵ), quv (s) = 0, and consequently, νuv (s) = 0, holds. Thus, for t = min{t − δtuv + ϵ, t} ∈ (t − δtuv , t) = (min{t − δtuv , t}, max{t − δtuv , t}), by (1) basic results from analysis, we obtain ˆ t′
(1)
νuv (s) ds = 0 t−δtuv | {z } =0
– a contradiction to condition (24). In combination with the findings above, therefore, quv (t − δtuv ) ≥ 0 needs to hold, and if the (unlikely) case quv (t − δtuv ) = 0 indeed applies, we observe that for all ϵ > 0 there exists an s ∈ [t − δtuv , t − δtuv + ϵ) such that quv (s) ̸= 0, and consequently, νuv (s) ̸= 0, holds. Along the same lines, one can also show that indeed, quv (s) > 0, and consequently, νuv (s) > 0, holds (otherwise, we would again obtain a contradiction to condition (24)). The other cases are analogous adaptations of this discussion and the findings above.
E.6
Proof of Lemma C.8
Lemma (Lemma C.8). Let euv ∈ E and t ∈ R hold. δtuv ∈ R is an inverse pass-through time w.r.t. euv ∈ E and t ∈ R iff δtvu := δtuv ∈ R is a pass-through time w.r.t. evu ∈ E and t ∈ R Proof. “⇒” If δtuv ∈ R is an inverse pass-through time w.r.t. euv ∈ E and t ∈ R, by (1) the fact that the geometric pipe features are symmetric in the sense that lvu = luv and αvu = αuv holds for all v ∈ V and q u ∈ N (v), (2) definition of the flow velocity νvu (·) = αvu(·) , (3) the conservation of flows (cf. Eq. (3)) and vu (4) condition (21), we obtain ˆ t ˆ t (1,2,3) (4) (1) νvu (s) ds = − νuv (s) ds = −(−luv ) = lvu , t−δtuv
t−δtuv
that is, condition (19) is satisfied for δtvu := δtuv and w.r.t. evu ∈ E and t ∈ R. Similarly, by additionally using (5) condition (22), we obtain ˆ t′ νvu (s) ds
(1,2,3)
=
ˆ t
(5)
t−δtuv
(1)
νuv (s) ds ∈ (0, luv ) = (0, lvu ),
− t−δtuv
that is, condition (20) is satisfied for δtvu := δtuv and w.r.t. evu ∈ E and t ∈ R. Therefore, δtvu := δtuv ∈ R is a pass-through time w.r.t. evu ∈ E and t ∈ R. “⇐” If δtvu := δtuv ∈ R is a pass-through time w.r.t. evu ∈ E and t ∈ R, by (1) the fact that the geometric pipe features are symmetric in the sense that lvu = luv and αvu = αuv holds for all v ∈ V and q , (3) the conservation of flows (cf. Eq. (3)) and u ∈ N (v), (2) definition of the flow velocity νuv (·) = αuv(·) uv (4) condition (19), we obtain ˆ t ˆ t (1,2,3) (4) (1) νuv (s) ds = − νvu (s) ds = −(lvu ) = −luv , t−δtuv
t−δtuv
that is, condition (21) is satisfied for δtuv := δtvu and w.r.t. euv ∈ E and t ∈ R. Similarly, by additionally using (5) condition (20), we obtain ˆ t′ νuv (s) ds t−δtuv
(1,2,3)
=
ˆ t
(5)
(1)
νvu (s) ds ∈ (−lvu , 0) = (−luv , 0),
− t−δtuv
56
that is, condition (22) is satisfied for δtuv := δtvu and w.r.t. euv ∈ E and t ∈ R. Therefore, δtuv := δtvu ∈ R is an inverse pass-through time w.r.t. euv ∈ E and t ∈ R. E.7
Proof of Lemma C.9
Lemma (Lemma C.9). Let euv ∈ E and t ∈ R hold. δtuv ∈ R is an inverse self-loop time w.r.t. euv ∈ E and t ∈ R iff δtvu := δtuv ∈ R is a self-loop time w.r.t. evu ∈ E and t ∈ R The proof of Lemma C.9 is very similar to the proof of Lemma C.8 and only requires a few changes. The interested reader can compare both proofs and investigate why the changes are required, and how they induce the different results in Lemma C.9 as compared to Lemma C.9. Proof. “⇒” If δtuv ∈ R is an inverse self-loop time w.r.t. euv ∈ E and t ∈ R, by (1) the fact that the geometric pipe features are symmetric in the sense that lvu = luv and αvu = αuv holds for all v ∈ V and q , (3) the conservation of flows (cf. Eq. (3)) and u ∈ N (v), (2) definition of the flow velocity νvu (·) = αvu(·) vu (4) condition (25), we obtain ˆ t ˆ t (1,2,3) (4) νvu (s) ds = − νuv (s) ds = −0 = 0, t−δtuv
t−δtuv
that is, condition (23) is satisfied for δtvu := δtuv and w.r.t. evu ∈ E and t ∈ R. Similarly, by additionally using (5) condition (26), we obtain ˆ t′ νvu (s) ds
(1,2,3)
=
ˆ t
(5)
t−δtuv
(1)
νuv (s) ds ∈ (0, luv ) = (0, lvu ),
− t−δtuv
that is, condition (24) is satisfied for δtvu := δtuv and w.r.t. evu ∈ E and t ∈ R. Therefore, δtvu := δtuv ∈ R is a self-loop time w.r.t. evu ∈ E and t ∈ R. “⇐” If δtvu := δtuv ∈ R is a self-loop time w.r.t. evu ∈ E and t ∈ R, by (1) the fact that the geometric pipe features are symmetric in the sense that lvu = luv and αvu = αuv holds for all v ∈ V and q u ∈ N (v), (2) definition of the flow velocity νuv (·) = αuv(·) , (3) the conservation of flows (cf. Eq. (3)) and uv (4) condition (23), we obtain ˆ t ˆ t (1,2,3) (4) νuv (s) ds = − νvu (s) ds = −(0) = 0, t−δtuv
t−δtuv
that is, condition (25) is satisfied for δtuv := δtvu and w.r.t. euv ∈ E and t ∈ R. Similarly, by additionally using (5) condition (24), we obtain ˆ t′ νuv (s) ds t−δtuv
(1,2,3)
=
ˆ t
(5)
(1)
νvu (s) ds ∈ (−lvu , 0) = (−luv , 0),
− t−δtuv
that is, condition (26) is satisfied for δtuv := δtvu and w.r.t. euv ∈ E and t ∈ R. Therefore, δtuv := δtvu ∈ R is an inverse self-loop time w.r.t. euv ∈ E and t ∈ R. E.8
Proof of Theorem C.10
Theorem (Theorem C.10). Let euv ∈ E and t ∈ R hold. If the function cuv obeys Equation (2) and there exists a pass-through time δtuv ∈ R̸=0 w.r.t. euv ∈ E and t ∈ R, we obtain luv cuv (t, luv ) = cuv (t − δtuv , 0) = cuv t − ,0 . ν uv (t) 57
Even more, the concentrations cuv (·, luv ) and cvu (·, lvu ) at the end z = luv = lvu of the edges euv ∈ E and evu ∈ E are connected to the concentrations cu and cv of the nodes of that edge, respectively, by ( cuv (t, luv ) = cu (t − δtuv ) if δtuv > 0 . cvu (t − δtuv , lvu ) = cv (t) if δtuv < 0 Proof. The main idea of this proof is to use the continuity of the function cuv ∈ C 0 (R × [0, luv ], R≥0 ) to transfer the results from the advective transport along edges, i.e., Theorem C.4, to the ends of the edge euv ∈ E, i.e., to the functions cu , cv ∈ C 1 (R, R≥0 ). The crucial part is that we would like to apply Theorem C.4 for t0 = t − δtuv , z0 = 0 and t ∈ R as given. However, the first problem is that z0 = 0 ̸∈ Ωuv = (0, luv ) does not satisfy the condition of Theorem C.4. Moreover, since for this choice, ˆ t0 ˆ t z0 + νuv (s) ds = z0 + 0 = 0 and z0 + νuv (s) ds = z0 + luv = luv t0
t0
hold, t ̸∈ Tt0 ,z0 does also not satisfy the condition of Theorem C.4. Intuitively, we first have to transfer the information from the boundary of the edge euv to the inside of the edge using the continuity of the function cuv ∈ C 1 (R×[0, luv ], R≥0 ). Afterwards, we can use Theorem C.4 for slightly modified z0 and t. To do so, we define the two functions ˆ t−ϵ z :(min{0, δtuv }, max{0, δtuv }) −→ R, ϵ 7−→ z(ϵ) :=
νuv (s) ds and t−δtuv
z̃ :(min{0, δtuv }, max{0, δtuv }) −→ R, ϵ 7−→ z̃(e) :=
luv − z(ϵ) . 2
We first have to prove some technical details regarding these functions (Step 1 ) in order to be able to apply Theorem C.4 as illustrated above (Step 2 ). Afterwards, we can use the results from Lemma C.6 about the sign of the flow at time t − δtuv in dependence of the sign of the pass-through time δtuv in combination with the local inflow boundary conditions defined through Equation (5) to be able to transfer the results from the advective transport along edges to the end of the pipes (Step 3 ). Step 1: The functions z and z̃ satisfy the following properties: 1. z, z̃ and z + z̃ are continuous. 2. im(z) ⊂ (0, luv ), im(z̃) ⊂ (0, 12 luv ) ⊂ (0, luv ) and im(z + z̃) ⊂ ( 21 , luv ) ⊂ (0, luv ). 3. limϵ→0 z(ϵ) = luv , limϵ→0 z̃(ϵ) = 0 and limϵ→0 z(ϵ) + z̃(ϵ) = luv . Note that the limit is to be understood for ϵ ∈ (min{0, δtuv }, max{0, δtuv }) in the domain of z and z̃, i.e., if δtuv > 0, (min{0, δtuv }, max{0, δtuv }) = (0, δtuv ) holds, and the limit ϵ → 0 is to be understood as the limit ϵ ↘ 0, ϵ < δtuv . In contrast, if δtuv < 0, (min{0, δtuv }, max{0, δtuv }) = (δtuv , 0) holds, and the limit ϵ → 0 is to be understood as the limit ϵ ↗ 0, ϵ > δtuv . ´t 1.: By the fundamental theorem of calculus, the function t 7→ t−δtuv νuv (s) ds is differentiable and therefore, continuous. Therefore, as a composition of continuous functions, the functions z, z̃ and z + z̃ are continuous. 2.: By (1) basic transformations, (2) definition of z, z̃ and z + z̃, respectively, and (3) condition (20), we obtain (1)
im(z) = {z(ϵ) | ϵ ∈ (min{0, δtuv }, max{0, δtuv })} 58
(2)
ˆ t−ϵ
(1)
t−δtuv ˆ t−ϵ
=
=
νuv (s) ds
ϵ ∈ (min{0, δtuv }, max{0, δtuv })
νuv (s) ds
t − ϵ ∈ (min{t, t − δtuv }, max{t, t − δtuv })
t−δtuv (3)
⊂ (0, luv ).
By additionally using that therefore, (4) z(ϵ) ∈ (0, luv ) holds for all ϵ ∈ (min{0, δtuv }, max{0, δtuv }), we obtain (1)
im(z̃) = {z̃(ϵ) | ϵ ∈ (min{0, δtuv }, max{0, δtuv }) luv − z(ϵ) (2) = ϵ ∈ (min{0, δtuv }, max{0, δtuv }) 2 (4)
⊂ (0, 21 luv )
and
(1)
im(z + z̃) = {z(ϵ) + z̃(ϵ) | ϵ ∈ (min{0, δtuv }, max{0, δtuv })} luv − z(ϵ) (2) = z(ϵ) + ϵ ∈ (min{0, δtuv }, max{0, δtuv }) 2 luv + z(ϵ) (1) ϵ ∈ (min{0, δtuv }, max{0, δtuv }) = 2 (4)
⊂ ( 21 luv , luv ).
3.: By (1) continuity of z, z̃ and z + z̃, respectively, (2) definition of z, z̃ and z + z̃, respectively and (3) condition (19), we obtain ˆ t (1) (2) (3) lim z(ϵ) = z(0) = νuv (s) ds = luv , ϵ→0
t−δtuv
luv − z(0) (2,3) lim z̃(ϵ) = z̃(0) = = 0, ϵ→0 2 (1)
(2)
(1)
(2,3)
lim z(ϵ) + z̃(ϵ) = z(0) + z̃(0) =
ϵ→0
luv .
luv Step 2: cuv (t, luv ) = cuv (t − δtuv , 0) = cuv (t − ν uv (t) , 0) holds.
For each ϵ ∈ (min{0, δtuv }, max{0, δtuv }), we can now choose t0 = t − δtuv ∈ R and – as shown in Step 1 – z0 = z̃(ϵ) ∈ (0, luv ). Consequently, for each ϵ ∈ (min{0, δtuv }, max{0, δtuv }), by (1) the choices of t0 = t − δtuv ∈ R and z0 = z̃(ϵ) ∈ (0, luv ), (2) the definition of z and z̃ from Step 1 and (3) the results from Step 1, we observe that ˆ t−ϵ ˆ t−ϵ (3) (1) (2) z0 + νuv (s) ds = z̃(ϵ) + νuv (s) ds = z̃(ϵ) + z(ϵ) ∈ (0, luv ), t0
t−δtuv
and therefore, t − ϵ ∈ Tt0 ,z0 = Tt−δtuv ,z̃(ϵ) holds.18 Therefore, by (1) the results from Step 1, (2) the continuity of the function cuv ∈ C 1 (R × [0, luv ], R≥0 ), (3) definition of z and (4) Theorem C.4 applied to t0 = t−δtuv , z0 = z̃(ϵ) ∈ (0, luv ) and t−ϵ ∈ Tt0 ,z0 = Tt−δtuv ,z̃(ϵ) , we obtain (1) cuv (t, luv ) = cuv lim t − ϵ, lim z(ϵ) + z̃(ϵ) ϵ→0
ϵ→0
18 Note that all the technical work of constructing z and z̃ in Step 1 is because we would actually like to choose z
t ∈ R as given, however, this z0 and this t do not satisfy the requirements of Theorem C.4.
59
0 = 0 and
(2)
= lim cuv (t − ϵ, z(ϵ) + z̃(ϵ)) ϵ→0 ˆ t−ϵ (3) νuv (s) ds = lim cuv t − ϵ , z̃(ϵ) + ϵ→0 |{z} t−δtuv | {z } (1) ∈ (0,luv )
(1)
=z(ϵ) ∈ (0,luv )
{z
|
}
(1)
=z(ϵ)+z̃(ϵ) ∈ (0,luv ) (4)
= lim cuv (t − δtuv , z̃(ϵ)) ϵ→0 (2) = cuv t − δtuv , lim z̃(ϵ) ϵ→0
(1)
= cuv (t − δtuv , 0).
luv Finally, as seen in Lemma C.6.3, δtuv = ν uv (t) holds, and we can conclude
cuv (t − δtuv , 0) = cuv t −
luv ν uv (t)
.
Step 3: If δtuv > 0 holds, cuv (t, luv ) = cu (t − δtuv ) holds. If δtuv < 0 holds, cvu (t − δtuv , lvu ) = cv (t) holds. If δtuv > 0 holds, by (1) Step 2 and (2) Equation (5) together with the fact that by Lemma C.6.4, quv (t − δtuv ) ≥ 0 holds, we obtain (1)
(2)
cuv (t, luv ) = cuv (t − δtuv , 0) = cu (t − δtuv ). If instead, δtuv < 0 holds, by (1) Equation (18) for z = lvu , (2) Step 2 and (3) Equation (5) together with the fact that by the conservation of flows (Eq. (3)) and Lemma C.6.4, qvu (t) = −quv (t) ≥ 0 holds, we obtain (1)
(2)
(1)
(3)
cvu (t − δtuv , lvu ) = cuv (t − δtuv , 0) = cuv (t, luv ) = cvu (t, 0) = cv (t).
E.9
Proof of Theorem C.11
Theorem (Theorem C.11). Let euv ∈ E and t ∈ R hold. If the function cuv obeys Equation (2) and there exists a self-loop time δtuv ∈ R̸=0 w.r.t. euv ∈ E and t ∈ R, we obtain cuv (t, 0) = cuv (t − δtuv , 0). Even more, the concentration cvu (·, lvu ) at the end z = lvu = luv of the edge evu ∈ E is connected to the concentration cu of the node of that edge by ( cvu (t, lvu ) = cu (t − δtuv ) if δtuv > 0 . cvu (t − δtuv , lvu ) = cu (t) if δtuv < 0 The proof of Theorem C.11 is very similar to the proof of Theorem C.10 and only requires a few changes. The interested reader can compare both proofs and investigate why the changes are required, and how they induce the different results in Theorem C.11 as compared to Theorem C.10. Proof. The main idea of this proof is to use the continuity of the function cuv ∈ C 1 (R × [0, luv ], R≥0 ) to transfer the results from the advective transport along edges, i.e., Theorem C.4, to the end of the edge euv ∈ E, i.e., to the function cu ∈ C 1 (R, R≥0 ). 60
The crucial part is that we would like to apply Theorem C.4 for t0 = t − δtuv , z0 = 0 and t ∈ R as given. However, the first problem is that z0 = 0 ̸∈ Ωuv = (0, luv ) does not satisfy the condition of Theorem C.4. Moreover, since for this choice, ˆ t ˆ t0 νuv (s) ds = z0 + 0 = 0 νuv (s) ds = z0 + 0 = 0 and z0 + z0 + t0
t0
hold, t ̸∈ Tt0 ,z0 does also not satisfy the condition of Theorem C.4. Intuitively, we first have to transfer the information from the boundary of the edge euv to the inside of the edge using the continuity of the function cuv ∈ C 1 (R×[0, luv ], R≥0 ). Afterwards, we can use Theorem C.4 for slightly modified z0 and t. To do so, we define the two functions
ˆ t−ϵ
z :(min{0, δtuv }, max{0, δtuv }) −→ R, ϵ 7−→ z(ϵ) :=
νuv (s) ds and t−δtuv
z̃ :(min{0, δtuv }, max{0, δtuv }) −→ R, ϵ 7−→ z̃(e) :=
z(ϵ) . 2
We first have to prove some technical details regarding these functions (Step 1 ) in order to be able to apply Theorem C.4 as illustrated above (Step 2 ). Afterwards, we can use the results from Lemma C.7 about the sign of the flow at time t − δtuv in dependence of the sign of the self-loop time δtuv in combination with the local inflow boundary conditions defined through Equation (5) to be able to transfer the results from the advective transport along edges to the end of the pipe (Step 3 ). Step 1: The functions z and z̃ satisfy the following properties: 1. z, z̃ and z + z̃ are continuous. 2. im(z) ⊂ (0, luv ), im(z̃) ⊂ (0, 12 luv ) ⊂ (0, luv ) and im(z + z̃) ⊂ (0, 23 luv ) ⊂ (0, luv ). 3. limϵ→0 z(ϵ) = 0, limϵ→0 z̃(ϵ) = 0 and limϵ→0 z(ϵ) + z̃(ϵ) = 0. Note that the limit is to be understood for ϵ ∈ (min{0, δtuv }, max{0, δtuv }) in the domain of z and z̃, i.e., if δtuv > 0, (min{0, δtuv }, max{0, δtuv }) = (0, δtuv ) holds, and the limit ϵ → 0 is to be understood as the limit ϵ ↘ 0, ϵ < δtuv . In contrast, if δtuv < 0, (min{0, δtuv }, max{0, δtuv }) = (δtuv , 0) holds, and the limit ϵ → 0 is to be understood as the limit ϵ ↗ 0, ϵ > δtuv . ´t 1.: By the fundamental theorem of calculus, the function t 7→ t−δtuv νuv (s) ds differentiable and therefore, continuous. Therefore, as compositions of continuous functions, the functions z, z̃ and z + z̃ are continuous. 2.: By (1) basic transformations, (2) definition of z, z̃ and z + z̃, respectively, and (3) condition (24), we obtain (1)
im(z) = {z(ϵ) | ϵ ∈ (min{0, δtuv }, max{0, δtuv })} ˆ t−ϵ (2) = νuv (s) ds ϵ ∈ (min{0, δtuv }, max{0, δtuv }) t−δtuv
(1)
=
ˆ t−ϵ
νuv (s) ds
t − ϵ ∈ (min{t, t − δtuv }, max{t, t − δtuv })
t−δtuv (3)
⊂ (0, luv ).
By additionally using that therefore, (4) z(ϵ) ∈ (0, luv ) holds for all ϵ ∈ (min{0, δtuv }, max{0, δtuv }), we obtain (1)
im(z̃) = {z̃(ϵ) | ϵ ∈ (min{0, δtuv }, max{0, δtuv }) 61
(2)
=
ϵ ∈ (min{0, δtuv }, max{0, δtuv })
z(ϵ) 2
(4)
⊂ (0, 21 luv )
and
(1)
im(z + z̃) = {z(ϵ) + z̃(ϵ) | ϵ ∈ (min{0, δtuv }, max{0, δtuv })} z(ϵ) (2) ϵ ∈ (min{0, δtuv }, max{0, δtuv }) = z(ϵ) + 2 3 (1) = z(ϵ) ϵ ∈ (min{0, δtuv }, max{0, δtuv }) 2 (4)
⊂ (0, 32 luv ).
3.: By (1) continuity of z, z̃ and z + z̃, respectively, (2) definition of z, z̃ and z + z̃, respectively and (3) condition (23), we obtain ˆ t (1) (2) (3) lim z(ϵ) = z(0) = νuv (s) ds = 0, ϵ→0
t−δtuv
z(0) (2,3) lim z̃(ϵ) = z̃(0) = = 0, ϵ→0 2 (1)
(2)
(1)
(2,3)
lim z(ϵ) + z̃(ϵ) = z(0) + z̃(0) =
ϵ→0
0.
Step 2: cuv (t, 0) = cuv (t − δtuv , 0) holds. For each ϵ ∈ (min{0, δtuv }, max{0, δtuv }), we can now choose t0 = t − δtuv ∈ R and – as shown in Step 1 – z0 = z̃(ϵ) ∈ (0, luv ). Consequently, for each ϵ ∈ (min{0, δtuv }, max{0, δtuv }), by (1) the choices of t0 = t − δtuv ∈ R and z0 = z̃(ϵ) ∈ (0, luv ), (2) the definition of z and z̃ from Step 1 and (3) the results from Step 1, we observe that ˆ t−ϵ ˆ t−ϵ (3) (1) (2) z0 + νuv (s) ds = z̃(ϵ) + νuv (s) ds = z̃(ϵ) + z(ϵ) ∈ (0, luv ), t0
t−δtuv
and therefore, t − ϵ ∈ Tt0 ,z0 = Tt−δtuv ,z̃(ϵ) holds.19 Therefore, by (1) the results from Step 1, (2) the continuity of the function cuv ∈ C 1 (R × [0, luv ], R≥0 ), (3) definition of z and (4) Theorem C.4 applied to t0 = t−δtuv , z0 = z̃(ϵ) ∈ (0, luv ) and t−ϵ ∈ Tt0 ,z0 = Tt−δtuv ,z̃(ϵ) , we obtain (1) cuv (t, 0) = cuv lim t − ϵ, lim z(ϵ) + z̃(ϵ) ϵ→0
ϵ→0
(2)
= lim cuv (t − ϵ, z(ϵ) + z̃(ϵ)) ϵ→0 ˆ t−ϵ (3) = lim cuv t − ϵ , z̃(ϵ) + νuv (s) ds ϵ→0 |{z} t−δtuv {z } | (1) ∈ (0,luv )
(1)
=z(ϵ) ∈ (0,luv )
|
{z
(1)
}
=z(ϵ)+z̃(ϵ) ∈ (0,luv ) (4)
= lim cuv (t − δtuv , z̃(ϵ)) ϵ→0 (2) = cuv t − δtuv , lim z̃(ϵ) ϵ→0
19 Note that all the technical work of constructing z and z̃ in Step 1 is because we would actually like to choose z
t ∈ R as given; however, this z0 and this t do not satisfy the requirements of Theorem C.4.
62
0 = 0 and
(1)
= cuv (t − δtuv , 0).
Step 3: If δtuv > 0 holds, cvu (t, lvu ) = cu (t − δtuv ) holds. If δtuv < 0 holds, cvu (t − δtuv , lvu ) = cu (t) holds. If δtuv > 0 holds, by (1) Equation (18) for z = lvu , (2) Step 2 and (3) Equation (5) together with the fact that by Lemma C.7.4, quv (t − δtuv ) ≥ 0 holds, we obtain (1)
(2)
(3)
cvu (t, lvu ) = cuv (t, 0) = cuv (t − δtuv , 0) = cu (t − δtuv ). If instead, δtuv < 0 holds, by (1) Equation (18) for z = lvu , (2) Step 2 and (3) Equation (5) together with the fact that by Lemma C.7.4, quv (t) ≥ 0 holds, we obtain (1)
(2)
(3)
cvu (t − δtuv , lvu ) = cuv (t − δtuv , 0) = cuv (t, 0) = cu (t).
E.10
Proof of Theorem C.12
Theorem (Theorem C.12). Let euv ∈ E and t ∈ R hold. If νuv (t) ̸= 0 holds and the set ˆ t Auv (t) := δt ∈ (0, ∞) νuv (s) ds ∈ {−luv , 0, luv } t−δt = δt ∈ (0, ∞) ý(δt) ∈ {−luv , 0, luv } is not empty, the following properties hold: 1. δtuv := min Auv (t) = min{δt ∈ (0, ∞)
ý(δt) ∈ {−luv , 0, luv }} ∈ (0, ∞) exists.
2. δtuv is a positive pass-through time, inverse pass-through time, self-loop time or inverse self-loop time. More specifically, δtuv is (always w.r.t. euv and t) ... ... an inverse pass-through time
⇐⇒ ý(δtuv ) = −luv ,
... a self-loop time or an inverse self-loop time
⇐⇒ ý(δtuv ) = 0,
... a pass-through time
⇐⇒ ý(δtuv ) = luv .
3. There exists no other positive pass-through time, inverse pass-through time, self-loop time or inverse self-loop time than δtuv . ´ t′ ′ Proof. By the fundamental theorem of calculus, the function t 7→ t νuv (s) ds is differentiable and therefore, continuous. Therefore, as a composition of continuous functions, the function ý is continuous. 1: By assumption, Auv (t) ̸= ∅ holds. Moreover, by definition, the value 0 ∈ R is a lower bound of Auv (t) ⊂ (0, ∞). Therefore, it suffices to show that Auv (t) is closed in R; in this case, the minimum of Auv (t) exists. To do so, let (δtn )n∈N be a sequence in Auv (t) such that δtn → δt0 holds for n → ∞. We need to show that δt0 ∈ Auv (t) holds. Step 1.1: ý(δt0 ) ∈ {−luv , 0, luv } holds. On the one hand, by the continuity of ý, we obtain that ý(δtn ) → ý(δt0 ) holds for n → ∞. One the other hand, since (δtn )n∈N is a sequence in Auv (t), (ý(δtn ))n∈N is a sequence in {−luv , 0, luv }. Consequently, since {−luv , 0, luv } is closed in R, the limit ý(δt0 ) is also an element of {−luv , 0, luv }. 63
Step 1.2: δt0 > 0 holds. Since (δtn )n∈N is a sequence in Auv (t), δtn > 0 holds for a n ∈ N. Consequently, the limit δt0 = limn→∞ δtn satisfies δt0 ≥ 0, and it suffices to show that δt0 ̸= 0 holds. Let us assume that δt0 = 0 holds. By the continuity of νuv and the assumption that νuv (t) ̸= 0 holds, there exists an ϵ > 0 such that νuv (s) > 0 or νuv (s) < 0 holds for all s ∈ (t − ϵ, t + ϵ). W.l.o.g., we can assume that ϵ < luv holds (else, we choose a smaller ϵ that naturally satisfies this condition). Consequently, (1) by basic results from analysis, ˆ t (1) (36) νuv (s) ds ̸= 0 ý(δt) = t−δt | {z } >0 or <0
holds for all δt ∈ (0, ϵ). Specifically, since δtn → δt0 = 0 holds for n → ∞, for the ϵ > 0 from above, there exists an n0 ∈ N such that for all n ≥ n0 , |δtn | < ϵ holds; and since (δtn )n∈N is a sequence in Auv (t), even more, δtn ∈ (0, ϵ) holds for all n ≥ n0 . Consequently, for all n ≥ n0 , Equation (36) applies to δtn ; and since again (δtn )n∈N is a sequence in Auv (t), even more, ý(δtn ) ∈ {−luv , luv }
(37)
holds for all n ≥ n0 . At the same time, by the continuity of ý, we obtain that for the ϵ > 0 from above, there exists an n1 ∈ N such |ý(δtn ) − ý(δt0 )| < ϵ
(38)
holds for all n ≥ n1 . Bringing these findings together, by (1) basic transformations, (2) Equation (37) (3) the definition of ý (cf. Equation (7)), (4) the assumption that δt0 = 0 holds, (5) basic results from analysis, (6) Equation (38) and (7) the choice of ϵ, we obtain (1)
luv = | ± luv − 0|
(2,3,4,5)
=
(6)
(7)
|ý(δtn ) − ý(δt0 )| < ϵ < luv
for all n ≥ max{n0 , n1 } – a contradiction. 2: For δtuv := min Auv (t), we show an even stronger claim (namely Lemma C.14.1) in several steps. Step 2.1: im ý|(0,δtuv ) ⊂ (−luv , 0) ⊔ (0, luv ) holds. Let us assume that there exists a δt ∈ (0, δtuv ) such that ý(δt) ∈ (−∞, −luv ] ⊔ {0} ⊔ [luv , ∞) holds. Case 1: If ý(δt) ∈ {−luv , 0, luv } holds, we obtain δt ∈ Auv (t) with δt < δtuv – a contradiction to the choice of δtuv = min Auv (t) as the minimum of Auv (t). Case 2: If ý(δt) ∈ (−∞, −luv ) holds, by the fact that ý(0) = 0 holds and the intermediate value theorem, ′ ′ ′ there exists a δt ∈ (0, δt), such that ý(δt ) = −luv holds. Consequently, we obtain δt ∈ Auv (t) with ′ δt < δtuv – a contradiction to the choice of δtuv = min Auv (t) as the minimum of Auv (t). Case 3: If ý(δt) ∈ (luv , ∞) holds, by the fact that ý(0) = 0 holds and the intermediate value theorem, there ′ ′ ′ ′ exists a δt ∈ (0, δt), such that ý(δt ) = luv holds. Consequently, we obtain δt ∈ Auv (t) with δt < δtuv – a contradiction to the choice of δtuv = min Auv (t) as the minimum of Auv (t). Step 2.2: Either im z|(0,δtuv ) ⊂ (−luv , 0) or im z|(0,δtuv ) ⊂ (0, luv ) holds. By step Step 2.1, im z|(0,δtuv ) ⊂ (−luv , 0) ⊔ (0, luv ) holds. Let us assume that there exist δt1 , δt2 ∈ (0, δtuv ) such that ý(δt1 ) ∈ (−luv , 0) and ý(δt2 ) ∈ (0, luv ) holds. Again by the mean value theorem, there exists a 64
δt ∈ (δt1 , δ2 ), such that ý(δt) = 0 holds. Consequently, we obtain δt ∈ Auv (t) with δt < δt2 < δtuv – a contradiction to the choice of δtuv = min Auv (t) as the minimum of Auv (t). Step 2.3: If additionally, ý(δtuv ) = −luv holds, then im z|(0,δtuv ) ⊂ (−luv , 0) holds. If in contrast, additionally, ý(δtuv ) = luv holds, then im z|(0,δtuv ) ⊂ (0, luv ) holds. By step Step 2.2, either im z|(0,δtuv ) ⊂ (−luv , 0) or im z|(0,δtuv ) ⊂ (0, luv ) holds. If additionally, ý(δtuv ) = −luv holds, let us assume that there exists a δt ∈ (0, δtuv ) such that ý(δt) ∈ (0, luv ) holds. Again by the ′ ′ intermediate value theorem, there exists a δt ∈ (δt, δtuv ), such that ý(δt ) = 0 holds. Consequently, we ′ ′ obtain δt ∈ Auv (t) with δt < δtuv – a contradiction to the choice of δtuv = min Auv (t) as the minimum of Auv (t). If in contrast, ý(δtuv ) = luv holds, let us assume that there exists a δt ∈ (0, δtuv ) such that ý(δt) ∈ (−luv , 0) ′ ′ holds. Again by the intermediate value theorem, there exists a δt ∈ (δt, δtuv ), such that ý(δt ) = 0 holds. ′ ′ Consequently, we obtain δt ∈ Auv (t) with δt < δtuv – a contradiction to the choice of δtuv = min Auv (t) as the minimum of Auv (t). Step 2.4: For δtuv := min Auv (t), we obtain im z|(0,δtuv ) ⊂ (−luv , 0) im z|(0,δtuv ) ⊂ (−luv , 0) or im z|(0,δtuv ) ⊂ (0, luv ) im z|(0,δtuv ) ⊂ (0, luv )
if ý(δtuv ) = −luv , if ý(δtuv ) = 0, if ý(δtuv ) = luv .
The claim follows immediately from Step 2.2 and Step 2.3. Step 2.5: δtuv is (always w.r.t. euv and t) ... ... an inverse pass-through time
⇐⇒ ý(δtuv ) = −luv ,
... a self-loop time or an inverse self-loop time
⇐⇒ ý(δtuv ) = 0,
... a pass-through time
⇐⇒ ý(δtuv ) = luv .
“⇒” If δtuv is an inverse pass-through time, by (1) definition of ý and (2) condition (21), we obtain ˆ t (1) (2) ý(δtuv ) = νuv (s) ds = −luv . t−δtuv
Analogously, if δtuv is a self-loop time or an inverse self-loop time, we use condition (23) or (25) to conclude that ý(δtuv ) = 0 holds. If δtuv is a pass-through time, we use condition (19) to conclude that ý(δtuv ) = luv holds. “⇐” If ý(δtuv ) = −luv holds, by (1) definition of ý and (2) this assumption, we obtain ˆ t (1) (2) νuv (s) ds = ý(δtuv ) = −luv . t−δtuv
that is, condition (21) is satisfied. ′
′
Consequently, for any t ∈ (min{t−δtuv , t}, max{t−δtuv , t}) = (t−δtuv , t), or equivalently, for any t = t−δt with δt ∈ (0, δtuv ), by (1) by basic results from analysis, (2) condition (21), (3) definition of ý and (4) Step 2.4, we obtain ˆ t′
(1)
νuv (s) ds = t−δtuv
ˆ t (1) νuv (s) ds ∈ (−luv , 0), νuv (s) ds − t−δt t−δtuv {z } {z } | |
ˆ t
(2)
(3)
= −luv
(4)
= ý(δt) ∈ (−luv ,0)
65
that is, condition (22) is satisfied, too, and thus, δtuv is an inverse pass-through time. Analogously, if ý(δtuv ) = 0 holds, we use this assumption to show that condition (23), which is equal to condition (25), is satisfied. If ý(δtuv ) = luv holds, we use this assumption to show that condition (19) is satisfied. Consequently, we use these conditions and Step 2.4 to show that condition (24), (26) and (20) are satisfied in alignment with the conditions in Step 2.4. We leave the details as an exercise to the reader, however, a summary of the cases to be considered is given as follows: ˆ t
ˆ t t−δtuv
|
= = = =
(1)
ˆ t′ =t−δt
ˆ t′ =t−δt
νuv (s) ds = νuv (s) ds − νuv (s) ds t−δt t−δtuv {z } {z } | | {z } −luv 0 0 luv
∈ ∈ ∈ ∈
∈ ∈ ∈ ∈
(−luv , 0) (−luv , 0) (0, luv ) (0, luv )
(−luv , 0) (0, luv ) (−luv , 0) (0, luv )
t−δtuv
|
νuv (s) ds {z }
if ý(δtuv ) = −luv if ý(δtuv ) = 0 and im z|(0,δtuv ) ⊂ (−luv , 0) if ý(δtuv ) = 0 and im z|(0,δtuv ) ⊂ (0, luv ) if ý(δtuv ) = luv
3.: Similar to the proof of Lemma C.6.2 and Lemma C.7.1, let us assume that there exist two different and positive times δt1uv , δt2uv ∈ R̸=0 which are each either a pass-through time, an inverse pass-through time, a self-loop time or an inverse self-loop time. W.l.o.g., let 0 < δt1uv < δt2uv hold. Consequently, for ′ t = t − δt1uv ∈ (t − δt2uv , t) = (min{t − δt2uv , t}, max{t − δt2uv , t}), by (1) by basic results from analysis and (2) one or two of the conditions (19), (21), (23) or (25)20 , we obtain ˆ t′
(1)
ˆ t
ˆ t
νuv (s) ds = t−δt2uv
(2)
νuv (s) ds − t−δt2uv
νuv (s) ds ∈ {−2luv , −luv , 0, luv , 2luv } t−δt1uv
– a contradiction to one of the conditions (20), (22), (24) or (26)21 . E.11
Proof of Lemma C.14
Lemma (Lemma C.14). Let euv ∈ E and t ∈ R hold. In the setting of Theorem C.12, the following property holds: 1. δtuv is (always w.r.t. euv and t) ... ... a pass-through time
⇐⇒ ý(δtuv ) ∈ {−luv , luv } and quv (t) > 0,
... an inverse pass-through time
⇐⇒ ý(δtuv ) ∈ {−luv , luv } and quv (t) < 0,
... a self-loop time
⇐⇒ ý(δtuv ) = 0 and quv (t) < 0,
... an inverse self-loop time
⇐⇒ ý(δtuv ) = 0 and quv (t) > 0.
Proof. “⇒” If δtuv is a pass-through time, similarly to Step 2.5 in the proof of Theorem C.12.2, by (1) definition of ý and (2) condition (19), we obtain (1)
ˆ t
ý(δtuv ) =
(2)
νuv (s) ds = luv ∈ {−luv , luv }. t−δtuv
Moreover, by Lemma C.6 and the choice of δtuv > 0 , quv (t) ≥ 0 holds. Since quv (t) ̸= 0 holds by assumption, quv (t) > 0 follows. 20 Note that we have 3 · 3 = 9 possible combinations: If δt 2uv is a pass-through time, we work with condition (19). If δt2uv is an inverse pass-through time, we work with condition (21). If δt2uv is a(n) (inverse) self-loop time, we work with condition (23), which is equal to condition (25). Each of these three options need to be combined with the three analogous options for the second time δt1uv . 21 Note that we have four possible combinations: If δt 2uv is a pass-through time, we work with condition (20). If δt2uv is an inverse pass-through time, we work with condition (22). If δt2uv is a self-loop time, we work with condition (24). If δt2uv is an inverse self-loop time, we work with condition (26).
66
Similarly, if δtuv is an inverse pass-through time (w.r.t. euv and t), by (1) definition of ý and (2) condition (21), we obtain ˆ t (1) (2) ý(δtuv ) = νuv (s) ds = −luv ∈ {−luv , luv }. t−δtuv
Additionally, we use Lemma C.8 to conclude that δtuv is a pass-through time w.r.t. evu and t. Therefore, by (1) the conservation of flows (Eq. (3)), (2) Lemma C.6 and the choice of δtuv > 0 considered as a pass-through time w.r.t. evu and t, we obtain (1)
quv (t) = − qvu (t) ≤ 0. | {z } (2)
≥0
Since quv (t) ̸= 0 holds by assumption, quv (t) < 0 follows. Analogously, if δtuv is a self-loop time, we use condition (23) to conclude that ý(δtuv ) = 0 holds and Lemma C.7 to conclude that quv (t) ≤ 0 and thus, quv (t) < 0 holds. If δtuv is an inverse self-loop time, we use condition (25) to conclude that ý(δtuv ) = 0 holds and Lemma C.9, the conservation of flows and Lemma C.7 to conclude that quv (t) = −qvu (t) ≥ 0 and thus, quv (t) > 0 holds. “⇐” If ý(δtuv ) ∈ {−luv , luv } holds, we know by Theorem C.12.2 that δtuv is either a pass-through time or an inverse pass-through time. If additionally, quv (t) > 0 holds, let us assume that δtuv is an inverse pass-through time. Consequently, by the same arguments as in the inclusion above, we can conclude that in this case, quv (t) = −qvu (t) ≤ 0 holds – a contradiction. If instead, additionally, quv (t) < 0 holds, let us assume that δtuv is a pass-through time. Consequently, by the same arguments as in the inclusion above, we can conclude that in this case, quv (t) ≥ 0 holds – a contradiction. If ý(δtuv ) = 0 holds, we know by Theorem C.12.2 that δtuv is either a self-loop time or an inverse self-loop time. If additionally, quv (t) < 0 holds, let us assume that δtuv an inverse self-loop time. Consequently, by the same arguments as in the inclusion above, we can conclude that in this case, quv (t) ≥ 0 holds – a contradiction. If instead, additionally, quv (t) > 0 holds, let us assume that δtuv is a self-loop time. Consequently, by the same arguments as in the inclusion above, we can conclude that in this case, quv (t) ≤ 0 holds – a contradiction. E.12
Proof of Theorem 3.5 and C.15
Theorem (Theorem 3.5 and C.15). Let v ∈ V and t ∈ R hold. We assume that Assumption 3.4 holds and that for each u ∈ N (v), Auv (t) is not empty (that is, the OP (8) has a solution). If the function cv obeys the mixing at nodes (Eq. (4)) and if for each u ∈ N (v), the function cuv obeys the advection edge dynamics (Eq. (2)) with inflow continuity condition (Eq. (5)), then X X X cv (t) = wuv (t) · cu (t − δtuv ) + wuv (t) · cv (t − δtuv ) = wuv (t) · φ(t, euv , c, ν) u∈N− (v,t) ý(δtuv )=±luv
u∈N− (v,t) ý(δtuv )=0
u∈N (v)
holds, where δtuv = min Auv (t) is the transport time as defined in Theorem C.12. In other words, Equation (10) with ϕ as defined in Equation (9) holds. Proof. Well-definedness: Since by assumption, for each u ∈ N (v), the set Auv (t) as defined in Theorem C.12 is not empty, according to Theorem C.12.1, a transport time δtuv = min Auv (t) ∈ (0, ∞) exists for each 67
u ∈ N (v). Even more, by Theorem C.12.3, is it uniquely determined by the condition that it is positive, allowing to define the map N (v) → (0, ∞), u 7→ δtuv that associates a unique transport time with each neighbor u ∈ N (v). Especially to mention, for each u ∈ N (v) and by definition of the transport time δtuv (Theorem C.12), ý(δtuv ) ∈ {−luv , 0, luv } holds. Equality: By (1) definition of the in-, out- and no-flow neighborhoods and (2) the conservation of flows (Eq. (3)), we obtain (1)
(2)
(1)
(2)
N± (v, t) : = {u ∈ N (v) | sgn(qvu (t)) = ±1} = {u ∈ N (v) | sgn(quv (t)) = ∓1} and N0 (v) : = {u ∈ N (v) | sgn(qvu (t)) = ±0} = {u ∈ N (v) | sgn(quv )(t) = ±0}. Clearly, these three sub-neighborhoods define a partition of the neighborhood N (v), that is, they are pairwise disjoint and N (v) = N− (v, t) ⊔ N+ (v, t) ⊔ N0 (v, t) holds. Together with the definition of the ReLU function ReLU(q) = max{0, q} for all q ∈ R, we observe that for each u ∈ N (v), quv (t) if u ∈ N− (v, t) ReLU(quv (t)) = 0 (39) if u ∈ N+ (v, t) 0 if u ∈ N0 (v, t) holds. Consequently, by (1) Equation (4), (2) Equation (39), (3) the definition of the weights wuv (t) (Eq. (6)), (4) the well-definedness investigations from above and (5) basic transformations, we obtain P (1) u∈N− (v,t) quv (t) · cuv (t, luv ) P cv (t) = u∈N− (v,t) quv (t) P (2) u∈N (v) ReLU(quv (t)) · cuv (t, luv ) P = u∈N (v) ReLU(quv (t)) X (3) = wuv (t) · cuv (t, luv ) u∈N (v) (3)
=
X
wuv (t) · cuv (t, luv )
u∈N− (v,t) (4,5)
=
X
wuv (t) · cuv (t, luv ) +
X
wuv (t) · cuv (t, luv ).
u∈N− (v,t) ý(δtuv )=0
u∈N− (v,t) ý(δtuv )=±luv
This proves the first equality. Since the definition of N− (v, t) indicates that quv (t) > 0 is positive, by Lemma C.14.1, the transport times δtuv > 0 in the first sum are pass-through times, while the transport times δtuv > 0 in the second sum are inverse self-loop times. Since again by Lemma C.14.1, the other direction also holds, we can conclude the calculations as follows: By (1) applying Theorem C.10 to cuv (t, luv ) in the first sum and Lemma C.8 and Theorem C.11 to cuv (t, luv ) the second sum, (2) the definition of ϕ (Equation (9)), (3) basic transformations, (4) the definition of the transport time δtuv (Theorem C.12) and again, (5) the definition of the weights wuv (t) (Eq. (6)), we obtain X X (1) cv (t) = wuv (t) · cu (t − δtuv ) + wuv (t) · cv (t − δtuv ) u∈N− (v,t) ý(δtuv )=±luv (2)
=
X
u∈N− (v,t) ý(δtuv )=0
wuv (t) · ϕ(t, euv , c, ν) +
u∈N− (v,t) ý(δtuv )=±luv (3)
=
X
X u∈N− (v,t) ý(δtuv )=0
wuv (t) · ϕ(t, euv , c, ν)
u∈N− (v,t) ý(δtuv )∈{−luv ,0,luv }
68
wuv (t) · ϕ(t, euv , c, ν)
(4)
=
X
wuv (t) · ϕ(t, euv , c, ν)
u∈N− (v,t) (5)
=
X
wuv (t) · ϕ(t, euv , c, ν).
u∈N (v)
This proves the second equality. E.13
Proof of Theorem C.16
Theorem (Theorem C.16). Let T = (A, E) be a finite tree with root node a0 ∈ A, internal nodes Aint and leaf nodes Aleaf (that is, A = {a0 } ⊔ Aint ⊔ Aleaf holds). If each node a ∈ A is associated with a magnitude xa ∈ R that satisfies X xa = ωab · xb for all a ∈ {a0 } ⊔ Aint , (40) b∈C(a)
where the set C(a) ⊂ A consists of the children of a ∈ A and wab is a weight of the edge eab ∈ E that depends on both the parent a ∈ A and its child b ∈ C(a), then the magnitude x0 := xa0 of the root node a0 ∈ A is determined by the magnitudes xlp := xalp of the leaf nodes alp ∈ Aleaf by X x0 = ω(a0 , ..., alp ) · xlp , (41) p=(a0 ,...,alp )∈P(a0 ,Aleaf )
where P(a0 , Aleaf ) be the set of paths connecting the root node a0 to one of its leaf nodes Aleaf and lp −1
ω(a0 , ..., alp ) :=
Y
(42)
ωal al+1
l=0
is the product of weights along a(n) (arbitrary) path p = (a0 , ..., alp ) ∈ P with length lp ∈ N. Proof. For the sake of brevity, in this proof, we use the abbreviation P = P(a0 , Aleaf ). Step 1: For each l ∈ N, we obtain X xa0 =
X
ω(a0 , ..., alp ) · xalp +
p=(a0 ,...,alp )∈P: lp ≤l
ω(a0 , ..., al ) · xal .
(43)
p=(a0 ,...,alp )∈P: lp >l
We prove the claim by induction to l ∈ N. Induction base: If l = 1, by (1) Equation (40) applied to xa0 , (2) a change in notation, (3) basic transformations, (4) the definition of ω (cf. (42)) and (5) the choice of l = 1, we obtain X (1) xa0 = ωa0 b · xb b∈C(a0 ) (2,3)
=
(3)
=
X a1 ∈C(a0 )∩Aleaf
a1 ∈C(a0 )∩Aint
X
X
p=(a0 ,a1 )∈P: lp =1 (4)
=
X
ωa0 a1 · xa1 +
X p=(a0 ,a1 )∈P: lp =1
ωa0 a1 · xa1 +
ωa0 a1 · xa1 wa0 a1 · xa1
p=(a0 ,a1 ,...,alp )∈P: lp >1
X
ω(a0 , a1 ) · xa1 +
p=(a0 ,...,alp )∈P: lp >1
69
ω(a0 , a1 ) · xal
(3,5)
X
=
X
ω(a0 , ..., alp ) · xalp +
p=(a0 ,...,alp )∈P: lp ≤l
ω(a0 , ..., al ) · xal .
p=(a0 ,...,alp )∈P: lp >l
Induction hypothesis: Equation (43) holds for all ˜l ≤ l ∈ N. Induction step: By (1) the induction hypothesis, we obtain X (1) xa0 = ω(a0 , ..., alp ) · xalp + p=(a0 ,...,alp )∈P: lp ≤l
X
ω(a0 , ..., al ) · xal .
p=(a0 ,...,alp )∈P: lp >l
Similar to the induction base, by (1) Equation (40) applied to xal , (note that since l < lp , al ∈ Aint must hold), (2) a change in notation, (3) basic transformations and (4) the definition of ω (cf. (42)), we transform the second summand to X ω(a0 , ..., al ) · xal p=(a0 ,...,alp )∈P: lp >l
(1)
X
=
X
ω(a0 , ..., al ) ·
p=(a0 ,...,alp )∈P: lp >l
ωal b · xb
b∈C(al )
(2,3)
X
=
ω(a0 , ..., al ) ·
X
X
p=(a0 ,...,alp )∈P: lp ≥l+1
al+1 ∈C(al )∩Aleaf
X
X
p=(a0 ,...,alp )∈P: lp ≥l+1
al+1 ∈C(al )∩Ain
=
+
(3)
X
=
ωal al+1 · xal+1 +
al+1 ∈C(al )∩Aleaf
p=(a0 ,...,alp )∈P: lp ≥l+1 (3)
X
X
ωal al+1 · xal+1
al+1 ∈C(al )∩Ain
ω(a0 , ..., al ) · ωal al+1 · xal+1
ω(a0 , ..., al ) · ωal al+1 · xal+1
ω(a0 , ..., al ) · ωal al+1 · xal+1
p=(a0 ,...,alp )∈P: lp =l+1
X
+
ω(a0 , ..., al ) · ωal al+1 · xal+1
p=(a0 ,...,alp )∈P: lp >l+1 (4)
X
=
X
ω(a0 , ..., al , al+1 ) · xal+1 +
p=(a0 ,...,alp )∈P: lp =l+1
ω(a0 , ..., al , al+1 ) · xal+1 .
p=(a0 ,...,alp )∈P: lp >l+1
Bringing it all together, again by (1) basic transformations, we obtain X X xa0 = ω(a0 , ..., alp ) · xalp + ω(a0 , ..., al ) · xal p=(a0 ,...,alp )∈P: lp ≤l
=
X
p=(a0 ,...,alp )∈P: lp >l
ω(a0 , ..., alp ) · xalp
p=(a0 ,...,alp )∈P: lp ≤l
+
X
ω(a0 , ..., al , al+1 ) · xal+1
p=(a0 ,...,alp )∈P: lp =l+1
70
X
+
ω(a0 , ..., al , al+1 ) · xal+1
p=(a0 ,...,alp )∈P: lp >l+1 (1)
X
=
ω(a0 , ..., alp ) · xalp
p=(a0 ,...,alp )∈P: lp <l+1
X
+
ω(a0 , ..., alp ) · xalp
p=(a0 ,...,alp )∈P: lp =l+1
X
+
ω(a0 , ..., al+1 ) · xal+1
p=(a0 ,...,alp )∈P: lp >l+1 (1)
X
=
ω(a0 , ..., alp ) · xalp
+
X
ω(a0 , ..., al+1 ) · xal+1 ,
p=(a0 ,...,alp )∈P: lp >l+1
p=(a0 ,...,alp )∈P: lp ≤l+1
which concludes the induction step. Step 2: Equation (41) holds. Since T is finite, that is, T an acyclic graph with a finite number of nodes |A| < ∞, also each path p ∈ P has to be finite. Therefore, we can choose l = maxp∈P lp ∈ N. Consequently, the second sum in Step 1 is empty and the condition in the first sum in Step 1 is naturally satisfied. This yields Equation (41).
E.14
Proof of Theorem C.19
Theorem (Theorem C.19). Let v ∈ V and t ∈ R hold. Let T (v, t) be the inflow tree of v at time t as defined in Definition C.17. If cv (t) is defined by Equation (10), then cv (t) =
lp −1
Y
wul+1 vl (tl ) · cvlp (tlp )
X p=((v0 ,t0 )=(v,t),(u1 ,v1 ,t1 ),...,(ulp ,vlp ,tlp ))∈P(v,t)
l=0
with weights wul+1 vl (tl ) according to Equation (6) holds.
Proof. For each layer l ∈ N0 of the tree, we now associate each node al = (ul , vl , tl ) ∈ A of that layer with the magnitude xl := xal := cvl (tl ). Step 1: Together with these magnitudes, T (v, t) defines a tree in terms of Theorem C.16. In order to prove the claim, we first of all investigate on the term ϕ(tl , eul+1 vl , c, ν) for a given parent node al = (ul , vl , tl ) in the l-th layer and its child node al+1 = (ul+1 , vl+1 , tl+1 ) ∈ C(al ).22 Indeed, by (1) Equation (9) and (2) the definition of the inflow tree T (v, t) (Definition C.17), for any l ∈ N0 , ϕ(tl , eul+1 vl , c, ν) is
22 Especially to mention, by definition of the inflow tree T (v, t) (Definition C.17), u l+1 ∈ N− (vt , tl ) holds, that is, the edge eul+1 vl ∈ E exists.
71
equal to ϕ(tl , eul+1 vl , c, ν) ( cul+1 (tl − ù(tl , eul+1 vl , ν)) if ý(ù(tl , eul+1 vl , ν)) ∈ {−luv , luv }, (1) = if ý(ù(tl , eul+1 vl , ν)) = 0 +1 cvl (tl − ù(tl , eul+1 vl , ν)) ( cul+1 (tl+1 ) if ý(ù(tl , eul+1 vl , ν)) ∈ {−luv , luv }, (2) = c (tl+1 ) if ý(ù(tl , eul+1 vl , ν)) = 0 +1 vl
(44)
(2)
= cvl+1 (tl+1 )
(2)
= xl+1 .
Consequently, again by (1) the definition of the inflow tree T (v, t) (Definition C.17), (2) Equation (10) and (3) Equation (44), for any al ∈ {a0 } ⊔ Aint , we obtain (1)
xal := xl = cvl (tl ) X (2) = wul+1 vl (tl ) · ϕ(tl , eul+1 vl , c, ν) ul+1 ∈N− (vl ,tl ) (3)
=
X
wul+1 vl (tl ) · xl+1
ul+1 ∈N− (vl ,tl ) (1)
=
X
ωal al+1 · xal+1
al+1 ∈C(al )
with weights ωal al+1 := wul+1 vl (tl ).
(45)
Consequently, T (v, t) is a tree in terms of Theorem C.16. Step 2: Equation (29) holds. Again by (1) the definition of the inflow tree T (v, t) (Definition C.17), (2) Theorem C.16 (which we can apply according to Step 1 ) and (3) Equation (45), we obtain (1)
(1)
cv (t) = cv0 (t0 ) = x0 X (2) =
ω(a0 , ..., alp ) · xlp
p=(a0 ,...,alp )∈P(a0 ,Aleaf )
(1,3)
lp −1
Y
wul+1 vl (tl ) · cvlp (tlp ),
X
=
p=((v0 ,t0 )=(v,t),(u1 ,v1 ,t1 ),...,(ulp ,vlp ,tlp ))∈P(v,t)
l=0
which was to show. E.15
Proof of Theorem C.20
Theorem (Theorem C.20). Let v ∈ V and t ∈ T = Thi ∪ Tfu hold. Let T̂ (v, t) be the interpolation inflow (k) tree of v at time t as defined in Definition C.17. Let k ≥ maxp∈P̂(v,t) lp hold. If ĉv (t) is defined by Equation (11), then ! lp −1 X Y (k−l ) (k) wul+1 vl (tl ) · a(tl+1 )
ĉv (t) =
p=((v0 ,t0 )=(v,t),(u1 ,v1 ,t1 ),...,(ulp ,vlp ,tlp ))∈P̂(v,t)
72
l=0
· ĉvlp p (tlp )
with weights wul+1 vl (tl ) and a(tl+1 ) according to Equation (6) and (16), respectively, holds. Proof. Similar to the proof of Theorem C.19, for each layer l ∈ N0 of the tree, we now associate each node (k−l) al = (ul , vl , tl ) ∈ A of that layer with the magnitude xl := xal := ĉvl (tl ). Step 1: Together with these magnitudes, T̂ (v, t) defines a tree in terms of Theorem C.16. In order to prove the claim, we first of all investigate on the term φ(tl , eul+1 vl , ĉ(k−(l+1)) , ν) for a given parent node al = (ul , vl , tl ) in the l-th layer and its child node al+1 = (ul+1 , vl+1 , tl+1 ) ∈ C(al ). Indeed, by (1) Equation (17), for any l ∈ N0 , φ(tl , eul+1 vl , ĉ(k−(l+1)) , ν) is equal to φ(tl , eul+1 vl , ĉ(k−(l+1)) , ν) ( (k−(l+1)) (k−(l+1)) a(tj ) · ĉul+1 (tj ) + a(tj+1 ) · ĉul+1 (tj+1 ) (1) = (k−(l+1)) (k−(l+1)) a(tj ) · ĉvl (tj ) + a(tj+1 ) · ĉvl (tj+1 )
if ý(ù(tl , eul+1 vl , ν)) ∈ {−luv , luv }, if ý(ù(tl , eul+1 vl , ν)) = 0
for the j ∈ I for which tl − ù(tl , eul+1 vl , ν) ∈ [tj , tj+1 ) holds. Following (1) the definition of the interpolation inflow tree T̂ (v, t) (Definition C.17), we obtain φ(tl , eul+1 vl , ĉ(k−(l+1)) , ν)
P (k−(l+1)) (tl+1 ) a(tl+1 ) · ĉul+1 tl+1 ∈{tj ,tj+1 }: (1) tl −ù(tl ,eul+1 vl ,ν)∈[tj ,tj+1 ) P = (k−(l+1)) (tl+1 ) a(tl+1 ) · ĉvl tl+1 ∈{tj ,tj+1 }:
if ý(ù(tl , eul+1 vl , ν)) ∈ {−luv , luv }, if ý(ù(tl , eeul+1 vl , ν)) = 0
tl −ù(tl ,eul+1 vl ,ν)∈[tj ,tj+1 )
(1)
=
X
(46)
(tl+1 ) a(tl+1 ) · ĉ(k−(l+1)) vl+1
tl+1 ∈{tj ,tj+1 }: tl −ù(tl ,eul+1 vl ,ν)∈[tj ,tj+1 ) (1)
=
X
a(tl+1 ) · xl+1 .
tl+1 ∈{tj ,tj+1 }: tl −ù(tl ,eul+1 vl ,ν)∈[tj ,tj+1 )
Consequently, again by (1) the definition of the interpolation inflow tree T (v, t) (Definition C.17), (2) Equation (11), (3) Equation (46) and (4) basic transformations, for any al ∈ {a0 } ⊔ Aint , we obtain (1)
xal := xl = ĉ(k−l) (tl ) vl X (2) = wul+1 vl (tl ) · φ(tl , evl+1 vl , ĉ(k−(l+1)) , ν) ul+1 ∈N− (vl ,tl ) (3)
=
! X
wul+1 vl (tl ) ·
ul+1 ∈N− (vl ,tl ) (4)
=
(1)
=
X
X
X
ul+1 ∈N− (vl ,tl )
tl+1 ∈{tj ,tj+1 }: tl −ù(tl ,e(l+1)l ,ν)∈[tj ,tj+1 )
X
a(tl+1 ) · xl+1
tl+1 ∈{tj ,tj+1 }: tl −ù(tl ,eul+1 vl ,ν)∈[tj ,tj+1 )
wul+1 vl (tl ) · a(tl+1 ) · xl+1
ωal al+1 · xal+1
al+1 ∈C(al )
with weights ωal al+1 := wul+1 vl (tl ) · a(tl+1 ). 73
(47)
Consequently, T̂ (v, t) is a tree in terms of Theorem C.16. Step 2: Equation (30) holds. Again by (1) the definition of the interpolation inflow tree T̂ (v, t) (Definition C.17), (2) Theorem C.16 (which we can apply according to Step 1 ) and (3) Equation (47), we obtain
(1)
(1)
(k−0) ĉ(k) (t0 ) = xa0 v (t) = ĉv0 X (2) =
ω(a0 , ..., alp ) · xlp
p=(a0 ,...,alp )∈P(a0 ,Aleaf )
lp −1
(1,3)
X
=
Y
p=((v0 ,t0 )=(v,t),(u1 ,v1 ,t1 ),...,(ulp ,vlp ,tlp ))∈P̂(v,t)
p) (tlp ), wul+1 vl (tl ) · a(tl+1 ) · ĉv(k−l lp
l=0
which was to show. E.16
Proof of Theorem 3.6 (k)
Theorem (Theorem 3.6). There exists a κ ∈ N, such that ĉv (t) as defined in Equation (11) is equal to zero for all k > κ, v ∈ V and t ∈ T = Thi ∪ Tfu .
Proof. Let v ∈ V = Vb ⊔ (V \ Vb ) and t ∈ Thi ⊔ Tfu hold. We first identify a node- and time-dependent κ(v, t). (k)
Case 1: If v ∈ Vb or t ∈ Thi holds, by Equation (11), ĉv (t) is equal to zero if k > 0 and equal to (k) cv (t) if k = 0. Therefore, in this case, κ(v, t) = 0. (k)
Case 2.1: If v ∈ V \ Vb , t ∈ Tfu and k ≥ maxp∈P̂(v,t) lp holds, by Theorem C.20, we can express ĉv (t) as (k−l )
a function of values ĉvlp p (tlp ) for paths p = ((v0 , t0 ) = (v, t), (u1 , v1 , t1 ), ..., (ulp , vlp , tlp )) ∈ P̂(v, t). By definition of P̂(v, t) (Definition C.17), vlp ∈ Vb or tlp ∈ Thi holds (and even more, vl ∈ V \ Vb and tl ∈ Tfu holds for all l = 0, ..., lp − 1). Therefore, again by Equation (11), for each path p ∈ P̂(v, t) with length (k−l ) lp ≤ k, ĉvlp p (tlp ) is equal to zero if k > lp and equal to cvlp (tlp ) if k = lp . (k)
Case 2.2: If v ∈ V \ Vb , t ∈ Tfu and k < maxp∈P̂(v,t) lp holds, we can express ĉv (t) as a function of values (k−l )
ĉvlp p (tlp ) for paths p = ((v0 , t0 ) = (v, t), (u1 , v1 , t1 ), ..., (ulp , vlp , tlp )) ∈ P̂(v, t) which satisfy lp ≤ k, and (k−l)
(0)
values ĉvl (tl ) = ĉvl (tl ) for l = k for paths p = ((v0 , t0 ) = (v, t), (u1 , v1 , t1 ), ..., (ul , vl , tl ), ..., (ulp , vlp , tlp )) ∈ P̂(v, t) which satisfy k < lp . For the former paths, the same findings as in Case 2.1 hold. For the latter paths, by definition of P̂(v, t) (Definition C.17), vl ∈ V \ Vb and tl ∈ Tfu holds. Therefore, again by Equation (k−l) (0) (11), for each path p ∈ P̂(v, t) with length lp > k, ĉvl (tl ) = ĉvl (tl ) is equal to zero. (k)
Case 2: In summary, if v ∈ V \ Vb and t ∈ Tfu holds, we can express ĉv (t) as a function of values (k−l ) ĉvlp p (tlp ) which are only non-zero for k ∈ {lp | p ∈ P̂(v, t)}. Therefore, in this case, κ(v, t) = maxp∈P̂(v,t) lp . Finally, we choose κ = maxv∈V maxt∈Thi ⊔Tfu κ(v, t) = maxv∈V maxt∈Thi ⊔Tfu maxp∈P̂(v,t) lp .
Then the
claim follows as long as lp ∈ P̂(v, t) is finite for all v ∈ V and t ∈ Thi ⊔ Tfu . We have discussed that this is the case in Remark C.18. 74
E.17
Proof of Theorem 3.7 and Theorem C.24
Theorem (Theorem 3.7 and Theorem C.24). Let v ∈ V and t ∈ T = Thi ∪ Tfu hold. Let T (v, t) and T̂ (v, t) be the inflow tree and the interpolation inflow tree of v at time t as defined in Definition C.17, respectively. Let Assumption C.21 and C.22 hold. Then the prediction error in-between cv (t) and ĉv (t) is given by |cv (t) − ĉv (t)|
≤
X p̃∈P̂(v,t)
a(t̂1 , ..., t̂lp ) · |cvlp (tlp ) − cvlp (t̂lp )| +|cvlp (tlp )| · |w(t0 , ..., tlp −1 ) − w(t̂0 , ..., t̂lp −1 )| | {z } {z } | Accumulated weighting error
Temporal interpolation error
with joint paths p̃ = ((v0 , t0 ) = (v, t), (u1 , v1 , t1 , t̂1 ), ..., (ulp , vlp , tlp , t̂lp )) ∈ P̂(v, t) according to Definition C.23, lp −1
a(t̂1 , ..., t̂lp ) :=
Y
lp −1
a(t̂l+1 )
w(t̂0 , ..., t̂lp −1 ) :=
and
l=0
Y
wul+1 vl (t̂l )
l=0
with weights wul+1 vl (tl ) and a(tl+1 ) according to Equation (6) and (16), respectively. If the temporal interpolation error is zero, we obtain |cv (t) − ĉv (t)| = 0. Proof. By Theorem C.19 and C.20, lp −1
X
Y
p=((v0 ,t0 )=(v,t),(u1 ,v1 ,t1 ),...,(ulp ,vlp ,tlp ))∈P(v,t)
l=0
cv (t) =
! wul+1 vl (tl )
· cvlp (tlp ) and
lp −1
X
Y
p=((v0 ,t0 )=(v,t),(u1 ,v1 ,t1 ),...,(ulp ,vlp ,tlp ))∈P̂(v,t)
l=0
ĉ(k) v (t) =
! (k−l )
· ĉvlp p (tlp )
wul+1 vl (tl ) · a(tl+1 )
holds for all k ≥ maxp∈P̂(v,t) lp . As already discussed in the proof of Theorem 3.6, for a k < maxp∈P̂(v,t) lp , (k−l )
(k)
paths p ∈ P̂(v, t) with lp ≤ k contribute with the same contribution ĉvlp p (tlp ) = cvlp (tlp ) to ĉv (t) as (k−l)
they would for a larger k, while the ones with lp > k contribute with ĉvl (tl ) = 0 (for l = k < lp ). More (k) specifically, considering a flexible k ∈ N, each path contributes only once to ĉv (t), namely for k = lp . (k−lp ) Therefore, ĉvlp (tlp ) = δklp · cvlp (tlp ) holds for all k ≥ lp and also gives a meaning for cases k < lp . Consequently, by (1) Equation (11) and (2) the findings above, we can express ĉv (t) by (1)
ĉv (t) =
∞ X (k)
ĉv (t)
k=0 (2)
=
lp −1
∞ X
X
Y
k=0 p=((v0 ,t0 )=(v,t),(u1 ,v1 ,t1 ),...,(ul ,vl ,tl ))∈P̂(v,t) p p p
lp −1 (2)
=
X
Y
p=((v0 ,t0 )=(v,t),(u1 ,v1 ,t1 ),...,(ulp ,vlp ,tlp ))∈P̂(v,t)
l=0
! wul+1 vl (tl ) · a(tl+1 )
(k−l )
· ĉvlp p (tlp )
l=0
! wul+1 vl (tl ) · a(tl+1 )
· cvlp (tlp ).
Now we want to use the shared paths p̃ from Equation (31) to express cv (t) and ĉv (t) using the same summands. To achieve this, it is important to note that as an artifact of the interpolation, there are multiple paths p̂ ∈ P̂(v, t) that are identifiable with a single path p ∈ P(v, t). However, due to the fact that certain 75
interpolation weights a(tl+1 ) add up to one, a single contribution cvlp (tlp ) of a path p ∈ P(v, t) to cv (t), weighted by the weight lp −1 Y wul+1 vl (tl ) , l=0
can be distributed to the multiple contributions cvlp (tlp ) of a path p ∈ P(v, t) to cv (t), weighted by the weights lp −1 Y wul+1 vl (tl ) · a(tl+1 ) l=0
associated with corresponding paths p̂ ∈ P̂(v, t) that each are identifiable with the path p ∈ P(v, t). Using the shared paths p̃ notation from Equation (31), this allows to express both cv (t) and ĉv (t) as ! lp −1 X Y wul+1 vl (tl ) · a(t̂l+1 )
cv (t) =
lp −1
X
Y
p̃=((v0 ,t0 )=(v,t),(u1 ,v1 ,t1 ,t̂1 ),...,(ulp ,vlp ,tlp ,t̂lp ))∈P̂(v,t)
l=0
ĉ(k) v (t) =
· cvlp (tlp ) and
l=0
p̃=((v0 ,t0 )=(v,t),(u1 ,v1 ,t1 ,t̂1 ),...,(ulp ,vlp ,tlp ,t̂lp ))∈P̂(v,t)
! wul+1 vl (t̂l ) · a(t̂l+1 )
· cvlp (t̂lp ),
making the sums perfectly comparable. Therefore, by (1) basic transformations and triangle-inequality and (2) the fact that |w1 c1 − w2 c2 | = |w1 c1 + w2 c1 − w2 c1 − w2 c2 | = |w2 c1 − w2 c2 + c1 w1 − c1 w2 | ≤ |c1 − c2 | + |c1 ||w1 − w2 | holds for any w1 , w2 ∈ [0, 1] and c1 , c2 ∈ R, we obtain
(1)
≤
|cv (t) − ĉv (t)| lp −1 lp −1 lp −1 X Y Y Y wul+1 vl (tl ) · cvlp (tlp ) − a(t̂l+1 ) · wul+1 vl (t̂l ) · cvlp (t̂lp ) | | {z } {z } l=0 l=0 l=0
p̃∈P̂(v,t)
∈[0,1]
(2)
≤
lp −1
Y
X p̃∈P̂(v,t)
∈[0,1]
a(t̂l+1 ) ·
|cvlp (tlp ) − cvlp (t̂lp )|
l=0
+
lp −1
|cvlp (tlp )| ·
Y l=0
lp −1
wul+1 vl (tl ) −
Y
! wul+1 vl (t̂l )
.
l=0
If the temporal interpolation error is zero, that is, if ϕ = φ holds, all interpolation times t̂l corresponding to left interpolation bound will be equal to tl hand have interpolation weights a(t̂l ) = 1. Interpolation times t̂l corresponding to right interpolation bounds have interpolation weights a(t̂l ) = 0 (cf. Eq. (16). Therefore, in this case, the only non-zero summands in the error bound are the ones where tl = t̂l holds for all l = 0, ..., lp , in which case |cv (t) − ĉv (t)| = 0 holds. This proves the claim.
76