An Efficient Lagrangian Interpolation Scheme for Computing Flow Maps and Line Integrals using Discrete Velocity Data

Abstract

We propose a new class of Lagrangian approaches for constructing the flow maps of given dynamical systems. In the case when only discrete velocity data at mesh points is available, an interpolation step will be required. However, in our proposed approaches all particle trajectories share a common global interpolation at each time step and therefore interpolation operations will not increase the overall computational complexity. The old Lagrangian approaches propose to solve the corresponding ordinary differential equations (ODEs) backward in time to obtain the backward flow map. It is inconvenient and not natural, especially when incorporated with certain computational fluid dynamic solvers, because the velocity field needs to be loaded from the terminal time backward to the initial time. In contrast, our proposed approaches for computing the backward flow map propose to solve the corresponding ODEs forward in time which is more practical. We will also extend the proposed approaches to compute line integrals along any particle trajectory. This paper gives a detailed analysis on the computational complexity and error estimate of the proposed Lagrangian approach. Finally, a wide range of applications of our approaches will be given, including the so-called coherent ergodic partition and the high frequency wave propagations based on geometric optic.

1 Introduction

Motivated by the internal wealth and the exceptional practical importance, a growing number of researchers in various fields have been studying efficient tools to help them understand complicated dynamical systems. Most of them, however, have to first finely solve the following two fundamental problems.

Problem F: Given an initial particle location in the phase space at the initial time, how to determine the terminal particle location at a later time and how to compute the line integral of a given integrable function along the particle trajectory.

Problem B: Given a terminal particle location in the phase space at the final time, how to determine its take-off location at the initial time and how to compute the line integral of a given integrable function along the particle trajectory.

These two related yet different problems are of significant importance in many fields such as the level set equation based on the level set method [9, 26, 27, 30], traveltime tomography in seismology [20], high frequency wave propagation [6, 23] and semiclassical solutions to the Schrödinger equations [21, 22]. They are also essential blocks of various recently developed numerical tools for better flow visualizations including the finite size Lyapunov exponent (FSLE) [2, 3, 7, 15, 17], the finite time Lyapunov exponent [12,13,14, 16, 19, 31], the Lagrangian coherent structures [1, 10, 11, 14, 31], the coherent ergodic partition [34], and some studies on chaotic mixing in atmospheric modelings [24, 25, 32].

In general, there are two classes of approaches to solve these two problems. The first class of approaches are the traditional Lagrangian ray tracing approaches which numerically solve the underlying ordinary differential equations (ODEs) governing the dynamic of both the Lagrangian particle and the corresponding line integral along the trajectory. Numerically, these ODEs need to be solved multiple times to obtain the solutions at all mesh points. Let’s take the two-dimensional case for example. Suppose that the computational domain is initially discretized by MN mesh points $\mathbf {x}_{i,j}$ at $t=t_0$ where M, N are the numbers of mesh points in the x direction and y direction, respectively. Then for each mesh point $\mathbf {x}_{i,j}$ with $i=1,2,\ldots ,M$ and $j=1,2,\ldots ,N$, one needs to solve the ODEs one time from the initial time $t=t_0$ to the terminal time $t=T$ to obtain the solutions to Problem F at $\mathbf {x}_{i,j}$. As a result, one needs to solve the ODEs from $t=t_0$ to $t=T$ totally MN times in order to collect the solutions to Problem F at all mesh points. This is fine when the underlying velocity field is smooth enough and defined everywhere in the computational domain, because there are some existing well-developed high order ODE solvers such as ode45 in MATLAB. Sometimes the velocity field might only be discretely defined on a mesh, e.g. the velocity data is from a CFD solver. Under this circumstance, an extra high order interpolation on the discrete velocity data is required at each time step. There are generally two ways of Lagrangian implementation with different orders of complexity as will be discussed in Sect. 2. However, all the old ways of Lagrangian implementation propose to obtain solutions to Problem B by discretizing the corresponding ODEs backward in time. This is numerically inconvenient, especially when incorporated with some CFD solvers, since the velocity field is loaded from the terminal time $t=T$ backward in time to the initial time. This implies that the whole field at all time steps has to be stored in the hard disk, which might not be practical at all.

Another class of approaches are the Eulerian approaches which propose to embed the solutions to Problem F and also Problem B in a higher dimensional space [18, 19, 34]. Based on the level set method [26, 27, 30], they construct a Liouville equation so that the solution to this equation is exactly the flow map. One can solve the corresponding hyperbolic partial differential equations (PDEs) by any well-established high order accurate and robust numerical methods. These Eulerian approaches have an advantage that they only need velocity data defined on the uniform Cartesian mesh and therefore no interpolation to the velocity data is required. However, they have suggested to solve the PDEs backward in time to obtain solutions to Problem F. In particular, one needs to reverse the time and solve the Liouville equations backward from $t = T$ to $t = 0$ with the given terminal condition at $t = T$. Also this implementation is inconvenient and not natural. Furthermore, due to the CFL condition, the constraint for the temporal step size $\varDelta t$ of the Eulerian approaches could be quite strict.

Considering these disadvantages in the original Lagrangian approaches and the Eulerian approaches, in this paper we are proposing and analyzing a new class of Lagrangian interpolation approaches to solve both Problem F and Problem B. There are several main advantages of our proposed new Lagrangian interpolation approaches. Since these approaches are still Lagrangian based methods, the constraint for the temporal step size is much more relaxed than that of the Eulerian approaches. Secondly, the overall computational complexity of our proposed approaches is only $O(KN^2)$ which is optimal in the sense that each grid point is visited for only O(1) time where N is number of mesh points in each spatial direction and K is the number of time steps. Moreover, the convergence behaviour of our proposed approaches can be explicitly given as in Sect. 4. Most importantly, the proposed approaches propose to solve the corresponding ODEs forward in time to obtain the solutions to Problem F and also Problem B which is more natural to implement in practice.

In Sect. 2, we summarize the original two Lagrangian approaches for solving Problem F and Problem B. The proposed Lagrangian interpolation approaches are given in Sect. 3. Some accuracy analysis will be given in Sect. 4. Section 5 will give various examples to demonstrate the effectiveness and significance of the proposed Lagrangian interpolation approaches for computing the flow maps and corresponding line integrals along particle trajectories.

2 Background

In the Lagrangian formulation, the problems are formulated into a system of ordinary differential equations (ODEs) as

$$\begin{aligned} \dot{\mathbf {x}}(t)= & {} \mathbf {u}(\mathbf {x}(t),t) \end{aligned}$$

(1)

$$\begin{aligned} \dot{f}(t)= & {} g(\mathbf {x}(t)), \end{aligned}$$

(2)

where $\mathbf {u}:{\mathbb {R}}^d \times {\mathbb {R}} \rightarrow {\mathbb {R}}^d$ is a given smooth enough velocity field and $g:{\mathbb {R}}^d \rightarrow {\mathbb {R}}$ is a given integrable function. Numerically, one first discretizes the computational domain $[x_{min},x_{max}]\times [y_{min},y_{max}]$ uniformly by $M\times N$ mesh points. To obtain the solutions to Problem F at a particular mesh point $\mathbf {x}_{i,j}$, one needs to first solve ODEs (1) and (2) together from the initial time $t=0$ up to the final time $t=T$ given the initial conditions $\mathbf {x}(0)=\mathbf {x}_{i,j}$ and $f(0)=0$ respectively. The forward flow map of the particle starting from the mesh point $\mathbf {x}_{i,j}$ at $t=0$ to the final time level $t=T$ is then given by $\varPhi _0^T(\mathbf {x}_{i,j})=\mathbf {x}(T)$ and f(T) is the corresponding line integral along this particle trajectory from $t=0$ to $t=T$. In general, one solves these ODEs from $t=0$ up to $t=T$ totally MN times in order to collect the solutions to Problem F at all mesh points $\mathbf {x}_{ij}$ with $i=1,2,\ldots ,M$ and $j=1,2,\ldots ,N$. There are some well-developed numerical methods to integrate these ODEs such as the second order total variation diminishing Runge–Kutta (TVD-RK2) method, the third order total variation diminishing Runge–Kutta (TVD-RK3) method, ODE45 developed in MATLAB, etc. In this section, we review two ways of implementation of the Lagrangian approach in the case when velocity data is only available at mesh points.

2.1 The First Way of Lagrangian Implementation

When velocity data is only available at mesh points, interpolation on the discrete velocity data is required for solving both Problem F and Problem B. The first and also the most natural way is that we use a global interpolation at each time step for each individual particle trajectory starting from a particular mesh point $\mathbf {x}_{i,j}$ at the initial time. Suppose the time interval [0, T] is discretized by $K+1$ points $\{t_n:\,\,\, n=0,1,\ldots ,K\}$ with $t_0=0$ and $t_K=T$. Consider the forward flow map $\varPhi _0^T(\mathbf {x}_{i,j})$ at each mesh point $\mathbf {x}_{i,j}$ where $i=1,2,\ldots ,M$ and $j=1,2,\ldots ,N$. In the first way of Lagrangian implementation, the flow map $\varPhi _0^T(\mathbf {x}_{i,j})$ at all mesh points are computed one by one. For example, one first solves ODEs (1) from $t=0$ to $t=T$ forward in time with the initial condition $\mathbf {x}(0)=\mathbf {x}_{1,1}$ and then assigns $\varPhi _0^T(\mathbf {x}_{1,1})=\mathbf {x}(T)$. After $\varPhi _0^T(\mathbf {x}_{1,1})$ is computed, ODEs (1) are solved again from $t=0$ to $t=T$ forward in time with a different initial condition $\mathbf {x}(0)=\mathbf {x}_{1,2}$ and then we have $\varPhi _0^T(\mathbf {x}_{1,2})=\mathbf {x}(T)$. Repeating this process MN times, we will completely collect the values of $\varPhi _0^T(\mathbf {x}_{i,j})$ for each $i=1,2,\ldots ,M$ and $j=1,2,\ldots ,N$. There exists an order in the computations of $\varPhi _0^T(\mathbf {x}_{i,j})$ for the MN mesh points $\mathbf {x}_{i,j}$. In other words, the flow map of one particle is not computed until the flow map of the former particle in the queue is obtained. This way is very inefficient due to its large computational complexity. In particular, to compute $\varPhi _0^T(\mathbf {x}_{1,1})$, one first determines $\varPhi _0^{t_1}(\mathbf {x}_{1,1})$ which will generally require the velocity data at a non-mesh point at the time level $t=t_1$ and this could be done by interpolating on the given discrete velocity data at $t=t_1$, i.e. $\mathbf {u}(\mathbf {x}_{i,j},t_1)$. This interpolation step takes $O(N^2)$ operations if a global interpolation scheme is used where we have supposed $M=O(N)$. After obtaining $\varPhi _0^{t_1}(\mathbf {x}_{1,1})$, one will continue to compute $\varPhi _{t_1}^{t_2}(\varPhi _0^{t_1}(\mathbf {x}_{1,1}))$, which will also take $O(N^2)$ operations due to interpolation on $\mathbf {u}(\mathbf {x}_{i,j},t_2)$, and therefore the flow map from $t=0$ to $t=t_2$ is obtained by assigning $\varPhi _0^{t_2}(\mathbf {x}_{1,1})=\varPhi _{t_1}^{t_2}(\varPhi _0^{t_1}(\mathbf {x}_{1,1}))$. Repeating this process in turn from $t=t_1$ to $t=t_K=T$, we will obtain the overall flow map $\varPhi _0^T(\mathbf {x}_{1,1})$ which totally takes $K\cdot O(N^2)$ operations because interpolation on $\mathbf {u}(\mathbf {x}_{i,j},t_n)$ takes $O(N^2)$ operations at each intermediate time level $t=t_n$ where $n=1,\ldots ,K$. Similarly, the computation of any other $\varPhi _0^T(\mathbf {x}_{i,j})$ also takes $K\cdot O(N^2)=O(KN^2)$ operations. Summing up them together, the overall computational complexity would be extremely large up to $MN\cdot O(KN^2)=O(KN^4)$.

