{'fraction_non_alphanumeric': 0.04381721504817495, 'fraction_numerical': 0.008215727821532803, 'mean_word_length': 4.219300991131977, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 14, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Infinite Time Solutions of Numerical Schemes for Advection Problems

This paper addresses the question whether there are numerical schemes for constant-coefficient advection problems that can yield convergent solutions for an infinite time horizon. The motivation is that such methods may serve as building blocks for long-time accurate solutions in more complex advection-dominated problems. After establishing a new notion of convergence in an infinite time limit of numerical methods, we first show that linear methods cannot meet this convergence criterion. Then we present a new numerical methodology, based on a nonlinear jet scheme framework. We show that these methods do satisfy the new convergence criterion, thus establishing that numerical methods exist that converge on an infinite time horizon, and demonstrate the long-time accuracy gains incurred by this property.


Introduction
The accurate numerical resolution of transport is a crucial ingredient in many computational mathematics/physics/engineering problems. This work focuses on the linear transport under a velocity field that is either given, or arises from other components in a multiphysics computation, such as: (i) the evolution of structures and/or interfaces defined by a level set function φ(x, t) [19,18]; (ii) the passive convection of density fields or marker distributions under a flow field [25]; or (iii) the free-streaming step in particle transport computations [3,5] solved via fractional steps. Moreover, while some applications admit Lagrangian approximations such as particle methods [17,10], this work focuses on the accurate capture of advection via Eulerian methods, i.e., fixed-grid approaches, possibly with adaptive mesh refinement (AMR) [2].
Modern (and future) compute infrastructure, combined with adaptive methodologies (such as AMR), enable the numerical resolution and evolution of small structures (such as droplets in a multi-phase fluid flow simulation, represented by a level set function) over millions of time steps. And with most existing Eulerian methods, even if the structure (level set function) itself remains unchanged (the droplet is merely moved by a locally constant velocity field, i.e., simple transport), the numerical scheme produces a spurious slow change of shape of the solution over time. While in some applications (but not all), overlapping grids [13] and related techniques can mitigate the problem, they do not fully remove the problem of slow deformation over time. Hence, a question of fundamental interest is: are there numerical methods that do not suffer from a slow deterioration of the numerical solution in time, but instead possess convergent infinite time limits? -with the rationale that such methods will be very accurate for long-time computations of transport.
To investigate this fundamental question, we start out at a conceptual level, by considering the constant-coefficient linear advection equation in one space dimension u t + au x = 0 (1) on the periodic domain x ∈ [0, 1], with initial conditions u(x, 0) = u 0 (x). This problem is to be understood as a model problem for the advection of a field quantity (or level set function), u t + a · ∇u = 0, under a locally constant (a = const.) or slowly varying (a(x, t)) velocity field; and the periodic boundary conditions stand as a proxy for a long evolution on a large or unbounded domain.
Despite the simplicity of eq. (1), developing numerical methods that accurately evolve the solution on a fixed grid is a non-trivial task. Of course, for eq. (1) itself, there are methods that are truly exact, namely most standard linear schemes when the CFL number is chosen exactly 1, so that the solution is moved exactly from one grid point to the next in one time step. However, because such a choice is not an option is most practical problems (e.g., due to varying or unknown velocity fields), we explicitly exclude this trivial case and restrict to CFL numbers truly less than 1.
The weaker objective to design numerical schemes that exhibit little deformation of the numerical solution over long times has been pursued, for instance via Fourier continuation methods [1,4]. Due to low dispersive errors, these methods achieve high accuracy for advection-dominated problems compared to traditional finite difference/volume methods. Yet, the methods are unable to generate convergent infinite time limits, for the simple reason that they are linear (see section 3).
To construct numerical methods that do possess convergent infinite time limits, we employ the jet scheme methodology [21], equipped with a special nonlinear interpolant. Jet schemes are semi-Lagrangian fixed-grid approaches that achieve high order by tracking function values and derivatives along characteristic curves, defined via optimally local Hermite interpolation at the foot points of the characteristic curves. This paper is organized as follows. In section 2, we define the new notion of convergence and formulate the mathematical problem. In section 3, we show that linear methods cannot achieve convergent infinite time limits. Then in section 4 we present the new nonlinear jet scheme and establish the numerical analysis used to prove that the method yields convergent infinite time limits. Moreover, it is demonstrated how the new method produces superior longtime accuracy relative to classical nonlinear schemes. We close with a broader discussion and outlook in section 5.

Notion of convergence and formulation of the problem
To study numerical schemes for eq. (1) in terms of their long-time, or even infinite time, accuracy, a departure from established numerical analysis notions of convergence of numerical approximations is needed, as follows. Traditionally, one chooses a fixed finite final time t f , and considers a sequence of meshes, parameterized by a mesh-size h. On each mesh, with time step ∆t = µh a (where 0 < µ < 1 is the fixed Courant-Friedrichs-Lewy (CFL) number [9]) one computes a numerical approximation U h,t f at the fixed final time t f . The numerical scheme is then called convergent if the error ε h, where u t f is the true solution at time t f , the norm · is a suitable function space norm [16], and I(U h,t f ) extends U h,t f to the full domain via a suitable interpolation (see section 3.2 and section 4.1).
The proposed new notion of convergence considers the error as a function of both h and time t n , where t n = n∆t, n ∈ N. For instance, given h with 1 h = m ∈ N and t n , compute the numerical approximation (with ∆t = µh a ) at time t n , denoted U h,tn . Then define the error ε h,tn = I(U h,tn ) − u tn , using the same norm as above. With this bivariate error, the limit h → 0 with nh fixed (i.e., point-wise convergence in time) recovers the traditional notion of convergence. However, one can also require more restrictive notions on the limit h → 0. The ideal situation, which is the strongest notion of convergence that we establish here, is the commuting limits property Here the limits lim sup n→∞ are to be understood in a moving frame of reference x − at, in line with the true solution of eq. (1) which is u(x, t n ) = u 0 (x − at n ), or u(x, t n ) = u 0 (mod(x − at n , 1)) on the periodic domain [0, 1]. Note that the lim sup n→∞ is used because a common occurrence in numerical methods is that they produce a quasi-periodic temporal evolution of the error. So an error limit may not exist in the strict sense, and the lim sup is the right measure of the maximum error over time in such a quasi-periodic case.
Our objective is to have a numerical scheme that produces accurate numerical solutions at any time step. In the case of eq. (1), we want a numerical scheme that mimics the translation property of the advection equation, i.e., the method correctly shifts a numerical approximation from some instance in time to the next time step. Given that on the periodic domain x ∈ [0, 1] the solution returns to its initial configuration every 1 a times, a numerical solution is said to be a fixed point if it agrees at time t+ n a , for some integer n, with its state at time t. Moreover, given an interpolant I that defines a numerical solution on the whole domain (see above), a numerical approximation is said to be a single step fixed point if the scheme precisely moves it from one time step to the next time step. Hence, having a single step fixed point for all steps is a special case of having a fixed point solution. Finally, we call a scheme an exact method if any smooth solution can be approximated by a convergent sequence of fixed points (see section 4.4 for establishing these properties).
We see that exact schemes are characterized by the commuting limits property in eq. (2). If a numerical scheme is convergent in the traditional sense, then for any fixed time with nh fixed, lim h→0 ε h,tn = 0 by definition. So lim sup n→∞ lim h→0 ε h,tn = 0. In order to check whether a convergent method is exact, it suffices to verify that lim h→0 lim sup n→∞ ε h,tn = 0 .

Linear methods
The first question is: can linear schemes have the commuting limits property eq. (2), or even more: do linear schemes exist that precisely shift the numerical solutions (to constantcoefficient linear advection problems) from one time step to next time step? Moreover, if such numerical solutions exist, then are those exact solutions rich enough (i.e., form a dense set) to allow the approximation of any arbitrary initial conditions from a class of functions? Given a grid with mesh size h, let U h,tn denote the state vector of the numerical approximation at the n-th time step. There is a wide variety of numerical methodologies, for example: finite difference methods [16] choose the components of U h,tn to be approximations to the solution at grid points; finite volume methods [11,15] use cell averages; discontinuous Galerkin methods [7,8] use polynomial moments; jet schemes [21,6] use point-wise function values and derivatives. However, what is common in linear methods (i.e., schemes without limiters [24] or other nonlinear components), is that the effects of many steps are fully captured by understanding a single step of the method. Here we consider two frameworks to understand the long time solutions of linear methods. Given h and ∆t related via a CFL number, we either take a finite domain with periodic boundary conditions (to demonstrate numerical examples), or we consider an infinite domain (to carry out a Von Neumann stability analysis [16] without boundary condition artifacts). We consider the recurrence relation to obtain numerical solutions U h,tn = U n j ∞ j=−∞ on an infinite grid at time t n , where the amplification factor g(ξ) is a function of the wave number ξ. Note that, for notational economy, we here omit the superscript h in the components of U h,tn and write n instead of t n . For a finite grid with periodic boundary conditions, X m = {x 0 , x 1 , . . . , x m−1 }, where x j = jh for j = 0, 1, . . . , m − 1 and h = 1 m , using the analogous notation, the update rule of the method can be written as Now g(ξ) assumes only discrete values, and they are exactly the eigenvalues of the matrix M h . We also assume that the scheme uses the same stencil for each grid point. Then the scheme is defined by a vector c = (c 0 , c 1 , . . . , c m−1 ) ∈ R m , and M h becomes a circulant matrix. The scheme's behavior is then determined by the spectral properties of M h . The matrix M h and its j-th eigenvector are given by We obtain the numerical approximation after n steps, T . All the components corresponding to λ j with |λ j | < 1 will decay to zero as n → ∞, and all that remains are the components corresponding to λ j with |λ j | = 1. Let C = {j : |λ j | = 1}, and for all j ∈ C write λ j = e iθ j , θ j ∈ [0, 2π). Then the infinite time limit becomes lim Before we move on to providing a generic proof on an infinite grid, we demonstrate the long time solutions via three illustrative examples on a finite domain.

3.1.
Illustrative examples. Figure 1 shows the (well-known) behavior of the basic upwind method, whose behavior is typical for many other schemes as well. All eigenvalues are truly inside the unit disc, except for one at λ 1 = 1 whose eigenvector corresponds to a constant function (a consequence of having a consistent scheme). Hence, as n → ∞, all modes decay except for the constant mode, and the numerical solution always decays towards a constant function (the average of the initial condition), independent of h, thus lim h→0 lim sup n→∞ ε h,tn = const. Consequently, the scheme is not exact, and the limits h → 0 and lim sup n→∞ do not commute.  Our second example considers a method where the spatial derivative is approximated by standard second order central finite differences (u ), thus the corresponding method of lines matrix has purely imaginary eigenvalues. To make the scheme stable, we need to time-step via an ODE solver whose stability region contains at least some portion of the imaginary axis. We choose the Shu-Osher method [22] (a SSP scheme [12]), whose intercept on the imaginary axis renders the scheme conditionally stable. The eigenvalues of M h (shown in fig. 2) again have only λ 1 = 1 on the unit circle; all other modes decay to zero as n → ∞. So the n → ∞ limit is the same as with upwind. Note that for this scheme (and even more so for methods of higher order), some modes decay slowly (and slower and slower as h → 0). However, the limits h → 0 and lim sup n→∞ nevertheless do not commute.  The last example is a method where the spatial derivative is approximated by the fourth- 12h , leading to a method of lines matrix with all purely imaginary eigenvalues. Now, we time-step via Crank-Nicolson (which we use for the sake of argument -implicit methods are non-popular for advection problems [16]). Because the imaginary axis is exactly the stability boundary for Crank-Nicolson, the resulting scheme is neutrally stable. In contrast to the other examples above, the third scheme shown in fig. 3 has all eigenvalues on the unit circle, i.e., the scheme is purely dispersive. Consequently, all modes contained in the initial (t = 0) approximation persist for all time. However, they move at different speeds, and generally do not reproduce the initial conditions after integer multiples of 1/a, a fact that holds true broadly, see below.
, and eigenvalues (left) and numerical solution at different times (right). All eigenvalues are on the unit circle, hence no modes decay; however, they move at wrong speeds.
3.2. Qualitative analysis. The numerical results above highlight some principles how linear schemes tend to fail to yield commuting limits. To prove that (generic) linear methods in fact do not have commuting limits property eq. (2), we consider problem eq. (1) on the whole real line and assume the following about a linear scheme: Assumption 1.
(i) Consider a sequence of schemes (parameterized by h) with fixed CFL number µ = a∆t h . Moreover, assume 0 < µ < 1 to exclude trivial schemes. (ii) The scheme is defined with constant coefficients c h j for each grid point. (iii) The scheme uses same update stencil for each grid point at every time step, and (for fixed µ) the stencil is independent of h, i.e., c h j = c j ∀ h. (iv) The scheme is stable and consistent.
We employ Von Neumann stability analysis to characterize the stability assumption of the scheme. The analysis is particularly simple for a Cauchy problem [16], i.e., we consider eq. (1) on the real line −∞ < x < ∞ with no boundaries. We discretize the whole domain and approximate the initial condition to get the grid function j is the evaluation of initial condition at grid points x j = jh for j ∈ Z. As above, here we discard the superscript h in the components of u h,0 , for notational economy. As we consider a sequence of meshes and approximate the initial condition on infinite grids hZ, here we appropriately use the scaled 2-norm, u h, that is consistent with the notation in [16]. Then from item v, it follows that u h,0 ∈ 2 (hZ). Now we can express u 0 j by using the semidiscrete Fourier transform [23], defined by the map and its inverse map Note that our physical variable is discrete and unbounded, and corresponds to the space 2 (hZ), . 2 , while the frequency variable is bounded and continuous, and refers to the space L 2 − π h , π h , . L 2 . By item ii and item iii, the one-step update rule for any linear scheme using the stencil [ As the initial condition u 0 ∈ C(R) ∩ L 2 (R), the approximated initial grid function u h,0 is in 2 (hZ). We decompose the initial grid function into Fourier modes by using the semidiscrete Fourier transform and express the grid values as u 0 We now look at the effect of the numerical method on each of the Fourier modes and combine all the modified modes to obtain the numerical solution at some given time. Using equation eq. (3) and the numerical scheme eq. (4) on the grid function U n j = e iξjh of single wave number ξ, we calculate g(ξ) = s ν=−s c ν e iνhξ . Using the recurrence relation eq. (3) on each mode and combining them all we obtain the numerical solution U h,tn at t n given by its components We seek to compare the numerical solution with the true solution in continuum L 2 norm (as we employ the Fourier modes), so it is natural to consider the trigonometric interpolant based on the grid function U n j ∞ j=−∞ that is defined by as our numerical solution function (defined for all x) at t n . A natural (yet not the only possible) way to initialize the numerical scheme is with the discrete data u 0 j ∞ j=−∞ , directly recovered from the initial profile. The numerical schemes now evolve this discrete solution over time, and the interpolant (5) defines a corresponding continuum function.
Here we are only interested in studying the schemes' evolution error, so we employ an error measure that yields zero error at time t = 0: rather than comparing the numerical approximation (5) to the true initial profile u 0 , we instead compare it to a proxy of the true solution, the trigonometric interpolant of the initial grid function, Moreover, like the true solution is just the shift of the initial profile, so is the proxy of the true solution at time t n obtained via shifting the function I(u h,0 )(x) by n∆t to yield Having a stable scheme requires that |g(ξ)| ≤ 1. All components corresponding to ξ with |g(ξ)| < 1 will decay to zero as n → ∞, and all that remains are the components corresponding to ξ with |g(ξ)| = 1. Let C = {ξ : |g(ξ)| = 1}, and for all ξ ∈ C write g(ξ) = e iθ(ξ) . This means that the exact numerical solutions of a linear scheme are precisely the Fourier modes that have their amplification on the unit circle and move with the correct speed as the true solution.
For a linear scheme to have the commuting limits property eq. (2), we would need C to possess increasingly many elements that move with the correct speed, as h → 0 (and thus the number of grid points goes to infinity). However, all linear schemes satisfying assumption 1 have the property that as h → 0, the Fourier modes approaches an algebraic curve in the complex plane. Consequently, one of two scenarios arises: either C has only finitely many elements, or all Fourier modes lie on unit circle. In the latter case, the scheme is purely dispersive; however, expect for trivial cases, the wave speeds of the different Fourier modes do not match up in a way to recover an exact solution. Regarding this result, we now present the following technical lemma, which we are going to use to prove our goal.

Lemma 1.
Under assumption 1, one has where ξ is a wave number, µ is the CFL number, and C is a positive constant.

Then for almost all
otherwise , then clearly d(ξ) ∈ L 2 (R), asû ∈ L 2 − π h , π h and |g(ξ)| ≤ 1, by item iv. Also by construction, d(ξ) is the Fourier transform of the function I(U h,tn ) − I(u h,tn ) (x), for x ∈ R, following from [23]. Using Parseval's identity we can write First taking lim sup n→∞ and then taking the limit as h → 0 of the error ε h,tn , we get

Quantitative scaling.
Having shown that linear schemes are not suitable to produce accurate solutions for really long times, we now provide scaling laws of the errors incurred by these schemes depending on the grid size h and the final time t f . In order to have a chance to obtain an accurate infinite time limit solution, we require that, for a fixed h, the numerical schemes produce non-constant numerical approximations as n → ∞. In practice, we cannot take arbitrarily small ∆t or h, but we are happy to accept h for which the numerical scheme will produce the desired error to exhibit its correct order of convergence. For a given final time t f , there exists a critical value of h = h * , for which the method will produce that desired error. Now, if we take a larger final time t f , then it will be reasonable to believe that we would require a smaller value of h * to produce the desired error. So, the critical value of h depends on the choice of final time t f , and it decreases as t f increases. To rigorously understand this fact we consider (1) and solve this problem by second-order accurate standard Lax-Wendroff method [16]. The modified equation [16] for this numerical method is given by The true solution of (7) captures the behavior of the numerical solutions produced by the Lax-Wendroff method well. As a∆t h = µ is a fixed number (CFL), we can write this equation as v t + av x = Ch 2 v xxx , where C = − a 6 (1−µ 2 ). The differential operators ∂ ∂x and ∂ 3 ∂x 3 commute, which implies that we can solve the equation by splitting into advection and dispersion problems in any order. To understand the error made by the scheme, we can ignore the advection part as it is also present in the original problem (1). So, it suffices to consider only v t = Ch 2 v xxx to understand the change in the numerical solution over time. We compare the solution of this equation, denoted by w h (x, t), with the solution of the original problem with zero velocity, i.e., u(x, t) = u 0 (x). By the time-scaling t → h 2 t, we can deduce that w h (x, t) = w 1 (x, h 2 t). We measure the error for a given final time t f and the corresponding critical resolution h * , by ε h * ,t f = ||w 1 (·, (h * ) 2 t f ) − u 0 (·)|| L 1 in the continuum L 1 norm. We want this error to be a constant which is the desired error to get second order accuracy for the Lax-Wendroff method. The desired error depends on h * and t f , and this forces (h * ) 2 t f = C 1 , for some constant C 1 . In fact, for a p-th order linear scheme, self-similarity arguments based on the modified equation reveal that the error generally scales like ε h,t f ∝ h p t f , i.e., doubling t f also doubles the error, and h * ∝ t −1/p f . Even though for a fixed t f the error decreases as h → 0, for a fixed h the error grows as t f → ∞, thus preventing the commuting limits property (2). The typical shape of the bivariate error function ε h,t f is shown in fig. 5 (left panel), here for the Lax-Wendroff method. For any fixed t f , the error ε h,t f is L-shaped, with a transition mesh size h * (t f ) below which a clean second order manifests.

Jet scheme with nonlinear interpolation
We now demonstrate that there are in fact numerical methods that are exact. To that end, we use the jet scheme formalism, established in [18,21,6] and extended, e.g., in [20]. Jet schemes are semi-Lagrangian approaches that track and evolve derivative information in addition to function values. In the original papers, the derivative information enables one to achieve high order, while maintaining optimally local update rules, based on Hermite interpolation. This concept can also be extended to higher derivatives and higher space dimensions [21]. The jet scheme methodology is as follows. The corresponding errors of only the jet scheme evolution (without i.c. approximation errors) are shown by the dark blue/red curves for (i)/(ii), respectively. For the jet scheme, the total error ε h,t f decreases with h and it is (essentially) constant w.r.t. t f . Hence the scheme has the potential to be exact. Even though there is no traditional convergence rates due to the structure of fixed points (see section 4.3 for a more detailed discussion of the convergence of jet schemes), the jet scheme's evolution error (dark blue) is observed to decrease (in h) more rapidly than the approximation error of the i.c.
1. Given function value (φ) and derivative (ψ) approximations at all grid points at time t, to obtain the state at time t + ∆t, first, for each grid point, solve the characteristic equatioṅ x = a(x, t) backwards from t + ∆t to t to find the associated foot point x ft . 2. Evaluate (based on the grid data) the Hermite interpolant P (x ft ) and its derivative P (x ft ) at this foot point. 3. Then, using these evaluations, evolve the characteristic equationsφ = 0 andψ = −a x (x, t)ψ (if a x = 0, this reduces toψ = 0) forward from t to t + ∆t to yield the new function value and derivative information at the grid points.
Suitable approximations of the initial condition are discussed later in this section.
While the jet scheme methodology works in principle with any Hermite interpolant, it has thus far been established only with interpolations that are linear operators (including in the context of the gradient-augmented level set method [18]). In contrast, here we propose a specific nonlinear interpolant that is tailored to produce the exactness properties requested in section 2. Thus the method is nonlinear as the update rule depends on the state vector nonlinearly. The nonlinear jet scheme interpolant is defined as follows (in 1D). Given data φ(x L ), φ(x R ) (function values) and ψ(x L ), ψ(x R ) (derivatives) we wish to define an interpolant P (x) on [x L , x R ]. We consider the two linear functions emanating from x L and x R : If ψ(x L ) = ψ(x R ) and the lines' intersection point lies inside (x L , x R ), we use the piecewise linear function formed by them (see section 4): φ(x R ) Figure  6. Genuine interpolant. x Otherwise, we use the simple linear interpolant as a fallback case (see fig. 7). This last situation can be interpreted as having an extra line between the segments L L and L R whose length is taken to the maximum. Clearly, the fact that this new interpolant possesses jumps in its derivatives requires particular care in how evaluations of derivatives are defined. We discuss these aspects below.

Two equivalent interpretations of jet schemes.
An important property of jet schemes is that they can be interpreted in two ways [21]: (a) an update rule in a finitedimensional state space; or (b) an update rule of functions. Here we establish that duality for the nonlinear scheme defined above. For a given m = 1/h, let U h,tn = [φ n 1 , ψ n 1 , φ n 2 , ψ n 2 , ..., φ n m , ψ n m ] T ∈ R 2m denote the state vector that defines the jet scheme approximation at time step n. Here φ n j and ψ n j approximate values and derivatives, respectively, of the true solution at the grid point x j . Moreover, let H denote the space of all continuous piecewise linear functions on [0, 1] with finitely many kinks. With these spaces, method 1 can be written as the update rule U h,t n+1 =Ñ U h,tn , whereÑ = E • A • I .
Here A : H → H is the true solution advection operator (which for problem (1) is simply the simple shift operator Au(x) = u(x − a∆t)). It is followed by the evaluation operator E : H → R 2m that evaluates values and derivatives at the grid points with the convention that left-sided derivative values are taken at kinks (note that the particular choice of convention does change the method, but it does not compromise its exactness property); and preceded by the interpolation operator I : R 2m → H, defined on each grid cell [x j , x j+1 ] by the nonlinear interpolant P (x) introduced above.
At t = 0, an initial state vector U h,0 must be defined that represents an approximation I(U h,0 ) to the true initial condition u 0 . Methodologies to do so are discussed in section 4.2. With this setup, conducting many steps of the jet scheme can be interpreted in two equivalent ways: (a) S n = · · · · · · EAI EAI Ñ EAI Ñ (U h,0 ) , or (b) S n = · · · · · · EA IEA N IEA N (I(U h,0 )) .
Here N = I • E • A is the induced jet scheme one-step operator acting between functions, i.e., the jet scheme is the true solution operator A followed by the operator IE. Note that for the linear jet schemes considered in [21], the operator IE is a projection. This last property is generally not given for the nonlinear jet scheme; however, it is not needed for the characterization of the t → ∞ limits of the jet scheme. To conclude, both the finitedimensional update rule (a), which represents how the scheme is implemented, as well as the update rule of functions, (b), can be employed for the numerical analysis carried out in section 4.4.

4.2.
Approximation of the initial condition. There are many ways how one can approximate a given initial condition u 0 via an initial function that is a jet scheme interpolant. For a given m = 1/h, denote the set of grid points by X m = {x 1 , x 2 , . . . , x m }. If u 0 ∈ C 1 , a natural approach, denoted (φ, ψ), is the direct evaluation of values and slopes on the grid, and the assignment φ j = u 0 (x j ) and ψ j = u 0 (x j ) to define the initial jet scheme approximation. This approach, discussed in section 4.5, turns out to be accurate, but it does not guarantee commuting limits.
For that reason, we now introduce a different approximation strategy of the initial condition, denoted (φ δ , ψ δ ), that only requires u 0 ∈ C 0 and that establishes the jet scheme's exactness. First, we evaluate u 0 on the grid points to obtain preliminary valuesφ j = u 0 (x j ) for j = 1 . . . , m. Without loss of generality, we can assume that no three consecutive points are collinear (if there is such a triple, then we can modify the position of the middle point with minimal error to make them non-collinear). Based on these data {(x 1 ,φ 1 ), (x 2 ,φ 2 ), · · · , (x m ,φ m )}, let I(x) be the standard piecewise linear interpolation (that connects the dots, forming kinks at the grid points). We now choose a parameter δ with δ h, where is the machine precision, and consider the function I(x + δ). This is the interpolant I(x) shifted to the left by δ, i.e., its kinks are at x j − δ for j = 1 . . . , m. Hence, I(x + δ) has well-defined slopes at the grid points, and we can assign Below we show that this particular approximation of the initial condition gives rise to the jet scheme's exactness, thus providing a rigorous way to construct methods that possess the commuting limits property. Afterwards, we discuss the pros and cons of other sampling strategies.

4.3.
Quantifying approximation errors and convergence. Conceptualizing convergence of numerical methods requires defining a distance between the true solution (a function) and the sequence of numerical approximations (as h → 0). This is commonly done via a function space norm (e.g., L p ) as in finite elements, or a scaled discrete norm that is induced from a continuum norm as in finite differences [16]. The connection between the h-scaled grid norm and the continuum norm is: (a) the former is a quadrature rule for the latter; and moreover: (b) the former can be defined by means of the latter via an interpolation.
Because jet schemes are fundamentally based on an interpolation step, this last interpretation comes naturally for them: rather than working with a discrete norm on the data (φ j , ψ j ), we directly measure errors in a continuum norm, based on the nonlinear Hermite interpolant defined above. In contrast to the Von Neumann analysis for linear methods in section 3, we do not employ Fourier modes for the analysis of jet schemes. Hence we use the L 1 norm.
An important consequence of measuring errors in a continuum norm is that in general there is a nonzero error between the true initial condition and the initial (t = 0) interpolant that approximates it. Hence, there are two ways to measure the error of the numerical approximation at some time t = t f : the total error, i.e., the L 1 -distance of the interpolant at t = t f to the true solution; and the evolution error, i.e., the L 1 -distance of the interpolant at t = t f to a proxy of the true solution, namely the interpolant from t = 0 translated via (1) to t f . We have already encountered this important distinction in section 3, and it is equally important for jet schemes.
The right panel of fig. 5 presents a convergence study of the jet scheme error ε h,t f as function of h and t f . The light blue/red curves show the total error (at different final times), while the dark blue/red curves display the pure evolution errors (blue vs. red are the two different ways to approximate the initial conditions, see section 4.2). The results show that the total error is generally dominated by the approximation error of the initial condition (the (φ, ψ) strategy yields a factor 2 smaller errors, but the (φ δ , ψ δ ) strategy does not require u 0 to be differentiable), while the scheme's evolution error is much smaller.
One can show that the best one can do to approximate a general smooth initial condition by a piecewise linear function is second order (in the number of kinks), so the jet scheme cannot be better than second order if the total error is considered. However, as we are interested in the (loss of) accuracy as t → ∞ of the scheme itself, we focus on the evolution error.
Because the evolution error is the L 1 -distance between two jet scheme interpolants that are piecewise linear, this error can be calculated up to machine precision by (i) forming subintervals bounded by the kinks (of both interpolants) and the intersection points between the interpolants, and (ii) integrating the difference on each sub-interval exactly. In contrast, the total error cannot be calculated exactly (for general u 0 ), but it can be approximated very accurately by forming similarly splitting grid cells into sub-intervals and then using Gaussian quadrature on each sub-interval (we find that a 3-point Gauss rule suffices here).
Considering the evolution error in the right panel of fig. 5 one can observe that the special (φ δ , ψ δ ) approximation strategy (8) turns out to yield an exact method by means of creating fixed point solutions (see section 4.4) right from the get go (the error (dark red) is essentially machine precision times the number of time steps taken). In contrast, the (φ, ψ) approximation strategy yields smaller total errors, but the evolution errors (dark blue) are larger than with (φ δ , ψ δ ). However, they do decay rapidly as h → 0. Moreover, the scheme appears to be exact nevertheless (the errors do not (up to machine precision) grow with t f ), albeit not via fixed point solutions right away. Comparing those results to the Lax-Wendroff method (left panel of fig. 5), we can see that the exactness property of the jet scheme results in substantially smaller errors, particularly for large final times (the black reference line is the same line in both panels).

4.4.
Numerical analysis for nonlinear jet scheme. Thus far, the exactness of the jet scheme has been based on numerical evidence. Now we establish the exactness of the methods with the strategy (8) via rigorous arguments. Moreover, we provide a pathway towards a novel numerical analysis for the exactness of numerical methods (in one space dimension). Here we assume that the CFL number µ ∈ Q, and let µ = p q with p and q co-prime. It is numerically observed that the restriction µ ∈ Q can be relaxed, but for the analysis of the jet schemes done here, we stick to this mild assumption. We exploit the connection between the two interpretations of jet schemes, established in section 4.1, and employ the function space formulation of the jet scheme to prove its exactness. On the periodic domain [0, 1] with CFL number 0 < µ < 1 and h = 1/m, we have ∆t = µh a . Hence, after N = 1 µh time steps (as µ ∈ Q, we can choose m such that N ∈ N) the true solution has traveled a distance 1 and thus has returned to its initial configuration, i.e., u(x, t + 1 a ) = u(x, t) considering periodicity. We consider the N -step mapping from the function space H defined in section 4.1 into itself that is an approximation to the identity map. Because we are now concerned with the h → 0 convergence of the family of methods, we now use the superscript h to denote the maps' and state vectors' dependence on the resolution. As motivated above, our goal is to characterize the functions that are fixed points of Ψ h , because any such function I(U h ), where U h ∈ R 2m , that satisfies Ψ h I(U h ) = I(U h ) constitutes an exact solution. Moreover, given the family of maps Ψ h , parameterized by h, we have an exact method if for a given true solution u 0 (x), a sequence {I(U h )} h→0 of fixed points (i.e., Ψ h I(U h ) = I(U h )) exists that converges (in L 1 ) to u 0 as h → 0. Now we provide a condition for a function I(U h ) to be a fixed point (of the map Ψ h ) in terms of its kinks' positions. The function I(U h ) can have at most one kink in the interior of each cell, and there could be kinks on the boundary of the cell. We divide each of the m cells into q sub-cells of size h/q, and we consider these resulting mq sub-cells from left to right. We can now determine the positions of all the kinks of the function I(U h ). For kinks that fall on the boundary between adjacent sub-cells, we use the consistent convention to associate them to the left of those two sub-cells.
As established in section 4.1, the one-step map N amounts to first shifting the function I(U h ) to the right by µh via the true solution operator A, followed by the operator IE that replaces the shifted function by the jet scheme interpolant based on its values and slopes evaluated at the grid points. Due to the construction of the jet scheme interpolant, the application of IE leaves the shifted interpolant I(U h )(x − µh) from the prior step unchanged, except on those grid cells that contain two kinks (it is impossible for the shifted interpolant to have more than two kinks per cell). On those cells, the application of IE replaces the shifted interpolant by a new genuine piecewise linear function with a single kink in the cell, or by the fallback function.
A key mechanism to finding fixed points of the N -step mapping Ψ h = N N is to characterize functions for which such discrete changes in the approximating function do not occur. We define: Definition 1. A piecewise linear jet scheme interpolant I(U h ), defined by a state vector U h , is called defective if there exists at least one adjacent pair of kinks that have less than q − 1 sub-cells in between them. Otherwise, it is called non-defective. Proof. As I(U h ) is non-defective, it has at least q −1 sub-cells in between any pair of adjacent kinks. Hence, it is impossible to ever generate a situation where two kinks fall within the same grid cell, if the function I(U h ) is only shifted by integer multiples of the sub-cell width h/q. And in fact, that is the case, because by construction each update step N shifts I(U h ) by µh = p h q , i.e., by exactly p sub-cells.
Based on this intuitive geometric characterization of an important class of fixed points of the jet scheme, we can easily analyze the two initial condition approximation strategies described in section 4.2. While the (φ, ψ) approach generally generates defective interpolants (e.g., two adjacent grid cells may have kinks very close to their shared grid point), the (φ δ , ψ δ ) method, defined in (8), always generates non-defective interpolants, because adjacent kinks are precisely a distance h apart from each other.
x ft x k δ We are now in a position to prove the main theorem, establishing the exactness of the jet scheme with the specific approximation of the initial condition eq. (8).

4.5.
Other initial condition approximation: good news and some caveats. With theorem 4 we have established the exactness of jet schemes, if the (φ δ , ψ δ ) initial function approximation (8) is employed. The underlying reason is that the initial function I(U h,0 ) is a fixed point of the jet scheme update step. However, other approximations of the initial conditions may yield smaller error constants. For instance, fig. 5 shows that the (φ, ψ) strategy, which is natural if u 0 ∈ C 1 as it is based on direct evaluations of the initial condition and its derivatives, turns out to generate second order convergent approximations with smaller errors than the (φ δ , ψ δ ) strategy. A key drawback is that (φ, ψ) strategy generally yields an initial vector U h,0 ∈ R 2m , or function I(U h,0 ), that is not a fixed point of the jet scheme. However, numerical evidence ( fig. 5 and below) suggests that, for a wide range of sufficiently smooth initial functions, the jet scheme does always assume a fixed point after a finite number of update steps, and thus nevertheless gives rise to the commuting limits property eq. (2).
To demonstrate this property, and its consequences, more distinctly, we present a numerical test of the jet scheme with the (φ, ψ) strategy up to a really large final time t f = 100000. fig. 9 shows the numerical results on m = 100 grid points with CFL number µ = 0.9 at t = 100, t = 1000, and the final time. One can see that the jet scheme solution persists up to infinite time. Moreover, we observe (not shown here) that if the computational grid is refined, one always obtains exact solutions after finitely many steps, and they converge (as h → 0) to the true solution. In other words, this scheme is exact for this problem.
Moreover, we compare the long-time evolution of jet schemes with other nonlinear (we already know from section 3 that linear methods have no chance of converging at t → ∞) classical high resolution methods: fifth-order WENO [14], Lax-Wendroff with van Leer limiter [15], and Lax Wendroff with Superbee limiter [15]. In fig. 9, the three snapshots in time are complemented (bottom right panel) by the temporal evolution of the maximum of the numerical solutions. Clearly, all classical schemes keep on deforming in time, with the WENO and van Leer solutions being essentially constant at t f , and the Superbee solution flattening on a very slow time scale. These results demonstrate the jet scheme's superiority over the classical methods in terms of its long-time numerical behavior for constant-coefficient advection problems.
While these observations are very promising, it is actually not true that any possible initial vector U h,0 will be transformed (by the jet scheme) into a fixed point. Here we establish a procedure to construct examples where the jet scheme does not assume a fixed point after a finite number of time steps. To find those, we seek an eigenvector of the mapping Ψ h with associated eigenvalue of modulus strictly less than 1. We do so by restricting the map Ψ h to be of the form M N , where M is a specific one-step update rule that is linear. Because in this case the jet scheme applies the same update rule in every time step, it suffices to analyze just one step of the scheme.
Specifically, consider m = 6 and µ = 3 4 , so that after N = 8 steps the solution has traveled a distance 1. Now assume that the updates on the 6 cells are (from left to right, where 'L' denotes the case that the foot point is left to the kink and 'F' denotes the fallback case): [L L F L L F]. It is straightforward to check that the resulting matrix M has a real eigenvalue with modulus less than 1. Hence, if we start with the associated eigenvector as initial condition, the numerical solution will exponentially decay to zero, and it will never turn into a fixed point in finite time, as shown in fig. 10.  Figure 10. Eigenvector of the one-step matrix M corresponding to an eigenvalue with modulus less than 1. Left: initial condition. Right: solution at t f = 5, which is exactly same as the initial profile, but multiplied by 10 −6 .
Based on this example, we can also construct an infinite family of corresponding examples for larger grid sizes m (multiples of 6), by simply concatenating the sequence [L L F L L F]. However, such functions are by construction highly oscillatory (on the scale h). In contrast, h → 0 sequences of approximations to a smooth initial function are the opposite of such gridscale oscillatory functions; and in fact, decaying eigenvectors like the ones constructed here never arose (for h sufficiently small) in any test cases we conducted. An intuitive explanation for this fact is as follows. As h becomes small, any smooth function locally is O(h 3 ) close to a parabola on a grid cell of size h. If one had an exact parabola (with nonzero quadratic part) on that cell, then the associated piecewise linear jet scheme interpolant would have its kink precisely in the center of the interval. Consequently, for the true function, kinks will be O(h 3 ) away from the center of each cell, except for those few cells in which the function has inflection points. As a consequence, for h sufficiently small, the approximation of a smooth initial condition will be away from those oscillatory eigenvectors. That being said, even for tiny h, one generally cannot expect to start with a fixed point of the jet scheme because near inflection points, fallback cases may arise. However, as discussed above, we do observe that for h sufficiently small, a fixed point is generally assumed after O(1) steps (independent of h).

Conclusions and outlook
The key takeaways from this study are that: (a) there is merit in conceptualizing infinite time limits of numerical methods for advection problems; and (b) schemes that possess convergent infinite time limits do exist. Conceptually, the focus on an infinite (or asymptotically large) time horizon is a departure from classical numerical analysis paradigms, and this work establishes new notions of numerical analysis and formalizes the infinite time convergence via the commuting limits property (2) and the concept of exact methods that possess convergent sequences of fixed point solutions. Moreover, in establishing these new concepts, traditional numerical analysis perspectives have been extended, such as considering error convergence in both mesh size h and final time t f , and the focus on long-time behavior of existing numerical methods.
The specific results and insights found in this work are as follows. First, it is formally shown that linear methods for advection problems cannot (except for trivial cases) yield exact solutions, and thus they are unable to produce approximations that do not deteriorate in time. Second, the formal proof is complemented by asymptotic arguments, based on the modified equation, that establish the scaling of the truncation error ε h,t f interpreted as a bivariate function of h and t f . Third, the jet scheme methodology with a specific nonlinear interpolant is introduced, including a duality between being an update rule in a finite-dimensional state space vs. in a function space. Fourth, exploiting that duality, we prove that the jet scheme, with a specific initial condition approximation, does indeed satisfy the stringent requirement of producing convergent infinite time solutions. And finally, the benefit of the jet schemes' infinite time convergence, relative to existing widely used nonlinear methods, is illustrated in a comparative long-term computation example.
There is a broad variety of questions that are triggered by the concepts and results presented here. One important (open) question is whether there are traditional fixed-grid advection schemes that are exact. While we have proven that exactness is impossible to achieve with linear methods, and have demonstrated numerically that various popular nonlinear methods all fail to be exact, this does not mean that there are no such methods. Certainly, the results in fig. 9 demonstrate that different finite volume limiters may result in vastly different long-time behavior of the numerical solutions. Specifically, while Lax-Wendroff with van Leer limiter clearly produces an essentially constant solution after some large but finite time (t = 10 5 ), Lax-Wendroff with Superbee limiter preserves a (deformed but) clearly non-constant profile even after such a long time. This property plausibly has to do with the fact that the Superbee limiter follows the uppermost boundary of the TVD region [15], and it might hint at the existence of limiters that generate exact methods.
Another important question is to what extent the concepts established herein generalize to more complicated transport problems. Clearly, it is a long route from problem (1) to real application problems that may involve (a) variable coefficients, (b) higher dimensions, and/or (c) coupling of transport to other equations. Regarding (a): while the jet scheme methodology can be applied to variable coefficient advection, the exactness properties fail to carry over. At the same time, simple numerical tests indicate that for non-rapidly varying coefficients, it does take a long time for the numerical solution to degrade. Regarding (b): while linear Hermite interpolants generalize naturally to higher space dimensions [21], there is no such canonical generalization of the piecewise linear interpolant.
Finally, regarding (c), an important question for future research is to what extent exact numerical methods can give rise to more accurate long-time solutions for problems where the accurate resolution of linear advection is a critical sub-problem, such as level set methods in fluid flow simulations or contact discontinuities in compressible gas dynamics. Likewise, one may ask whether some of the presented ideas and concepts can also be extended to other problems with nonlinear traveling waves (solitons, diffusion-reaction, advection-reaction).