One can simply reverse the above process to solve Problem B. In particular, for each mesh point $\mathbf {x}_{i,j}$, one solves ODEs (1) and (2) together from the final time $t=T$ backward to the starting time $t=0$ with the given terminal conditions $\mathbf {x}(T)=\mathbf {x}_{i,j}$ and $f(T)=0$ respectively. The backward flow map is then given by $\varPhi _T^0(\mathbf {x}_{i,j})=\mathbf {x}(0)$ and the corresponding line integral along this particle trajectory from $t=T$ back to $t=0$ is given by f(0). Again these ODEs need to be solved from $t=T$ back to $t=0$ totally MN times in order to obtain the solutions to Problem B at all mesh points $\mathbf {x}_{ij}$ with $i=1,2,\ldots ,M$ and $j=1,2,\ldots ,N$. If velocity data is only given on the mesh, solving Problem B will then have the same computational complexity as in solving Problem F. Furthermore, because this implementation proposes to solve ODEs backward in time, the velocity field needs to be loaded from the terminal time backward to the initial time which is not natural and inconvenient in practice.

2.2 The Second Way of Lagrangian Implementation

There is an alternative way of Lagrangian implementation which can greatly reduce the computational complexity. Unlike in the first way where the forward flow map at different mesh points are computed one by one in certain order, the second way computes simultaneously the flow map $\varPhi _0^T(\mathbf {x}_{i,j})$ for $i=1,2,\ldots ,M$ and $j=1,2,\ldots ,N$. In particular, one first computes the one step forward flow map $\varPhi _0^{t_1}(\mathbf {x}_{i,j})$ at those MN mesh points from $t=0$ to $t=t_1$ simultaneously, and then the one step flow map $\varPhi _{t_1}^{t_2}(\varPhi _0^{t_1}(\mathbf {x}_{i,j}))$ from $t=t_1$ to $t=t_2$ are computed simultaneously. That is, $\varPhi _0^{t_2}(\mathbf {x}_{i,j})$ at all mesh points are obtained at once. In general, as soon as the intermediate forward flow map $\varPhi _0^{t_n}(\mathbf {x}_{i,j})$ are obtained for all $\mathbf {x}_{i,j}$ where $n=1,2,\ldots ,K-1$, the one step forward flow map $\varPhi _{t_n}^{t_{n+1}}(\varPhi _0^{t_n}(\mathbf {x}_{i,j}))$ for all $\mathbf {x}_{i,j}$ are computed simultaneously and then $\varPhi _0^{t_{n+1}}(\mathbf {x}_{i,j})$ are immediately given by $\varPhi _0^{t_{n+1}}(\mathbf {x}_{i,j})=\varPhi _{t_n}^{t_{n+1}}(\varPhi _0^{t_n}(\mathbf {x}_{i,j}))$.

In the implementation, the one step forward flow map $\varPhi _{t_n}^{t_{n+1}}(\varPhi _0^{t_n}(\mathbf {x}_{i,j}))$ at all mesh points $\mathbf {x}_{i,j}$ are computed simultaneously at any particular time level $t=t_n$ where $n=0,1,\ldots ,K-1$. In particular, one first discretizes the ODEs (1) MN times from $t=t_n$ to $t=t_{n+1}$ with MN different initial conditions $\mathbf {x}(t_n)=\varPhi _0^{t_n}(\mathbf {x}_{i,j})$ and then assigns $\varPhi _{t_n}^{t_{n+1}}(\varPhi _0^{t_n}(\mathbf {x}_{i,j}))=\mathbf {x}(t_{n+1})$. This discretization step could be quickly done. For example, if we use the simple first order Euler method for the discretization, we have

$$\begin{aligned} \frac{\mathbf {x}(t_{n+1})-\mathbf {x}(t_n)}{\varDelta t} = \mathbf {u}(\mathbf {x}(t_n),t_n) \end{aligned}$$

where $\mathbf {x}(t_n)=\varPhi _0^{t_n}(\mathbf {x}_{i,j})$. After that $\varPhi _{t_n}^{t_{n+1}}(\varPhi _0^{t_n}(\mathbf {x}_{i,j}))$ are explicitly given as

$$\begin{aligned} \varPhi _{t_n}^{t_{n+1}}\left( \varPhi _0^{t_n}(\mathbf {x}_{i,j})\right) =\mathbf {x}(t_{n+1})=\varPhi _0^{t_n}(\mathbf {x}_{i,j})+ \varDelta t \, \mathbf {u}\left( \varPhi _0^{t_n}(\mathbf {x}_{i,j}),t_n\right) . \end{aligned}$$

(3)

The implementation of the above formula for all $\mathbf {x}_{i,j}$ can be simultaneously and easily done by a simple matrix operation in MATLAB if $\mathbf {u}(\varPhi _0^{t_n}(\mathbf {x}_{i,j}),t_n)$ is known. Fortunately, $\mathbf {u}(\varPhi _0^{t_n}(\mathbf {x}_{i,j}),t_n)$ can be obtained by a single interpolation on $\mathbf {u}(\mathbf {x}_{i,j},t_n)$. For example, one can use the interp2 function in MATLAB for the interpolation. The implementation of (3) and the interpolation on $\mathbf {u}(\mathbf {x}_{i,j},t_n)$ both require O(MN) operations at each time step $t=t_n$. As a result, the total complexity at each time step remains O(MN) and summing up K time steps together the overall computational complexity of the second way of Lagrangian implementation is only $K\cdot O(MN)=O(KN^2)$ which has greatly reduced the computational complexity in the first way. We will give the detailed discussion for computational complexities in Sect. 3.4. Higher order schemes, such as TVD-RK2 [28], can also be used for the discretization while keeping the computational complexity to be $O(KN^2)$.

To solve Problem B, we can also reverse the above process. In particular, we first compute $\varPhi _{t_K}^{t_{K-1}}(\mathbf {x}_{i,j})$ for all $\mathbf {x}_{i,j}$ simultaneously and then compute $\varPhi _{t_{K-1}}^{t_{K-2}}(\varPhi _{t_K}^{t_{K-1}}(\mathbf {x}_{i,j}))$ for all $\mathbf {x}_{i,j}$ simultaneously where interpolation on discrete velocity data might be required. Then we have actually obtained $\varPhi _{t_K}^{t_{K-2}}(\mathbf {x}_{i,j})$ for all $\mathbf {x}_{i,j}$ since $\varPhi _{t_K}^{t_{K-2}}(\mathbf {x}_{i,j})=\varPhi _{t_{K-1}}^{t_{K-2}}(\varPhi _{t_K}^{t_{K-1}}(\mathbf {x}_{i,j}))$. Continue this process till $t=0$ and we will obtain the overall backward flow map $\varPhi _{t_K}^0(\mathbf {x}_{i,j})$ from $t=t_K$ to $t=0$ for all $\mathbf {x}_{i,j}$.

3 Our New Lagrangian Interpolation Approaches

The second way of Lagrangian implementation is efficient in solving Problem F. When applied to Problem B, however, it needs to discretize the corresponding ODEs backward in time from $t=T$ to $t=0$ in order to obtain the corresponding solutions, e.g. the overall backward flow map $\varPhi _T^0(\mathbf {x}_{i,j})$. This kind of implementation is numerically inconvenient, especially when incorporated with some CFD solvers, since the velocity field is loaded from the terminal time $t=T$ backward in time to the initial time. This implies that the whole field at all time steps has to be stored in the hard disk, which might not be practical at all. In this section, we will propose a new Lagrangian interpolation approach for the Problem F. The proposed approach has two advantages. First it can be extended to solve Problem B more naturally since it discretizes the corresponding ODEs forward in time to obtain the solutions to Problem B. Furthermore, the convergence of the proposed approach can be theoretically guaranteed as will be given in Sect. 4.

3.1 A New Lagrangian Interpolation Approach for the Problem F

We first propose a new Lagrangian interpolation approach to efficiently compute the forward flow map $\varPhi _0^T(\mathbf {x}_{i,j})$ at all the MN mesh points $\mathbf {x}_{i,j}$ based on discrete velocity data given only on the mesh. Like the second way of Lagrangian implementation introduced in Sect. 2.2, we also propose to compute simultaneously the flow map $\varPhi _0^T(\mathbf {x}_{i,j})$ for all $i=1,2,\ldots ,M$ and $j=1,2,\ldots ,N$. However, at each time level $t=t_n$, the one step forward flow map we compute is now $\varPhi _{t_n}^{t_{n+1}}(\mathbf {x}_{i,j})$, rather than $\varPhi _{t_n}^{t_{n+1}}(\varPhi _0^{t_n}(\mathbf {x}_{i,j}))$. In particular, we first compute the one step forward flow map $\varPhi _0^{t_1}(\mathbf {x}_{i,j})$ at those totally MN mesh points from $t=0$ to $t=t_1$ simultaneously, and then the one step flow map $\varPhi _{t_1}^{t_2}(\mathbf {x}_{i,j})$ from $t=t_1$ to $t=t_2$ are computed simultaneously. After that the flow map $\varPhi _0^{t_2}(\mathbf {x}_{i,j})$ from $t=0$ to $t=t_2$ at MN mesh points can be constructed as a composition of those already computed one step flow maps $\varPhi _0^{t_2}(\mathbf {x}_{i,j})=\varPhi _{t_1}^{t_2}\circ \varPhi _0^{t_1}(\mathbf {x}_{i,j})$. We then subsequently compute the one step flow map $\varPhi _{t_2}^{t_3}(\mathbf {x}_{i,j})$ from $t=t_2$ to $t=t_3$ for all the MN mesh points $\mathbf {x}_{i,j}$ simultaneously and thus the forward flow maps from $t=0$ to $t=t_3$ can be obtained by $\varPhi _0^{t_3}(\mathbf {x}_{i,j})=\varPhi _{t_2}^{t_3}\circ \varPhi _0^{t_2}(\mathbf {x}_{i,j})$. In general, as soon as we obtain the intermediate forward flow map $\varPhi _0^{t_n}(\mathbf {x}_{i,j})$ for all the MN mesh points $\mathbf {x}_{i,j}$ where $n=1,2,\ldots ,K-1$, we will compute the one step forward flow map $\varPhi _{t_n}^{t_{n+1}}(\mathbf {x}_{i,j})$ and then use the composition to obtain $\varPhi _0^{t_{n+1}}(\mathbf {x}_{i,j})=\varPhi _{t_n}^{t_{n+1}}\circ \varPhi _0^{t_n}(\mathbf {x}_{i,j})$.

In the proposed new Lagrangian interpolation approach, the two key points are how to efficiently compute the one step forward flow map $\varPhi _{t_n}^{t_{n+1}}(\mathbf {x}_{i,j})$ at the MN mesh points simultaneously and how to implement the composition $\varPhi _0^{t_{n+1}}(\mathbf {x}_{i,j})=\varPhi _{t_n}^{t_{n+1}}\circ \varPhi _0^{t_n}(\mathbf {x}_{i,j})$ given $\varPhi _0^{t_n}(\mathbf {x}_{i,j})$ and $\varPhi _{t_n}^{t_{n+1}}(\mathbf {x}_{i,j})$ at each particular time level $t=t_n$ with $n\in \{1,2,\ldots ,K-1\}$. Actually at each time level $t=t_n$, the one step forward flow map $\varPhi _{t_n}^{t_{n+1}}(\mathbf {x}_{i,j})$ at each $\mathbf {x}_{i,j}$ can be easily computed simultaneously by discretizing the ODEs (1) from $t=t_n$ to $t=t_{n+1}$ and the composition $\varPhi _{t_n}^{t_{n+1}}\circ \varPhi _0^{t_n}(\mathbf {x}_{i,j})$ can be done by an interpolation following the idea in [5]. In particular, we first determine the one step forward flow map $\varPhi _0^{t_1}(\mathbf {x}_{i,j})$ at each $\mathbf {x}_{i,j}$ by discretizing the ODEs (1) from $t=0$ to $t=t_1$. In particular, if the first order Euler method is used for the discretization, then for each particular $\mathbf {x}_{i,j}$, we have

$$\begin{aligned} \frac{\varPhi _0^{t_1}(\mathbf {x}_{i,j})-\mathbf {x}_{i,j}}{\varDelta t} = \mathbf {u}(\mathbf {x}_{i,j},0) \end{aligned}$$

and therefore the flow map $\varPhi _0^{t_1}(\mathbf {x}_{i,j})$ are explicitly given by

$$\begin{aligned} \varPhi _0^{t_1}(\mathbf {x}_{i,j})= \mathbf {x}_{i,j}+ \varDelta t \, \mathbf {u}(\mathbf {x}_{i,j},0). \end{aligned}$$

It’s important to see that a simple implementation of the above formula, e.g. a simple matrix operation in MATLAB, will give us the one step forward flow map $\varPhi _0^{t_1}(\mathbf {x}_{i,j})$ for any $\mathbf {x}_{i,j}$ simultaneously where $i=1,2,\ldots ,M$ and $j=1,2,\ldots ,N$. In general, as soon as we have $\varPhi _0^{t_n}(\mathbf {x}_{i,j})$ where $n=1,2,\ldots ,K-1$, we will discretize the ODEs (1) from $t=t_n$ to $t=t_{n+1}$ to obtain the one step flow map $\varPhi _{t_n}^{t_{n+1}}(\mathbf {x}_{i,j})$. Also if we use the first order Euler method for the discretization, we have

$$\begin{aligned} \frac{\varPhi _{t_n}^{t_{n+1}}(\mathbf {x}_{i,j})-\mathbf {x}_{i,j}}{\varDelta t} = \mathbf {u}(\mathbf {x}_{i,j},t_n) \end{aligned}$$

and $\varPhi _{t_n}^{t_{n+1}}(\mathbf {x}_{i,j})$ are then explicitly given as

$$\begin{aligned} \varPhi _{t_n}^{t_{n+1}}(\mathbf {x}_{i,j})= \mathbf {x}_{i,j}+ \varDelta t \, \mathbf {u}(\mathbf {x}_{i,j},t_n). \end{aligned}$$

(4)

Then we propose to construct the forward flow map $\varPhi _0^{t_{n+1}}(\mathbf {x}_{i,j})$ from $t=0$ to $t=t_{n+1}$ for all the MN mesh points $\mathbf {x}_{i,j}$ by a composition of the computed flow maps $\varPhi _0^{t_n}(\mathbf {x}_{i,j})$ and $\varPhi _{t_n}^{t_{n+1}}(\mathbf {x}_{i,j})$, i.e. $\varPhi _0^{t_{n+1}}(\mathbf {x}_{i,j})=\varPhi _{t_n}^{t_{n+1}}\circ \varPhi _0^{t_n}(\mathbf {x}_{i,j})$. This composition can be easily implemented by a single interpolation. For example, one can use the interp2 function in MATLAB for the interpolation because the computed flow maps $\varPhi _0^{t_n}(\mathbf {x}_{i,j})$ and $\varPhi _{t_n}^{t_{n+1}}(\mathbf {x}_{i,j})$ are uniformly defined on mesh points and thus can be naturally represented in the matrix form. Continue these steps in order for $n=1,2,\ldots ,K-1$ and we will finally obtain the overall forward flow map $\varPhi _0^T(\mathbf {x}_{i,j})=\varPhi _{t_{K-1}}^{t_K} \circ \cdots \circ \varPhi _{t_1}^{t_2}\circ \varPhi _0^{t_1}(\mathbf {x}_{i,j})$ for all mesh points $\mathbf {x}_{i,j}$ simultaneously. In this proposed Lagrangian approach, the computation of one step flow map $\varPhi _{t_n}^{t_{n+1}}(\mathbf {x}_{i,j})$ requires only O(MN) operations at each time step $t=t_n$ because the computation is actually a single matrix operation (see (4) for example). The composition $\varPhi _0^{t_{n+1}}(\mathbf {x}_{i,j})=\varPhi _{t_n}^{t_{n+1}}\circ \varPhi _0^{t_n}(\mathbf {x}_{i,j})$ at each time step $t=t_n$ also requires only O(MN) operations because it is essentially a matrix interpolation operation. As a result, the total complexity at each time step remains O(MN) and summing up K time steps together the overall computational complexity of our proposed Lagrangian interpolation approach is $O(KMN)=O(KN^2)$ which is the same as that of the second way of Lagrangian implementation. We will give the detailed discussion for computational complexities in Sect. 3.4.

We can also use higher order schemes to compute those one step forward flow map $\varPhi _{t_n}^{t_{n+1}}(\mathbf {x}_{i,j})$ at each time step $t=t_n$. For example, if the TVD-RK2 method is used [28], the forward flow map $\varPhi _{t_n}^{t_{n+1}}(\mathbf {x}_{i,j})$ at each mesh point $\mathbf {x}_{i,j}$ is computed by

$$\begin{aligned} \frac{\hat{\mathbf {x}}(t_{n+1})-\mathbf {x}_{i,j}}{\varDelta t}= & {} \mathbf {u}(\mathbf {x}_{i,j},t_n) \\ \frac{\hat{\mathbf {x}}(t_{n+2})-\hat{\mathbf {x}}(t_{n+1})}{\varDelta t}= & {} \mathbf {u}(\hat{\mathbf {x}}(t_{n+1}),t_{n+1}) \\ \varPhi _{t_n}^{t_{n+1}}(\mathbf {x}_{i,j})= & {} \frac{1}{2} \left( \mathbf {x}_{i,j} +\hat{\mathbf {x}}(t_{n+2}) \right) . \end{aligned}$$

Here functions $\hat{\mathbf {x}}(t_{n+1})$ and $\hat{\mathbf {x}}(t_{n+2})$ are intermediate predicted solutions at $t_{n+1}$ and $t_{n+2}$, respectively. More explicitly, we have

$$\begin{aligned} \varPhi _{t_n}^{t_{n+1}}(\mathbf {x}_{i,j})=\mathbf {x}_{i,j}+\frac{\varDelta t}{2}[\mathbf {u}(\mathbf {x}_{i,j},t_n)+\mathbf {u}(\mathbf {x}_{i,j}+\varDelta t\mathbf {u}(\mathbf {x}_{i,j},t_n),t_{n+1})]\,. \end{aligned}$$

(5)

Also we propose to implement (5) in the matrix form to obtain $\varPhi _{t_n}^{t_{n+1}}(\mathbf {x}_{i,j})$ for all $\mathbf {x}_{i,j}$ simultaneously. Compared to formula (4) which has a complexity of O(MN), formula (5) seems a little more complex because the term $\mathbf {u}(\mathbf {x}_{i,j}+\varDelta t\mathbf {u}(\mathbf {x}_{i,j},t_n),t_{n+1})$ requires an interpolation on the discrete velocity data $\mathbf {u}(\mathbf {x}_{i,j},t_{n+1})$ given at $t=t_{n+1}$. However, this interpolation step requires only O(MN) operations and thus formula (5) has kept the computational complexity to be O(MN).

Remark

In the implementation of Algorithm 1, we have supposed that those one step forward flow maps are computed using the TVD-RK2 scheme. Other discretization schemes can also be used and the discussion is similar.

Concerning the line integral along the particle trajectory, we consider the following ODE given by

$$\begin{aligned} \frac{d \hat{f}(\mathbf {x}_{i,j},t)}{dt} = g(\mathbf {x}(t)) \end{aligned}$$

(6)

with an integrable function $g(\mathbf {x})$. Numerically, we couple ODE (6) with the ODEs (1) and solve them from $t=t_n$ to $t=t_{n+1}$ forward in time with the initial conditions given by $\hat{f}(\mathbf {x}_{i,j},t_n)=0$ and $\mathbf {x}(t_n)=\mathbf {x}_{i,j}$, respectively. The solution $\hat{f}(\mathbf {x}_{i,j},t_{n+1})$ represents the line integral of the given function g along the particle trajectory starting from $\mathbf {x}_{i,j}$ at time $t=t_{n}$ reaching at the time level $t=t_{n+1}$. We can also use any well-developed numerical integrator to discretize ODE (6). In particular, if the ODE is discretized with the TVD-RK2 scheme, we have

$$\begin{aligned} \hat{f}(\mathbf {x}_{i,j},t_{n+1}) =\frac{\varDelta t}{2} [g(\mathbf {x}_{i,j})+g(\mathbf {x}_{i,j}+\varDelta t\mathbf {u}(\mathbf {x}_{i,j},t_n))]\,. \end{aligned}$$

(7)

Once we have the forward flow map, we interpolate and use the recursion to obtain

$$\begin{aligned} f(\mathbf {x}_{i,j},t_{n+1}) = \hat{f}\left( \varPhi _0^{t_n}(\mathbf {x}_{i,j}),t_{n+1}\right) + f(\mathbf {x}_{i,j},t_n) \end{aligned}$$

(8)

which represents the line integral of the given integrable function g along the particle trajectory reaching the time $t=t_{n+1}$ starting from the location $\mathbf {x}_{i,j}$ at time $t=0$.

3.2 A New Lagrangian Interpolation Approach for the Problem B

No matter in the first way or the second way of Lagrangian implementation as introduced in Sect. 2, solutions to Problem B are obtained by discretizing the corresponding ODEs backward in time. For example, one needs to solve the ODEs (1) backward in time from $t=T$ to $t=0$ in order to obtain the overall backward flow map $\varPhi _T^0(\mathbf {x})$ from $t=T$ to $t=0$. As we said before, this implementation is numerically inconvenient, especially when incorporated with some CFD solvers, since the velocity field is loaded from the terminal time $t=T$ backward in time to the initial time. This implies that the whole field at all time steps has to be stored in the hard disk, which might not be practical at all. In this part, following the ideas of our proposed new approach for Problem F, we propose a new Lagrangian approach for Problem B. In this new appraoch, we are solving the corresponding ODEs forward in time to obtain the solutions to Problem B, which is more natural and convenient to numerically implement.

We first propose a new Lagrangian interpolation approach to compute the backward flow map on the fly. That is, we will discretize the ODEs (1) forward in time to obtain the backward flow map. In particular, we first determine the one step backward flow map $\varPhi ^0_{t_1}(\mathbf {x}_{i,j})$ for all $\mathbf {x}_{i,j}$ simultaneously and then compute the one step backward flow map $\varPhi ^{t_1}_{t_2}(\mathbf {x}_{i,j})$ for all $\mathbf {x}_{i,j}$. Once we have $\varPhi ^0_{t_1}(\mathbf {x}_{i,j})$ and $\varPhi ^{t_1}_{t_2}(\mathbf {x}_{i,j})$ for all $\mathbf {x}_{i,j}$, the backward flow map for all mesh points from $t=t_2$ to $t=0$ can be easily constructed by a composition $\varPhi ^0_{t_2}(\mathbf {x}_{i,j})=\varPhi _{t_1}^0\circ \varPhi ^{t_1}_{t_2}(\mathbf {x}_{i,j})$. Generally, as soon as we obtain the intermediate backward flow map $\varPhi ^0_{t_n}(\mathbf {x}_{i,j})$ for any mesh point $\mathbf {x}_{i,j}$ where $n=1,2,\ldots ,K-1$, we will compute the one step backward flow map $\varPhi ^{t_n}_{t_{n+1}}(\mathbf {x}_{i,j})$ and then use the composition to obtain $\varPhi ^0_{t_{n+1}}(\mathbf {x}_{i,j})=\varPhi _{t_n}^0\circ \varPhi ^{t_n}_{t_{n+1}}(\mathbf {x}_{i,j})$. Finally we can obtain the overall backward flow map $\varPhi ^0_T(\mathbf {x}_{i,j})$ from $t=T$ to $t=0$ for all mesh points $\mathbf {x}_{i,j}$. In this new approach, the whole computational process is now marching forward in time from $t=0$ to $t=T$ while the computational complexity is kept to be $O(KN^2)$.

Also we have to properly deal with two points when implementing this proposed new approach for computing the backward flow map. The first point is how to efficiently compute the one step backward flow map $\varPhi _{t_{n+1}}^{t_n}(\mathbf {x}_{i,j})$ at all mesh points simultaneously for any $n\in \{0,1,\ldots ,K-1\}$. Actually the one step backward flow map $\varPhi _{t_{n+1}}^{t_n}(\mathbf {x}_{i,j})$ at each $\mathbf {x}_{i,j}$ can be computed simultaneously by easily discretizing the ODEs (1) backward from $t=t_{n+1}$ to $t=t_n$. If the first order Euler method is used for the discretization, then we have

$$\begin{aligned} \frac{\mathbf {x}_{i,j}-\varPhi _{t_{n+1}}^{t_n}(\mathbf {x}_{i,j})}{\varDelta t} = \mathbf {u}(\mathbf {x}_{i,j},t_{n+1}) \end{aligned}$$

and $\varPhi _{t_{n+1}}^{t_n}(\mathbf {x}_{i,j})$ can be explicitly given as $\varPhi _{t_{n+1}}^{t_n}(\mathbf {x}_{i,j})= \mathbf {x}_{i,j}- \varDelta t \mathbf {u}(\mathbf {x}_{i,j},t_{n+1})$. If the TVD-RK2 method is used, we compute the backward flow map $\varPhi _{t_{n+1}}^{t_n}(\mathbf {x}_{i,j})$ at each mesh point $\mathbf {x}_{i,j}$ by

$$\begin{aligned} \frac{\hat{\mathbf {x}}(t_n)-\mathbf {x}_{i,j}}{\varDelta t}= & {} -\mathbf {u}(\mathbf {x}_{i,j},t_{n+1}) \\ \frac{\hat{\mathbf {x}}(t_{n-1})-\hat{\mathbf {x}}(t_{n})}{\varDelta t}= & {} -\mathbf {u}(\hat{\mathbf {x}}(t_{n}),t_{n}) \\ \varPhi _{t_{n+1}}^{t_n}(\mathbf {x}_{i,j})= & {} \frac{1}{2} \left( \mathbf {x}_{i,j} +\hat{\mathbf {x}}(t_{n-1}) \right) \,. \end{aligned}$$

Here Functions $\hat{\mathbf {x}}(t_n)$ and $\hat{\mathbf {x}}(t_{n-1})$ are intermediate predicted solutions at $t_n$ and $t_{n-1}$, respectively. The flow map is then explicitly given as

$$\begin{aligned} \varPhi _{t_{n+1}}^{t_n}(\mathbf {x}_{i,j})=\mathbf {x}_{i,j}-\frac{\varDelta t}{2}[\mathbf {u}(\mathbf {x}_{i,j},t_{n+1})+\mathbf {u}(\mathbf {x}_{i,j}-\varDelta t\mathbf {u}(\mathbf {x}_{i,j},t_{n+1}),t_n)]\,. \end{aligned}$$

(9)

The second point is how to implement the composition $\varPhi ^0_{t_{n+1}}(\mathbf {x}_{i,j})=\varPhi _{t_n}^0\circ \varPhi ^{t_n}_{t_{n+1}}(\mathbf {x}_{i,j})$ for any $n\in \{1,2,\ldots ,K-1\}$. Similarly this composition can be done by a single interpolation. For example, we can use the interp2 function in MATLAB for the interpolation based on the computed flow maps $\varPhi ^{t_n}_{t_{n+1}}(\mathbf {x}_{i,j})$ and $\varPhi _{t_n}^0(\mathbf {x}_{i,j})$ defined on uniform mesh points.

Concerning the line integral along the particle trajectory, we follow a similar approach introduced in Sect. 3.1. Numerically, we couple ODE (6) with the ODEs (1) and solve them from $t=t_{n+1}$ to $t=t_n$ backward in time with the terminal conditions given by $\hat{f}(\mathbf {x}_{i,j},t_{n+1})=0$ and $\mathbf {x}(t_{n+1})=\mathbf {x}_{i,j}$, respectively. The solution $\hat{f}(\mathbf {x}_{i,j},t_n)$ represents the negative of the line integral of the given function g along the particle trajectory reaching $\mathbf {x}_{i,j}$ at the time level $t=t_{n+1}$ starting from the time $t=t_n$. If ODE (6) is discretized with the TVD-RK2 scheme, we have

$$\begin{aligned} \hat{f}(\mathbf {x}_{i,j},t_n) =-\frac{\varDelta t}{2} [g(\mathbf {x}_{i,j})+g(\mathbf {x}_{i,j}-\varDelta t\mathbf {u}(\mathbf {x}_{i,j},t_{n+1}))]\,. \end{aligned}$$

(10)

This implies

$$\begin{aligned} f(\mathbf {x}_{i,j},t_{n+1}) = f\left( \varPhi _{t_{n+1}}^{t_n}(\mathbf {x}_{i,j}),t_n\right) - \hat{f}(\mathbf {x}_{i,j},t_n) \end{aligned}$$

(11)

which represents the line integral of the funciton g along the particle trajectory reaching $\mathbf {x}_{i,j}$ at the final time $t=t_{n+1}$ starting from time $t=0$.

3.3 More Discussions on the Proposed Lagrangian Interpolation Scheme

Actually we are not the first to use the idea of decomposing the overall flow map into a composition of a sequence of short time flow maps. A similar idea has been proposed in [5] for FTLE computations. In particular, they decompose the overall time-T forward flow map $\varPhi _0^T$ into k time-h forward flow maps $\varPhi _{nh}^{(n+1)h}$ where $n=0,1,\ldots ,k-1$ and $T=kh$. After that they compute these k short time flow maps and store them in the memory and then obtain the the overall flow map by composition and interpolation $\varPhi _0^T=\varPhi _{(k-1)h}^{kh}\circ \cdots \circ \varPhi _h^{2h}\circ \varPhi _0^h$. However, there are several differences between our proposed approach and the approach in [5]. First, in [5] the underlying velocity field is supposed to be defined everywhere in the computational domain, so actually they can choose any Lagrangian approach to compute the required flow maps including those short time flow maps $\varPhi _{nh}^{(n+1)h}$ where h could be much larger than $\varDelta t$. In contrast, our proposed approach is developed especially for discrete velocity data defined only on mesh where we have fixed $h=\varDelta t$ and thus the short time flow maps $\varPhi _{nh}^{(n+1)h}$ actually become one step forward flow maps $\varPhi _{t_n}^{t_{n+1}}$ which can then be simply computed by formula (4) or (5). Furthermore, our proposed approach does not require to store those intermediate flow map data in the memory unlike in [5]. Most importantly, we have extended our approach to efficiently compute the backward flow map while solving the corresponding ODEs forward in time. This makes the numerical implementation more convenient, since the velocity field does not need to be loaded from the terminal time $t=T$ backward in time to the initial time.

3.4 Computational Complexities

We discuss the computational complexities of the two Lagrangian approaches introduced in Sect. 2 and also our proposed approaches here. Suppose that the computational domain $[x_{min},x_{max}]\times [y_{min},y_{max}]$ is discretized uniformly by $M\times N$ mesh points where $M=O(N)$ and the time interval [0, T] is discretized by K points just as we did before. Without loss of generality, we consider the computations of the forward flow map $\varPhi _0^T(\mathbf {x}_{i,j})$ at all mesh points.

In the first way of Lagrangian implementation as introduced in Sect. 2.1, to obtain the forward flow map $\varPhi _0^T(\mathbf {x}_{i,j})$ at a particular mesh point $\mathbf {x}_{i,j}$, one solves ODEs (1) forward from $t=0$ to $t=T$ with the initial condition $\mathbf {x}(0)=\mathbf {x}_{i,j}$. Actually when solving the ODEs, we will first obtain $K-1$ intermediate solutions $\varPhi _0^{t_1}(\mathbf {x}_{i,j})$, $\varPhi _0^{t_2}(\mathbf {x}_{i,j}), \cdots , \varPhi _0^{t_{K-1}}(\mathbf {x}_{i,j})$ in turn and finally the solution $\varPhi _0^T(\mathbf {x}_{i,j})$. If the velocity data is only available at mesh points, however, interpolation on the discrete velocity is required at each time step $t=t_n$ to obtain $\varPhi _0^{t_n}(\mathbf {x}_{i,j})$ which takes O(MN) operations and therefore the computational effort to obtain $\varPhi _0^T(\mathbf {x}_{i,j})$ is O(KMN). As i and j vary in $\{1,2,\ldots ,M\}$ and $\{1,2,\ldots ,N\}$ respectively, the overall computational complexity of computing the forward flow map $\varPhi _0^T(\mathbf {x}_{i,j})$ for all MN mesh points $\mathbf {x}_{i,j}$ is up to $MN\cdot O(KMN) =O(KN^4)$.

In the second way of Lagrangian implementation, the computations of the forward flow map $\varPhi _0^T(\mathbf {x}_{i,j})$ at all mesh points are pushed forward together. In particular, at each time level $t=t_n$ when we have already obtained the intermediate flow map $\varPhi _0^{t_n}(\mathbf {x}_{i,j})$ for all $\mathbf {x}_{i,j}$, we compute the one step forward flow map $\varPhi _{t_n}^{t_{n+1}}(\varPhi _0^{t_n}(\mathbf {x}_{i,j}))$ for all the MN mesh points $\mathbf {x}_{i,j}$ simultaneously by first discretizing the ODEs (1) MN times from $t=t_n$ to $t=t_{n+1}$ with MN different initial conditions $\mathbf {x}(t_n)=\varPhi _0^{t_n}(\mathbf {x}_{i,j})$ and then assigning $\varPhi _{t_n}^{t_{n+1}}(\varPhi _0^{t_n}(\mathbf {x}_{i,j}))=\mathbf {x}(t_{n+1})$. As we discussed in Sect. 2.2, this discretization step is essentially a simple matrix operation in MATLAB which takes only O(MN) operations. Interpolation on $\mathbf {u}(\mathbf {x}_{i,j},t_n)$ is required in the discretization process which fortunately also takes O(MN) operations. As a result, the total complexity at each time step $t=t_n$ remains O(MN) and summing up K time steps together the overall computational complexity of the second way of Lagrangian implementation is only $K\cdot O(MN)=O(KN^2)$ which has greatly reduced the computational complexity in the first way and is optimal in the sense that each grid point is visited for only O(1) time.

Then let’s consider our proposed approach for computing the forward flow map $\varPhi _0^T(\mathbf {x}_{i,j})$ at all mesh points. That is, the computational complexity of Algorithm 1 is studied here. In Algorithm 1, it’s easy to see that the computational complexities of Step 1, Step 2 and Step 4 are negligible compared to that of Step 3 and thus we here only consider Step 3. There are K iterations in Step 3 which exactly correspond to the K time steps and each iteration takes the same order of operations. In particular, at each time step $t=t_n$, the computation of those one step flow maps $\varPhi ^{t_{n+1}}_{t_n}(\mathbf {x}_{i,j})$ for all mesh points $\mathbf {x}_{i,j}$ takes O(MN) operations using (5). After that, interpolating on $\varPhi ^{t_{n+1}}_{t_n}(\mathbf {x}_{i,j})$ to obtain $\varPhi ^{t_{n+1}}_{t_n}(\varPhi _0^{t_n}(\mathbf {x}_{i,j}))$ also takes O(MN) operations. Therefore, in Step 3 each iteration takes $O(MN)+O(MN)=O(MN)$ operations. Summing up the operations in all the K iterations gives the overall computational complexity $K\cdot O(MN)=O(KN^2)$ which is same as that of the second way of Lagrangian implementation.

4 Analysis of the Proposed Approaches

In this section, we study the overall accuracy of numerical solutions computed based on the proposed Lagrangian interpolation approaches. First we would like to list the following two lemmas which extend the corresponding results in [33] from autonomous velocity fields to time-dependent flows as discussed in the current work. These extensions are straight-forward and we omit the proof.

Lemma 1

If the velocity field $\mathbf {u}(\mathbf {x},t)$ is sufficiently smooth on the computational domain M with the Lipschitz constant L, then we have

$$\begin{aligned} \left| \varPhi _0^t(\mathbf {x}_1)-\varPhi _0^t(\mathbf {x}_2)\right| \le e^{Lt}|\mathbf {x}_1-\mathbf {x}_2|,\,\,\,\forall \mathbf {x}_1, \mathbf {x}_2\in M. \end{aligned}$$

Lemma 2

If the velocity field $\mathbf {u}(\mathbf {x},t)$ is sufficiently smooth on the computational domain M with the Lipschitz constant L, then there exists a constant $C_s$ for each $s\ge 2$ such that for any multi-index $\gamma$ satisfying $|\gamma |=s$, we have

$$\begin{aligned} \left| \partial ^{\gamma }\varPhi _0^t(\mathbf {x})\right| \le C_s \, t \, e^{(2s-1)Lt} ,\,\,\, \forall t>0 \,, \forall \mathbf {x}\in M. \end{aligned}$$

Upon these two lemmas, we can now state and prove the following theorems on the accuracy of the proposed Lagrangian interpolation algorithms.

Theorem 1

If $\varPhi _{t_n}^{t_{n+1}}$ is computed using (5) and the cubic spline scheme is used for all required interpolation, then the forward flow map $\varPhi _0^T$ computed with the proposed approach in Sect. 3.1 is second order accurate in both spatial and temporal directions.

Proof

To have a clear explanation, we here use $\tilde{\cdot }$ to denote the numerical approximation and ${\mathcal {I}}$ to denote the interpolation operator. Formula (5) can then be rewritten as

$$\begin{aligned} \tilde{\varPhi }_{t_n}^{t_{n+1}}(\mathbf {x}_{i,j})= \mathbf {x}_{i,j} +\frac{\varDelta t}{2} \left[ \mathbf {u}(\mathbf {x}_{i,j},t_n)+{\mathcal {I}}\mathbf {u}(\mathbf {x}_{i,j}+\varDelta t\mathbf {u}(\mathbf {x}_{i,j},t_n),t_{n+1}) \right] . \end{aligned}$$

As a result,

$$\begin{aligned}&\left| \tilde{\varPhi }_{t_n}^{t_{n+1}}(\mathbf {x}_{i,j})-\varPhi _{t_n}^{t_{n+1}}(\mathbf {x}_{i,j})\right| \le C_1\varDelta t^3+\frac{\varDelta t}{2}|\mathbf {u}(\mathbf {x}_{i,j}+\varDelta t\mathbf {u}(\mathbf {x}_{i,j},t_n),t_{n+1})\\&\quad -{\mathcal {I}}\mathbf {u}(\mathbf {x}_{i,j}+\varDelta t\mathbf {u}(\mathbf {x}_{i,j},t_n),t_{n+1})|\,. \end{aligned}$$

The term $C_1\varDelta t^3$ is due to the discretization to the ODEs and the constant $C_1$ depends only on the velocity field $\mathbf {u}$. Since the cubic spline interpolation is second order accurate,

$$\begin{aligned}&|\mathbf {u}(\mathbf {x}_{i,j}+\varDelta t\mathbf {u}(\mathbf {x}_{i,j},t_n),t_{n+1})-{\mathcal {I}}\mathbf {u}(\mathbf {x}_{i,j}+\varDelta t\mathbf {u}(\mathbf {x}_{i,j},t_n),t_{n+1})|\\&\quad \le C_2^{\prime }\varDelta x^2\max \limits _{|\gamma |=2}\sup \limits _{\mathbf {x}_{i,j}}|\partial ^{\gamma }\mathbf {u}(\mathbf {x}_{i,j},t_{n+1})|. \end{aligned}$$

Since the velocity field $\mathbf {u}$ is supposed to be smooth enough, there exists a positive constant $C_2^{\prime \prime }$ such that

$$\begin{aligned} \max \limits _{|\gamma |=2}\sup \limits _{\mathbf {x}_{i,j}}|\partial ^{\gamma }\mathbf {u}(\mathbf {x}_{i,j},t_{n+1})|\le C_2^{\prime \prime } \end{aligned}$$

and therefore

$$\begin{aligned} \left| \tilde{\varPhi }_{t_n}^{t_{n+1}}(\mathbf {x}_{i,j})-\varPhi _{t_n}^{t_{n+1}}(\mathbf {x}_{i,j})\right| \le C_1\varDelta t^3+C_2\varDelta t\varDelta x^2 \end{aligned}$$

(12)

where $C_2\triangleq \frac{1}{2}C_2^{\prime }C_2^{\prime \prime }$.

Then $\varPhi _0^{t_{n+1}}$ is approximated by

$$\begin{aligned} \tilde{\varPhi }_0^{t_{n+1}}(\mathbf {x}_{i,j})={\mathcal {I}}\tilde{\varPhi }_{t_n}^{t_{n+1}}\circ \tilde{\varPhi }_0^{t_n}(\mathbf {x}_{i,j}) \end{aligned}$$

with the error computed as

$$\begin{aligned} \left| \tilde{\varPhi }_0^{t_{n+1}}(\mathbf {x}_{i,j})-\varPhi _0^{t_{n+1}}(\mathbf {x}_{i,j})\right|\le & {} \left| {\mathcal {I}}\tilde{\varPhi }_{t_n}^{t_{n+1}}\left( \tilde{\varPhi }_0^{t_n}(\mathbf {x}_{i,j})\right) -{\mathcal {I}}\varPhi _{t_n}^{t_{n+1}}\left( \tilde{\varPhi }_0^{t_n}(\mathbf {x}_{i,j})\right) \right| \\&+\left| {\mathcal {I}}\varPhi _{t_n}^{t_{n+1}}\left( \tilde{\varPhi }_0^{t_n}(\mathbf {x}_{i,j})\right) -\varPhi _{t_n}^{t_{n+1}}\left( \tilde{\varPhi }_0^{t_n}(\mathbf {x}_{i,j})\right) \right| \\&+\left| \varPhi _{t_n}^{t_{n+1}}\left( \tilde{\varPhi }_0^{t_n}(\mathbf {x}_{i,j})\right) -\varPhi _{t_n}^{t_{n+1}}\left( \varPhi _0^{t_n}(\mathbf {x}_{i,j})\right) \right| \\\triangleq & {} I_1+I_2+I_3. \end{aligned}$$

Suppose the norm of the interpolation operator is $N_I$, then $I_1$ is bounded by

$$\begin{aligned} I_1\le N_I\cdot \max \limits _{\mathbf {x}_{i,j}}\left| \tilde{\varPhi }_{t_n}^{t_{n+1}}(\mathbf {x}_{i,j})-\varPhi _{t_n}^{t_{n+1}}(\mathbf {x}_{i,j})\right| \le C_1N_I\varDelta t^3+C_2N_I\varDelta t\varDelta x^2. \end{aligned}$$

According to Lemma 2, there exists a constant $C_3^{\prime }$ such that

$$\begin{aligned} \max \limits _{|\gamma |=2}\sup \limits _{\mathbf {x}_{i,j}}\left| \partial ^{\gamma }\varPhi _{t_n}^{t_{n+1}}(\mathbf {x}_{i,j})\right| \le C_{3}^{\prime }\cdot e^{3 L\varDelta t}\cdot \varDelta t\,. \end{aligned}$$

Then $I_2$ is estimated as

$$\begin{aligned} I_2\le C_3^{\prime \prime }\varDelta x^2\max \limits _{|\gamma |=2}\sup \limits _{\mathbf {x}_{i,j}}\left| \partial ^{\gamma }\varPhi _{t_n}^{t_{n+1}}(\mathbf {x}_{i,j})\right| \le C_3 e^{3L\varDelta t}\cdot \varDelta t\varDelta x^2 \end{aligned}$$

where $C_3\triangleq C_3^{\prime }C_3^{\prime \prime }$ is a constant. Finally, according to Lemma 1, $I_3$ is bounded by $I_3\le e^{L\varDelta t}|\tilde{\varPhi }_0^{t_n}(\mathbf {x}_{i,j})-\varPhi _0^{t_n}(\mathbf {x}_{i,j})|$. Supposing $E^{k}\triangleq \max \limits _{\mathbf {x}_{i,j}}|\tilde{\varPhi }_0^{t_k}(\mathbf {x}_{i,j})-\varPhi _0^{t_k}(\mathbf {x}_{i,j})|$, we have $E^{n+1}\le e^{L\varDelta t}E^n+C_4\varDelta t(\varDelta x^2+\varDelta t^2)$ where $C_4$ is a constant depending on $C_1$, $C_2$, $C_3$, $N_I$, L. With the use of recurrence, we can obtain

$$\begin{aligned} E^{n+1}\le e^{nL\varDelta t}E^1+C_4\cdot n e^{nL\varDelta t}\varDelta t(\varDelta t^2+\varDelta x^2)\le e^{Lt}E^1+C_4\cdot t e^{Lt}\cdot (\varDelta t^2+\varDelta x^2). \end{aligned}$$

On the other hand, it is obvious that $E^1\le C_1\varDelta t^3+C_2\varDelta t\varDelta x^2$. Omitting the higher order terms, we finally have $E^{n+1}\le C_5(\varDelta t^2+\varDelta x^2)$ where $C_5$ only depends on $C_4$ and T. $\square$

Theorem 2

Suppose that $g(\mathbf {x})$ is a sufficiently smooth function, $\hat{f}$ is approximated using (7) and the cubic spline scheme is used for all required interpolation. Then the line integral of $g(\mathbf {x})$ computed with (8) has second order accuracy in both the temporal and spatial directions.

Proof

Again we use ${\mathcal {I}}$ to denote the interpolation operator and $\tilde{\cdot }$ to denote the numerical approximation. We first rewrite formula (8) as

$$\begin{aligned} \tilde{f}(\mathbf {x}_{i,j},t_{n+1}) = \tilde{f}(\mathbf {x}_{i,j},t_{n}) + {\mathcal {I}}\hat{f}_a\left( \tilde{\varPhi }_0^{t_n}(\mathbf {x}_{i,j}), t_{n+1}\right) \end{aligned}$$

where $\hat{f}_a$ is the numerical solution of $\hat{f}$. According to (7), $\hat{f}(\mathbf {x}_{i,j},t_{n+1})$ is approximated as

$$\begin{aligned} \hat{f}_a(\mathbf {x}_{i,j},t_{n+1})= \frac{\varDelta t}{2} [g(\mathbf {x}_{i,j})+g(\mathbf {x}_{i,j}+\varDelta t\mathbf {u}(\mathbf {x}_{i,j},t_n))]\,. \end{aligned}$$

Since the function g is explicitly given, no interpolation is required in the above formula and the error of $\hat{f}_a(\mathbf {x}_{i,j},t_{n+1})$ only comes from the discretization to the ODE (6). As a result, we have

$$\begin{aligned} |\hat{f}_a(\mathbf {x}_{i,j},t_{n+1})-\hat{f}(\mathbf {x}_{i,j},t_{n+1})|\le C_1\varDelta t^3\,. \end{aligned}$$

Let $E^{n}(\mathbf {x}_{i,j})\triangleq |\tilde{f}(\mathbf {x}_{i,j},t_{n})-f(\mathbf {x}_{i,j},t_{n})|$, then

$$\begin{aligned} E^{n+1}(\mathbf {x}_{i,j})\le & {} E^n(\mathbf {x}_{i,j}) + \left| {\mathcal {I}}\hat{f}_a\left( \tilde{\varPhi }_0^{t_n}(\mathbf {x}_{i,j}), t_{n+1}\right) -\hat{f}\left( \varPhi _0^{t_n}(\mathbf {x}_{i,j}), t_{n+1}\right) \right| \\\le & {} E^n(\mathbf {x}_{i,j})+\left| {\mathcal {I}}\hat{f}_a\left( \tilde{\varPhi }_0^{t_n}(\mathbf {x}_{i,j}), t_{n+1}\right) -{\mathcal {I}}\hat{f}\left( \tilde{\varPhi }_0^{t_n}(\mathbf {x}_{i,j}), t_{n+1}\right) \right| \\&+\left| {\mathcal {I}}\hat{f}\left( \tilde{\varPhi }_0^{t_n}(\mathbf {x}_{i,j}), t_{n+1}\right) -\hat{f}\left( \tilde{\varPhi }_0^{t_n}(\mathbf {x}_{i,j}), t_{n+1}\right) \right| \\&+\left| \hat{f}\left( \tilde{\varPhi }_0^{t_n}(\mathbf {x}_{i,j}), t_{n+1}\right) -\hat{f}\left( \varPhi _0^{t_n}(\mathbf {x}_{i,j}), t_{n+1}\right) \right| \\\triangleq & {} E^n(\mathbf {x}_{i,j})+I_1+I_2+I_3\, . \end{aligned}$$

Suppose the interpolation has the norm $N_I$, then $I_1$ is estimated as $I_1\le N_I\cdot \max \limits _{\mathbf {x}_{i,j}}|\hat{f}_a(\mathbf {x}_{i,j}, t_{n+1})-\hat{f}(\mathbf {x}_{i,j}, t_{n+1})|\le C_1N_I\varDelta t^3$. On the other hand, there exist a constant $C_2^{\prime }$ such that

$$\begin{aligned} I_2\le C_2^{\prime }\cdot \varDelta x^2\max \limits _{|\gamma |=2}\sup \limits _{\mathbf {x}_{i,j}}|\partial ^{\gamma }\hat{f}(\mathbf {x}_{i,j},t_{n+1})| \, . \end{aligned}$$

Differentiating both sides of $\hat{f}(\mathbf {x},t_{n+1})=\int _{t_n}^{t_{n+1}}g(\varPhi _{t_n}^{\tau }(\mathbf {x}))d\tau$ to obtain $\partial ^{\gamma }\hat{f}=\int _{t_n}^{t_{n+1}}\partial ^{\gamma }g(\varPhi _{t_n}^{\tau }(\mathbf {x}))d\tau$ where

$$\begin{aligned} \partial ^{\gamma }g\left( \varPhi _{t_n}^{\tau }(\mathbf {x})\right) =\sum \limits _{p=1}^2\sum \limits _{(\gamma _1,\ldots ,\gamma _p):\gamma =\sum \limits _{l=1}^p\gamma _l}\left( \sum \limits _{k_1,\ldots ,k_p=1}^d g_{k_1\ldots k_p}\prod \limits _{j=1}^p\partial ^{\gamma _j}\phi ^{k_j}\right) \end{aligned}$$

in which $g_{k_1\ldots k_p}\triangleq \frac{\partial ^p g}{\partial x_{k_1}\ldots \partial x_{k_p}}$ and $\phi ^i$ represents the i-th component of $\varPhi _{t_n}^{\tau }(\mathbf {x})$. From Lemma 2, we obtain

$$\begin{aligned} \left| \prod \limits _{j=1}^p\partial ^{\gamma _j}\phi ^{k_j}\right| \le C_2 \end{aligned}$$

for any $p\ge 1$ where $C_2$ is a constant depending only on the Lipschitz constant L. Since all these terms are finite, we then have $|\partial ^{\gamma }g(\varPhi _{t_n}^{\tau }(\mathbf {x}))|\le C_3^{\prime }$ where $C_3^{\prime }$ depends on d, $C_2$ and $\max \limits _{|\gamma |\le 2}|\partial ^{\gamma }g|_{L_{\infty }}$. As a result, $|\partial ^{\gamma }\hat{f}|\le C_3^{\prime }\cdot \varDelta t$ and hence $I_2\le C_3\cdot \varDelta x^2\cdot \varDelta t$ where $C_3=C_2^{\prime }C_3^{\prime }$.

Furthermore,

$$\begin{aligned} I_3\le \sup \limits _{\mathbf {x}\in M}|\nabla \hat{f}(\mathbf {x},t_{n+1})|\left| \tilde{\varPhi }_0^{t_n} (\mathbf {x}_{i,j})-\varPhi _0^{t_n}(\mathbf {x}_{i,j})\right| . \end{aligned}$$

For all $|\gamma |\le 2$, it has been shown that $|\partial ^{\gamma }\hat{f}|\le C_3^{\prime }\cdot \varDelta t$, then $|\nabla \hat{f}(\mathbf {x},t_{n+1})|\le C_4\varDelta t$ for all $\mathbf {x}\in M$. According to Theorem 1, we immediately have $|\tilde{\varPhi }_0^{t_n}(\mathbf {x}_{i,j})-\varPhi _0^{t_n}(\mathbf {x}_{i,j})|\le C_5(\varDelta t^2+\varDelta x^2)$ for certain constant $C_5$. In summary, $I_3\le C_4\cdot C_5\cdot (\varDelta t^3+\varDelta t\varDelta x^2)$.

Denote $E^k=\max \limits _{\mathbf {x}_{i,j}}E^k(\mathbf {x}_{i,j})$, then

$$\begin{aligned} E^{n+1}\le & {} E^n+C_1N_I\varDelta t^3+C_3\cdot \varDelta x^2\cdot \varDelta t+C_4 C_5(\varDelta t^3+\varDelta t\varDelta x^2)\\\le & {} E^n+C\varDelta t(\varDelta t^2+\varDelta x^2) \end{aligned}$$

where C is a constant. Finally, we have $E^{n+1}\le E^0+C\cdot t\cdot (\varDelta t^2+\varDelta x^2)\le C\cdot T\cdot (\varDelta t^2+\varDelta x^2)$ by recursion. $\square$

Theorem 3

Given a sufficiently smooth integrable function $g(\mathbf {x})$. If $\varPhi _{t_{n+1}}^{t_n}$ and $\hat{f}$ are estimated using (9) and (10) respectively and the cubic spline scheme is used for all required interpolation, then the line integral of $g(\mathbf {x})$ computed with (11) has second order accuracy in both the temporal and the spatial directions.

Proof

The proof is quite similar to that of Theorem 2 and we omit it here. $\square$

Remark

In the proof of all the above theorems, the cubic spline scheme is used for all required interpolation. Since the underlying mesh points have been made equally spaced, the cubic spline scheme is stable in terms of $L^{\infty }$ norm [4]. Indeed any other interpolation scheme can be used as long as its boundedness is guaranteed. However, we are definitely not saying that the higher the order of the interpolation scheme, the more accurate the numerical solutions. Actually some higher order interpolation schemes are not even guaranteed to be bounded, e.g. Runge’s counterexample [29], which will eventually lead to bad numerical results.

Fig. 1

(Section 5.1) The proposed Lagrangian interpolation approach to Problem B. a The line integral f of the function $g(x,y)=\cos (12\pi x)\cos (2\pi y)$ along the particle trajectory starting from $t=0$ reaching the location $\mathbf {x}_{i,j}$ at $t=10$. b The $L_2$ error of the numerical solution f using TVD-RK2 with different $\varDelta x$’s and $\varDelta t$’s

Fig. 2

(Section 5.1) The proposed new Lagrangian interpolation approach to Problem F. The solutions of $\phi$ and $\psi$ are shown respectively in (a, b) computed using TVD-RK2 and cubic spline

Fig. 3

(Section 5.1) The solution to Problem F. The solution of the line integral of the function $g(x,y)=\cos (12\pi x)\cos (2\pi y)$ along the particle trajectory starting from $\mathbf {x}_{i,j}$ at $t=0$ to $t=10$. a The solution of our proposed Lagrangian approach using TVD-RK2 and cubic spline interpolation. b The Eulerian approach where the TVDRK2-WENO5 scheme is used. c The second way of Lagrangian implementation where the ode solver is TVD-RK2

Fig. 4

(Section 5.1) The proposed Lagrangian approach to Problem F. The $L_2$ errors of $\phi$, $\psi$ and f using TVD-RK2 and cubic spline interpolation where a $\varDelta t=0.0001$ and $\varDelta x$ varies from 1 / 32 to 1 / 512, b $\varDelta x=1/512$ and $\varDelta t$ varies from 1 / 5 to 1 / 80. The red solid line and the blue dash dot line correspond to the errors of $\phi$ and $\psi$ respectively while the error of f is represented by the black dashed line. Black solid lines are the reference lines with the slope 2 (Color figure online)

Fig. 5

(Section 5.1) The $L_2$ error of the forward flow map $\varPhi _0^{10}=(\phi ,\psi )$ computed using the two old ways of Lagrangian implementation and the proposed approach. a The $L_2$ error of $\phi$, b the $L_2$ error of $\psi$. Red solid lines: the first way of Lagrangian implementation. Blue dashed lines: the second way of Lagrangian implementation. Black dash dotted lines: the proposed approach (Color figure online)

5 Numerical Examples

5.1 The Double gyre Flow

As the first example, we consider a periodically varying double-gyre flow which is taken from [31]. The stream function of this flow is given by $\psi (x,y,t)=A \sin [ \pi k(x,t) ] \sin (\pi y)$, where

$$\begin{aligned} k(x,t)= & {} a(t) x^2 + b(t) x, \\ a(t)= & {} \epsilon \sin (\omega t), \\ b(t)= & {} 1- 2\epsilon \sin (\omega t). \end{aligned}$$

Here we follow [31] and use $A=0.1$, $\omega =2\pi /10$.

According to (11), we first compute the line integral of the integrable function $g(x,y)=\cos (12\pi x)\cos (2\pi y)$ along each particle trajectory starting from the initial time $t=0$ arriving at the mesh point $\mathbf {x}_{i,j}$ at the time $t=10$. The result is shown in Fig. 1a with $\varDelta x=\varDelta y=1/256$ and the cubic spline method is used wherever the interpolation is required. Furthermore, we have used (10) and (9) to compute $\hat{f}$ and $\varPhi _{t_{n+1}}^{t_n}$, respectively. That is, $\hat{f}$ and $\varPhi _{t_{n+1}}^{t_n}$ are both estimated by solving corresponding ODEs using the TVD-RK2 scheme. According to Theorem 3, the solution of f must be second order accurate with respect to both $\varDelta x$ and $\varDelta t$. To verify this point, in Fig. 1b we show the $L_2$ errors of f computed using different $\varDelta x$’s varying from 1 / 32 down to 1 / 512 while keeping $\varDelta t/\varDelta x$ constant. Since the exact solution for this flow is not available, the exact solution is computed using the traditional Lagrangian ray tracing method with a very small time step. We can see that the numerical solution of f computed using the proposed approach has approximately second order accuracy satisfying Theorem 3. It has to be emphasized that we have only used velocity data on the mesh in the whole implementation and the corresponding ODEs are solved forward in time from $t=0$ to $t=10$ which is impossible in the traditional Lagrangian approach.

We then use the proposed new Lagrangian approach in Sect. 3.1 to compute the forward flow map $\varPhi _0^{10}(\mathbf {x})$. Suppose $\phi$ and $\psi$ are the two components of $\varPhi _0^{10}(\mathbf {x})$ and then the solutions of $\phi$ and $\psi$ are shown respectively in Fig. 2a, b using $\varDelta x=\varDelta y=1/256$. Also we use the proposed approach to compute the line integral $f(\mathbf {x},10)$ of the function $g(x,y)=\cos (12\pi x)\cos (2\pi y)$ along the particle trajectory starting from $\mathbf {x}_{i,j}$ at the initial time $t=0$ reaching the final time $t=10$. In particular, $f(\mathbf {x},10)$ is computed using (8) in which both $\hat{f}$ and $\varPhi _{t_n}^{t_{n+1}}$ are approximated by solving the corresponding ODEs using the TVD-RK2 scheme. That is, we use (7) and (5) to compute $\hat{f}$ and $\varPhi _{t_n}^{t_{n+1}}$, respectively. Furthermore, all required interpolation is done with the cubic spline method with $\varDelta x=\varDelta y=1/256$. The solution of $f(\mathbf {x},10)$ is plotted in Fig. 3a. As a comparison, Fig. 3b, c show the corresponding line integral of the same g(x, y) computed using the Eulerian approach and the second way of Lagrangian implementation as in Sect. 2.2, respectively.

To verify Theorems 1 and 2, we compute the errors of the forward flow map and the corresponding line integral with respect to both $\varDelta t$ and $\varDelta x$. In particular, we first fix $\varDelta t=0.0001$ and let $\varDelta x$ vary from 1 / 32 to 1 / 512. Figure 4 (a) shows the corresponding $L_2$ errors of the two components $\phi$ and $\psi$ of the flow map and also the error of the line integral f. It can be easily seen that all solutions have approximately second order accuracy with respect to $\varDelta x$ which satisfies the claim in Theorems 1 and 2. Then in Fig. 4b we plot the $L_2$ errors of two components of the flow map and the line integral f using different $\varDelta t$’s ranging from 1 / 5 to 1 / 80 while fixing $\varDelta x=1/512$. As can be seen, those solutions are also approximately second order accurate with respect to $\varDelta t$ satisfying Theorems 1 and 2.

We then compare the accuracy and efficiency of the two old ways of Lagrangian implementation introduced in Sect. 2 and our proposed approach based on the velocity data given only at mesh points. Figure 5a, b show the $L_2$ errors of two components $\phi$ and $\psi$ of the flow map $\varPhi _0^{10}(\mathbf {x})$, respectively. In this figure, red solid lines, blue dashed lines and black dash dotted lines correspond to the solutions of the first way of Lagrangian implementation, the second way of Lagrangian implementation and our proposed approach, respectively. In the implementation, we always keep $\varDelta t/\varDelta x$ fixed and let $\varDelta x$ vary from 1 / 32 to 1 / 512 in the second way of Lagrangian implementation and in our proposed approach and from 1 / 32 to 1 / 256 in the first way of Lagrangian implementation. As can be seen, all approaches are second order accurate and they have almost the same accuracy. Table 1 gives a comparison of all the three approaches in the computational (CPU) time where we have also kept $\varDelta t/\varDelta x$ fixed and let $\varDelta x$ vary. We do not give the solution of the first way of Lagrangian implementation for $\varDelta x=1/512$ since it might cost nearly one month. The numerical results are obtained using a laptop computer with a 2.5 GHz Intel core i7 processor. From this table we can see that our proposed approach is a little more efficient than the second way of Lagrangian implementation and much more efficient than the first way of Lagrangian implementation.

Table 1 The computational time of our proposed approach and the two old ways of Lagrangian implementation with different $\varDelta x$’s and fixed $\varDelta t/\varDelta x$

In all the above implementations of our proposed Lagrangian interpolation approach, we have used the TVD-RK2 scheme to discretize the corresponding ODEs and the cubic spline method wherever interpolation is required. Finally, Fig. 6 compares the errors of the proposed Lagrangian approach for Problem F using different discretization schemes and different interpolation methods. The bi-linear interpolation and cubic spline interpolation are used in Fig. 6a, b, respectively. In each of these two figures, we give the errors in the TVD-RK2 method, the TVD-RK3 method and also the first order Euler method. In the implementation, we still keep $\varDelta t/\varDelta x$ fixed and let $\varDelta x$ vary from 1 / 32 to 1 / 512. With the simple Euler method, the solution always shows approximately first order accuracy no matter the bi-linear interpolation or the cubic spline interpolation is used, as plotted using the red solid lines. This is because that the Euler method has first order accuracy in the temporal direction. When the TVD-RK2 scheme is used, we can observe second order accuracy in the solutions as plotted in blue dash dot lines which verifies the claim in Theorem 2. With black dash lines, we finally show the errors in the solution using the TVD-RK3 scheme. If the bi-linear interpolation is used, the errors in the solution are dominated by the interpolation error which is second order accurate in space as plotted in Fig. 6a. With the cubic spline interpolation instead, the errors show a convergence of higher than third order, as shown in Fig. 6b.

5.2 An Application to Real Dataset

As discussed in Sect. 3, only discrete velocity data on the mesh is required in our proposed Lagrangian interpolation approaches for computing the forward (backward) flow map and the corresponding line integrals of certain given function along particle trajectories. In this subsection, to demonstrate its effectiveness, we focus on the Ocean Surface Current Analyses Real-time (OSCAR) dataset in which velocity data is only available at discrete locations. The OSCAR data were obtained from JPL Physical Oceanography DAAC and developed by ESR. It covers $0^{\circ }$ to $360^{\circ }$ longitude and $-80^{\circ }$ to $80^{\circ }$ latitude. The resolution is $1/3^{\circ }$ in each spatial direction and about 5 days in the temporal direction.

Fig. 6

(Section 5.1) The $L_2$ error of the line integral f computed with different $\varDelta x$’s using the proposed Lagrangian approach where the interpolation method is a bi-linear interpolation, b cubic spline. Red solid line: first order Euler method. Blue dash dot line: TVD-RK2. Black dash line: TVD-RK3. Black solid lines are the reference lines with the slopes explicitly given (Color figure online)

We have chosen an ocean region near the Line Islands as the computational region, which is enclosed by $S17^{\circ }$ to $N8^{\circ }$ latitude and $E180^{\circ }$ to $E230^{\circ }$ longitude. For a better visualization, we first interpolate the velocity data to obtain a finer resolution of $1/3^{\circ }$ in each spatial direction and 0.25 days in the temporal direction. Numerically, we simulate the ocean surface current in the first 50 days in year 2015 using Algorithm 2 proposed in Sect. 3.1. In particular, we compute the line integral of $g(x,y)=\cos (0.5\pi x)\cos (0.1\pi y)$ along the particle trajectories taking off from the point (x, y) at the time $t=0$ arriving at the time $t=50$, as shown in Fig. 7.

Fig. 7

(Section 5.2) The solution of the line integral of $f(x,y)=\cos (0.5\pi x)\cos (0.1\pi y)$ from $t=0$ to $t=50$ along particle trajectories of the OSCAR dataset

5.3 Application to the Coherent ergodic Partition

The so-called coherent ergodic partition was proposed in [34] in order to extract invariant sets of a given time-dependent dynamical system in the extended phase space. To obtain this kind of partition, one needs to compute the forward flow map $\varPhi _0^T(\mathbf {x})$ for a sufficiently large $T>0$ and also the time averages of given observables along particle trajectories. In particular, the time average $g^{\star }$ of the observable $g(\mathbf {x})$ is computed as

$$\begin{aligned} g^{\star }(\mathbf {x})=\lim \limits _{T\rightarrow \infty }\frac{1}{T}\int _0^T\left( g\circ \varPhi _0^\tau \right) (\mathbf {x})d\tau . \end{aligned}$$

We here use Algorithm 2 proposed in Sect. 3.1 to approximate these quantities and then give the coherent ergodic partition. We still consider the double-gyre flow and use $A=0.1$, $\epsilon =0.1$, $\omega =2\pi /10$. The domain $[0,2]\times [0,1]$ is discretized such that $\varDelta x=\varDelta y=1/256$. First we solve the corresponding ODEs from $t=t_0=0$ to $t=T_m=10$ using Algorithm 2 to obtain $\varPhi _0^{T_m}$ and also $\int _0^{T_m}(g\circ \varPhi _0^\tau )(\mathbf {x})d\tau$. Then we iterate the flow map and also the line integral 5 times to obtain $\varPhi _0^T(\mathbf {x})$ and $\int _0^{T}(g\circ \varPhi _0^\tau )(\mathbf {x})d\tau$ for $T=T_m\cdot 2^5=320$.

Fig. 8

(Section 5.3) a The numerical approximation to the time average of the observable $g(\mathbf {x})=\cos (2\pi y)$ computed using Algorithm 2. b The Ergodic partition at the time $t=0$ computed using the four observables ($g_1$, $g_2$, $g_3$, $g_4$) and $T=320$

Figure 8a shows the quantity $\frac{1}{T}\int _0^{T}(g\circ \varPhi _0^\tau )(\mathbf {x})d\tau$ for $T=320$ using $g(\mathbf {x})=\cos (2\pi y)$. The whole computational domain, as can be seen, is roughly partitioned into five regions. They are two circular regions centered near (0.5, 0.5) and (1.5, 0.5), two arrow regions slightly below and slightly above these two circular regions respectively, and also the background region which contains lots of isolated dots. Within each of these regions, the values of $g^{\star }$ are similar. Then we use more test functions to better distinguish these regions. In particular, with the use of the four independent functions $g_1=\cos (2\pi y)$, $g_2=\cos (12\pi x)\cos (2\pi y)$, $g_3=\sin (2\pi y)$, and $g_4=\frac{x}{2}$, the domain is more clearly classified into five clusters as shown in Fig. 8b. These results match very well with those in [34].

5.4 Application to Geometrical Optics

Geometrical optics is an important class of asymptotic approximation to high frequency wave propagation. In the high frequency regime when $\omega \rightarrow \infty$, we approximate the phase function by the eikonal equation $|\nabla T| = \frac{1}{c}$. To obtain the multivalued solution to the traveltime field T in two dimensions, we reformulate the problem in the phase space by solving the eikonal equation using the Lagrangian formulation by the method of characteristics [8],

$$\begin{aligned} \frac{dx}{dt}= & {} c\sin \theta ,\nonumber \\ \frac{dy}{dt}= & {} c\cos \theta ,\nonumber \\ \frac{d\theta }{dt}= & {} c_y\sin \theta -c_x\cos \theta . \end{aligned}$$

(13)

For simplicity, we transform the above system to a reduced system by using depth as the running parameter

$$\begin{aligned} \frac{dx}{dy}= & {} \tan \theta , \nonumber \\ \frac{d\theta }{dy}= & {} \frac{c_y}{c}\tan \theta -\frac{c_x}{c}. \end{aligned}$$

(14)

As an important application of the coherent ergodic partition, we have studied the ray spreading in [34]. The emitted rays were finally coupled into several clusters and rays would have close behaviors if they are classified into the same cluster. The result could help us refine the angle space adaptively near the boundary of each cluster.

In particular, we used $c=1+0.2\sin (0.5\pi y)\sin (3\pi (x+0.55))$ and $g=\sin x$ in [34]. Then we solved the reduced system (14) from $y=0$ to $y=8$ and computed the line integral using

$$\begin{aligned} f(x,\theta ,8)\triangleq \int ^{8}_{0}g\circ \varPhi _0^y(x,\theta )\,dy \end{aligned}$$

to obtain the coherent ergodic clusters of the system.

In this subsection, we propose to recompute the solution of $f(x,\theta ,8)$ using Algorithm 2 with another test function $g=\sin (5x)$. From Fig. 9a, we can see that the structures of our new solution is quite similar to that in [34] which implies that we have obtained almost the same clusters with [34]. For example, on the $x=0$ cross section we can observe three subintervals $\theta \in [0.1375,0.21875]$, [0.0625, 0.125] and $[-0.4625,-0.18125]$ which share the same color. This implies that rays emitting from the origin with departure angles in the same subinterval may travel as a patch, as shown in Fig. 9. This observation matches exactly with that in [34].

Fig. 9

(Section 5.4) The departure angle subintervals with the same color are shown in a [0.1375, 0.21875], c [0.0625, 0.125] and e $[-0.4625,-0.18125]$; b, d, f Rays emitting from the origin with the departure angles in (a), (c) and (e), respectively