Neural Integral Equations

We hereby report additional figures and tables that complement the results given in Section 3.

The model DeepOnet + UNET in Table 2 is implemented similarly to [DAK23].

Neural ODEs (NODEs) can be slow and have poor scalability [KBJD20]. As such, several methods have been introduced to improve their performance [KBJD20, KCL21, DKM⁺20, PMY⁺20, PMSR21]. Despite these improvements, NODE is still significantly slower than discrete methods such as LSTMs. We hypothesize that (A)NIE has significantly better scalability than NODE, comparable to fast but discrete LSTMs, despite being a continuous model. To test this, we compare NIE and ANIE to the latest optimized version of (latent) NODE [RCD19] and to LSTM on three different dynamical systems: Lotka-Volterra equations, Lorenz system, and IE-generated 2D spirals (see Appendix F for data generation details). During training, models were initialized with the first half of the data and were tasked to predict the second half. Training speed are reported in Table 5. While all models achieve comparable (good) fits to the data, we find that ANIE outperforms all models in 2 out of the 3 datasets in terms of speed. Furthermore, ANIE has better MSE compared to all other models.

For most deep learning models, including NODE, finding numerically stable solutions usually requires an extensive hyperparameter search. Since IE solvers are known to be more stable than ODE solvers, we hypothesize that (A)NIE is less sensitive to hyperparameter changes. To test this, we quantify model fit, for the Lotka-Volterra dynamical system, as a function of two different hyperparameters: learning rate and L2 norm weight regularization. We perform this experiment for three different models: LSTM, Latent NODE, and ANIE. As shown in Figure 10, we find that ANIE generally has lower validation error as well as more consistent errors across hyperparameter values, compared to LSTM and NODE, therefore validating our hypothesis.

To further test the ability of ANIE in modeling non-local systems, we benchmark ANIE, NODE, and LSTM on a dataset of 2D spirals generated by integral equations. This data consists of 500 2D curves of 100 time points each. The data was split in half for training and testing. During training, the first 20 points were given as initial condition and the models were tasked to predict the full 100 point dynamics. Details on the data generation are described in Appendix F. For ANIE, the initialization is given via the free function $f$, which assumes the values of the first 20 points and sets the remaining 80 points to be equal to the value of the 20th point. For NODE, the initialization is given as the reverse RNN on the first 20 points, which outputs a distribution corresponding to the first time point (see [CRBD18] for more details on their Latent ODE experiments). For LSTM, we input the data in segments of 20 points to predict the consecutive point of the sequence. The process is repeated with the output of the previous step until all points of the curve are predicted. During inference, we test the models’ performance on never before seen initial conditions. Table 6 shows the correlation between the ground truth curve and the model predictions. Figure 9 shows the correlation coefficients for the 500 curves. In summary, ANIE significantly outperforms the other tested methods in predicting IE generated non-local dynamics.

We now consider the convergence of the solver to a solution of an IE for a trained model. Our experiment here considers a model that has been trained with a number of iteration, and we explore whether the solver iterations converge to a solution at the end of the training. These results show that the model learns to converge to a solution of Equation 5 within the iterations that are fixed during training. They show that a fixed point for IE is obtained when outputting a prediction.

Figure 12 and Figure 3 show the convergence error (i.e. the value $||y_{n+1}-y_{n}||$), and the guesses produced by the solver during inference (i.e. $y_{n}$ for $n$ corresponding to the iteration index), respectively.

We first discuss IEs where the integral operator only involves a temporal integration (i.e. 1D), as discussed in Section 2. In analogy with the case of differential equations, this case can be considered as the one corresponding to ODEs.

These IEs are given by an equation of type

where $f$ is the free term, which does not depend on $\mathbb{y}$, while the unknown function $\mathbb{y}$ appears both in the LHS, and in the RHS under the sign of integral. The term $\int_{\alpha(t)}^{\beta(t)}G(\mathbb{y},t,s)ds$ is an integral operator $\mathcal{C}(D)\longrightarrow\mathcal{C}(D)$ from the space of integrable functions $\mathcal{C}(D)$ over some domain of $\mathbb{R}$, into itself. We observe that the variables $t$ and $s$ appearing in $G$ are both in $D$, and they are interpreted as time variables. We refer to them as global and local times, respectively, following the same convention of [ZFM⁺23]. The functions $\alpha$ and $\beta$ determine the extremes of integration for each (global) time $t$. Common choices for $\alpha$ and $\beta$ include $\alpha(t)=0$ and $\beta(t)=t$ (Volterra equations), or $\alpha(t)=a$ and $\beta(t)=b$ (Fredholm equations).

The fundamental question in the theory of IEs is whether solutions exists and are unique. It turns out that under relatively mild assumptions on the regularity of $G$ IEs admit unique solutions [Zem12]. Furthermore, the proofs in the first chapter of [Lak95] show the close relation between IEs and IDEs, as existence an uniqueness problems for IDEs are shown to be equivalent to analogous problems for IEs. Then the fixed point theorems of Schauder and Tychonoff are used to prove the results.

We now discuss the case of IEs where the integral operator involves integration over a multi-dimensional domain of $\mathbb{R}^{n}$. This is the IE version of PDEs, and they are commonly referred to as Partial Integral Equations (PIEs) when integration occurs on different components separately. We will consider equations where the multi-dimensional integral is obtained through multiple integrations. An equation of this type takes the form

where $\Omega\subset\mathbb{R}^{n}$ is a domain in $\mathbb{R}^{n}$, and $\mathbb{y}:\Omega\times\mathbb{R}\longrightarrow\mathbb{R}^{m}$. Here $m$ does not necessarily coincide with $n$.

PIEs and higher dimensional IEs have been studied in some restricted form since the 1800’s, as they have been employed to formulate the laws of electromagnetism before the unified version of Maxwell’s equations was published. In addition, early work on the Dirichelet’s problem found the IE approach proficuous, and it is well known that several problems in scattering theory (molecular, atomic and nuclear) are formulated in terms of (P)IEs. In fact, the Schrödinger equation can be recast as an IE [TF59]. See also [SB51], where bound-state problems are treated with the IE formalism.

The most striking difference between the procedure of solving an IE and an ODE, is that for an IE in order to evaluate at a single time point, one needs to know the solution for all time points. This is clearly an issue, since solving for one point requires that we already know a solution for all points. In order to better elucidate this point, we consider a simple comparison between the solution procedure of an ODE equation of type $\dot{\mathbf{y}}=f(\mathbf{y},t)$, and an IE of type $\mathbf{y}=f(t)+\int_{0}^{1}G(\mathbf{y},t,s)ds$.

Let us assume we are solving an ODE of type $\dot{\mathbf{y}}(t)=f(\mathbf{y},t)$ and that $\mathbf{y}$ is known at time points $t_{0},t_{1},\ldots,t_{k-1}$. Then one can obtain $\mathbf{y}$ at $t_{k}$ by means of the Euler method by using the known value at $t_{k-1}$ by taking small enough steps $\Delta t$ forward in time. In general, therefore, one starts by the initial condition $\mathbf{y}_{0}$ and determines the solution $\mathbf{y}$ at the points $t_{0},\ldots,t_{n}$ by taking small steps and representing the derivative as $\Delta f/\Delta t$ for small intervals $\Delta t$. Of course, more sophisticated methods are possible for the numerical solution of the ODE, but they essentially produce the next time point from the previous one in a sequential way. Let us now consider an analogous Fredholm IE to the ODE given above. This is a simple equation of type $\mathbf{y}=f(t)+\int_{0}^{1}G(\mathbf{y},t,s)ds$. Suppose we know $\mathbf{y}$ at time points $t_{0},\ldots,t_{k-1}$. In order to determine $\mathbf{y}(t_{k})$, we need to compute $f(t_{k})+\int_{0}^{1}G(\mathbf{y},t_{k},s)ds$, which requires us to know $\mathbf{y}$ over the full interval $[0,1]$, as $G$ is integrated over $[0,1]$. It is obvious that knowing a single time point for $\mathbf{y}$ (or a sequence of values) does not suffice anymore. In a Volterra type of equation the integral would be between $[0,t_{k}]$ (where the unknown value $t_{k}$ is included!), which does not really change the essence of the issue.

Although several methods can be employed to solve IEs, most (if not all) of them are based on the concept of iteration over some initial guess for the solution of the IE. Iterating on the initial guess produces a sequence of functions that then converges to a solution of the integral equation. More specifically, one can consider the von Neumann series of the integral operator as we now discuss. In fact, let us consider Equation 5, which can be rewritten as

where we assume for a moment that $T$ is a linear operator. Observe that if we can find the inverse of $\mathbb{1}-T$, then we obtain $\mathbf{y}$ as $(\mathbb{1}-T)^{-1}(f)$. This can be done by writing the von Neumann series for $(\mathbb{1}-T)^{-1}=\sum_{k=0}^{\infty}T^{k}$. This expression makes sense whenever the series $\sum_{k=0}^{\infty}T^{k}$ converges in operator norm, which is guaranteed in important cases such as when $\sum_{k=0}^{\infty}||T||^{k}$ converges (e.g. when $||T||<1$), while milder conditions on the convergence of the series exist as well. In such a situation when the von Neumann series is meaningful, we can then obtain $\mathbf{y}$ by iteratively applying $T^{k}$ to $f$. The nonlinear case is handled in a similar iterative procedure, which is called Method of Successive Approximations, or also Picard’s Iterations [Dav60]. It is in fact straightforward to convince oneself that under mild conditions, the method will output a solution of the integral equation. Conditions under which such succession is guaranteed to converge can be found in [Dav60]. A particularly well known case is when the integrand of the integral operator is contractive (i.e. Lipschitz with constant between $0$ and $1$) with respect to the variable $\mathbf{y}$. We give a proof of such approach for our setting, c.f. similar results found in [Dav60].

The following now follows easily.

In practice, the method of successive approximations is implemented as follows. The initial guess for the integral equation is simply given by the free function $f$ (i.e. $T^{0}(f)$), which is used to initialize the iterative procedure. Then, we apply $T$ to $\mathbf{y}^{0}:=T^{0}(f)$ to obtain a new solution $\mathbf{z}^{1}:=f(t)+T^{1}(\mathbf{y}^{0})$. We set $y^{1}:=r\mathbf{y}^{0}+(1-r)\mathbf{z}^{1}$ and apply $T^{2}$ to the solution $y^{1}$ and repeat. Here $r$ is a smoothing factor that determines the amount of contribution from the new approximation to consider at each step. As the iterations grow, the fractions of the contributions due to the smoothing factor $r$ tend to $1$. Observe that when we sum $r\mathbf{y}^{i}+(1-r)\mathbf{y}^{i+1}$ with $r=0$, we obtain the terms of the von Neumann series up to degree $i$ applied to $f$: $\sum_{0}^{i}T^{k}(f)$. The smoothing factor generally shows good empirical regularization for IE solvers, and we have set $r=1/2$ throughout our experiments, even though we have not seen any concrete difference between different values of $r$. This procedure is shown in Fugure 2.

In [DM88], Chapter 4, computations on the error bounds for the iterative procedure described above when the integrand function $G$ splits into the product of a kernel (see above) and a linear function $F$ are given. In the same chapter a detailed description of the Nyström approximation for the computation of the error is given as well. We describe a concrete realization of the iterative procedure discussed above in Section D, along with the learning steps for the training of our model. Moreover, we additionally observe that the procedure described above does not depend on $T$ being an integral operator or a general operator and, therefore, applying this methodology to the case where we have a transformer instead of $T$ is still meaningful, in the assumption that $T$ is such that the iterated series of approximations is convergent.

Depending on the specific IE that one is solving (e.g. Fredholm or Volterra, $1D$ or $(n+1)D$ etc.) the actual numerical procedure for finding a numerical solution can vary. For instance, a list of articles that showcase such a wide variety of specific methods for the solution of certain types of equations is [KR12, BRSS17, LR02, KGES19, PYD20]. Such variations upon the same theme of iterative procedure depend on finding the most efficient way of converging to a solution, finding the best error bounds, improving stability of the solver and substantially depend on the form of the integral operator. As our method is applied without the actual knowledge of the shape of the integral operator, but it actually aims at inferring (i.e. learning) the integral operator from data, we implement an iterative procedure which is fixed and depends only on a hyperparameter smoothing factor. This is described in detail in the next section. However, we point out that since the integrand, and therefore the integral operator itself, is learned during the training, one can assume that the model will optimize with respect to the procedure in a way that our iterations are in a sense “optimal” with respect to the target.

Thus far, our considerations on the implementation of IE solvers seem to point to a fundamental computational issue, since they entail a more sophisticated solving procedure than that of ODEs or PDEs. However, in various situations, even solving ODEs and PDEs through IE solvers presents significant advantages that are not necessarily obvious from the above discussions. The first advantage is that IE solvers are significantly more stable than ODE and PDE solvers, as shown for instance in the articles [KR12, Rok85, Rok90]. This in particular provides a new solution to the issue of underflowing during training of NODEs that does not consist of a regularization, but of a complete change of perspective. In addition, even though to solve an IE one needs to iterate, in general the number of iterations is not particularly high. In fact, in our experiments the total number of iterations turned out to be sufficient to be fixed between $4$ and $6$. However, when solving for instance an ODE, one needs to go sequentially through each time step. These can be in the order of the $100$ (as in some of our experiments). On the contrary, our IE solver processes the full time interval in parallel for each iteration. This results in a much faster algorithm compared to differential solvers, as shown in our experiments.

The solver procedure described in the previous subsection of course assumes that there exists a solution to start with. As mentioned at the beginning of the section, we treat equations of the second kind in this article also because the existence conditions are better behaved than for the equations of the first kind. We give now some theoretical considerations in this regard. We will also discuss when these solutions are uniquely determined. Existence and uniqueness are two fundamental parts of the well-posedness of IEs, the other being the continuity of solutions with respect to initial data.

A concise and relatively self-contained reference for the existence and uniqueness of solutions for (linear) Fredholm IEs is Section 4.5 in [Mor13]. In fact, it is shown in Theorem 4.43 that if a Fredholm equation has Hermitian kernel, then the IE has a unique solution whenever $\lambda$ is not an eigenvalue of the integral operator. For real coefficients, which is the case we are interested into, one can simply reduce the case to symmetric kernels, which are kernels for which $K(t,s)=K(s,t)$ for all $t,s$. Since in this article we have assumed $\lambda=1$, the condition becomes equivalent to saying that there is no function $\mathbf{z}$ such that $\int_{0}^{1}K(t,s)z(s)ds=z(t)$ for all $t$.

For more general (linear) integral operators (bounded on a Hilbert space), a similar result holds. In fact, from Theorem 4.44 in [Mor13] we have that a generalized Fredholm integral equation admits solutions if and only if the free function is orthogonal to each solution of the associated homogeneous adjoint equation. The latter admits the zero function as a solution (so the solution set is not empty), and is obtained from Equation 5 by deleting $f$, and by taking the adjoint of $T$ and the complex conjugate of $\mathbf{y}$. In the real case, the conjugate of $\mathbf{y}$ is $\mathbf{y}$ itself. Moreover, uniqueness is guaranteed if the associated homogeneous equation has only trivial solutions. In the case of nonlinear integral operators several existence and uniqueness conditions can be found in the literature. The reader is referred to [Dav60, Zem12, Lak95] for specific formulations. Generally speaking, such conditions are assumed on the integrand functions that determine the integral operator, in such a way that contractive theorems (such as Schauder’s and Tikhonoff’s) can be applied.

Observe that such formulations of the existence and uniqueness based on the contractive properties of the operator $T$ are particularly interesting in the case where the integral operator is replaced by a general neural network (between function spaces) which is not necessarily obtained through integration. In practice, when $T$ is a general neural network that is possibly nonlinear on all the entries, except with respect to the function $\mathbf{y}$, $T$ can be approximated by an integral equation using the following reasoning. It is known that Hilber-Schmidt operators on the Hilbert space of square integrable functions are approximated by integral operators – see e.g. [Mor13] Theorem 4.26. It is reasonable to assume that neural network operators are well behaved enough to be considered Hilbert-Schmidt operators. They therefore approximate some integral operator, and the training process therefore learns an integral equation.

More generally, for IEs of Urysohn or Hammerstein type, the existence and uniqueness problems are well known under much milder conditions, namely when the operator is completely continuous [Kra64, Kra]. In this situation, it is sufficient for the operator to have a nonzero topological index to guarantee that the corresponding IE admits a solution, and to study the problem of uniqueness one can determine the value of the topological index in a bounded subset of the Banach space in consideration, since this is directly related to the number of fixed points of the given IE.

The previous discussion, however, does not directly apply to the case when $T$ is a transformer, due to the fact that $T(\mathbf{y})$ is not linear in $\mathbf{y}$ in this case. Such equations can still be considered generalized Fredholm equations, and conditions on nonlinear operators $T$ being approximated by integral operators can be found in the literature, but the extent to which such equations are equivalent to integral equations is a fascinating question which will not be explicitly considered in this article.

As a word of warning, we mention that the general theory ensures the existence and uniqueness of solutions under some (mild) assumptions. Of course, in principle one should impose constraints to ensure that such assumptions are satisfied and that the results would apply. However, in our experiments we have observed good stability and good convergence without imposing any additional constraints. This does not apply in general, but we hypothesize that during optimization the model converges towards operators whose associated IE is well behaved, in order to avoid regimes of poor stability due to lack of solutions or lack of uniqueness of solutions. For different datasets such behavior might not be satisfied, and extra care in this regard might be needed.

NIE does not learn a dynamical system via the derivative of a function $\mathbb{y}$, as is the case for ODEs and IDEs. Therefore, we do not need to specify an initial condition in the solver during training and evaluation. In fact, the initial condition for IEs is encoded in the equation itself. For instance, taking $t=0$ in a Volterra or a Fredholm equation uniquely fixes $\mathbb{y}(\mathbb{x},0)$ for all $x$.

Therefore, we can specify a condition for IEs by constraining the free function $f(\mathbb{y},t)$. In what follows, we will make use of this paradigm several times. There are two immediate ways one could impose constraints on the free function. The simplest is to fix a value $\mathbb{y}_{0}$ and let $f(\mathbb{y},t)$ be fixed to be $\mathbb{y}_{0}$ for all $t$. Alternatively, one could choose an arbitrary function $f$ and keep this function fixed. In practice, the latter is conceptually more meaningful. For instance, in theoretical physics, when transforming the Schrödinger equation into an integral equation, on the RHS one can choose an arbitrary function $\psi(\mathbb{y},t)$ which corresponds to the wave function of free particles, i.e. without potential $V$. Applications of this procedure are found below in the experiments.

We consider two settings, where the integral operarator is modeled by a single hidden layer feed forward neural network of arbitrary width, or by an arbitrarily deep neural network.

We want to show that when we restrict ourselves to single hidden layer feed forward neural networks of arbitrary depth for our function $G_{\theta}$ in Equation 1, we can approximate a wide class of integral equations over a suitable subset of the space of functions. In the case of deep neural networks, we will argue that the NIE architecture can approximate any “regular enough” integral operator, where regularity will be described below. We restrict our considerations to the case of function spaces where the domain is $\mathbb{R}$, since the higher dimensional case is easily adapted from this discussion. We will therefore use $y$ instead of $\mathbf{y}$ to indicate the elements of the domain of the integral operators.

Let $T:C([0,1])\longrightarrow C([0,1])$ be an integral operator on the space of continuous functions, defined as $y\mapsto T(y)(t):=\int_{\alpha(t)}^{\beta(t)}G(y(s),t,s)ds$ for continuous functions $\alpha,\beta:[0,1]\longrightarrow[0,1]$, and a continuous $G:\mathbb{R}\times[0,1]\times[0,1]\longrightarrow\mathbb{R}$. In fact, for what follows we could consider $G$ as being Borel measurable, instead of imposing the more restrictive condition of being continuous. However, since in applications continuity is generally required, we impose this more restrictive conditions. Moreover, our discussion easily extends to the case when the definition intervals are $[a,b]$ instead of $[0,1]$ with simple modifications, and a similar approach also extends to higher dimensional integrals. We assume that $T$ is such that the corresponding integral equation of the second kind, i.e. Equation 1, admits a solution $y^{*}:[0,1]\longrightarrow\mathbb{R}$ in $C([0,1])$. Since $y^{*}$ is continuous, there exists a compact $K=[-k,k]$, for $k>0$, such that $y^{*}([0,1])\subset K$. Let us consider now a neighborhood $U_{K}$ of $y^{*}$ in the compact-open topology such that for all $y\in U_{K}$ we have the property $y([0,1])\subset K$. This could be, for instance, the space of functions $y$ mapping $[0,1]$ into the open $(-k,k)=K^{\circ}$. We can therefore restrict $G$ to the domain $K\times[0,1]\times[0,1]$, and we will still indicate this restriction by $G$ and the corresponding integral operator by $T$ (defined over the neighborhood $U_{K}$), for notational simplicity.

For an arbirtary chosen $\epsilon>0$, we want to show that we can approximate $T(y)$ with error at most $\epsilon$ in the metric induced by $C([0,1])$ on $U_{K}$ through an NIE integral operator $T_{\theta}(y)(t):=\int_{\alpha(t)}^{\beta(t)}G_{\theta}(y(s),t,s)ds$. To this purpose, let us set $Q=\sup_{[0,1]}|\beta(t)-\alpha(t)|$, and observe that by applying the Universal Approximation Theorem for single hidden layer feed forward neural networks (see [HSW89]) to the function $G:K\times[0,1]\times[0,1]\longrightarrow\mathbb{R}$, we can find a single hidden layer neural network $G_{\theta}:K\times[0,1]\times[0,1]\longrightarrow\mathbb{R}$ such that for all $t,s\in[0,1]$ we have $|G(y(s),t,s)-G_{\theta}(y(s),t,s)|<\epsilon/Q$. With such a $G_{\theta}$, for all functions $y\in U_{K}$ we have for any fixed $t^{*}$ in $[0,1]$

Therefore, uniformly on $t$ we have that

But this means that $d(T(y),T_{\theta}(y))<\epsilon$ with the metric $d$ on $U_{K}$ induced by that of $C([0,1])$.

We observe that while this approximation does not hold in complete generality, it is valid for a class of integral operators of importance, since we are usually interested in operators whose corresponding IE is admits continuous solutions, and we are interested in modeling the operator in a neighborhood of a solution. Moreover, under mild assumptions (e.g. discussed in Appendix B.4 the dependence of the solution on the initial data is continuous, and therefore the solutions to the equation for perturbed $f$ lie in a neighborhood of a solution $y^{*}$ obtained for $f$. So, our results apply in such important cases. Lastly, we point out that throughout the previous reasoning we have implicitly assumed that numerical integration is performed with infinite precision. Of course this is not the case in practice, but since we can reduce the numerical error in the integration procedure arbitrarily by choosing densely enough samples for a choice of the integration scheme, the error due to numerical integration can be rendered small enough so that the previous inequalities hold.

We now consider the case where we allow deep neural networks as in [LPW⁺17], Theorem 1. In this case, we argue that for any IE of the second kind as in Equation 1 where we set $T(y)(t):=\int_{\alpha(t)}^{\beta(t)}G(y(s),t,s)ds$ for a Lebesgue integral function $G$, we can approximate the integral operator $T$ with arbitrary precision. As a consequence, there is a NIE model which realizes any IE as in Equation 1 with arbitrary accuracy. We can proceed as for the case of single hidden layers neural networks above, with the main difference that applying Theorem 1 from [LPW⁺17] we do not need to restrict ourselves to a neighborhood $U_{K}$ of a solution $y^{*}$ of the integral equation, and the neural integral operator $\int_{\alpha(t)}^{\beta(t)}G_{\theta}(y(s),t,s)ds$ approximates $T$ uniformly with respect to $t$ for any $y\in C([0,1])$. Observe that in order to use Theorem 1 in [LPW⁺17] we need to precompose $G$ and $G_{\theta}$ by a characteristic function $\chi_{[0,1]}$, which does not affect the result.

We give some comments on the approximation properties of ANIE with respect to generalized Fredholm equations. For simplicity, we consider the case where the integration is performed only over time, even though the same reasoning can be extended to spatiotemporal domains. Let $T:C([0,1])\longrightarrow C([0,1])$ denote a Fredholm integral operator defined through the assignment $T(y)(t)=\int_{0}^{1}G(\mathbf{y}(t),t,\mathbf{y}(s),s)ds$. Observe that this integral form is more general than considered in Equation 1, and it follows the interpretation of integration in terms of self-attention (c.f. Appendix D.2, where the integration approximation used in this article is given in more detail).

Let us assume that the integral equation $y=f^{*}+T(y)$ admits a unique continuous solution $y^{*}\in C^{2}([0,1])$, and that $G$ is regular enough so that the equation admits unique solution in $C([0,1])$ for given functions $f$ in a neighborhood of $f^{*}$ in the compact open topology. Observe, that such well-posedness conditions are usually relatively mild (see for instance Appendix B.4 and references therein), and this is the main situation of interest in applications. Then, there exists a compact $K=[-k,k]$ such that $y^{*}([0,1])\subset K$ and we can choose a neighborhood $U_{K}$ of $y^{*}$ in the compact-open topology of $C([0,1])$ such that $y([0,1])\subset K$ for all $y\in U_{K}$. In fact, one can simply choose $U_{K}:=\{y\in C([0,1])\ |\ |y([0,1])|<k\}$. Under such hypothesis, there are numerical integration schemes that can approximate the integral $\int_{0}^{1}G(\mathbf{y}(t),t,\mathbf{y}(s),s)ds$ for any fixed choice of $t$ with arbitrary precision, upon choosing a number of points for evaluation that is sufficiently large. For instance, for any fixed $t$, the error for trapezoidal rules is bound by a term that goes to zero as $n$ grows, where $n$ is the number of points chosen in $[0,1]$ for approximating the integral, see Equation 2.1.6 in [DR07]. This term is the modulus of continuity

For each choice of $n$, there exists a compact $K_{n}:=[-k_{n},k_{n}]$ such that $y^{*}$ maps into $K_{n}$, and $G$ restricted to $K_{n}\times[0,1]\times K_{n}\times[0,1]$ has $\omega_{t}(1/n)<1/n$ for all $t\in K_{n}$. In this situation we can choose a neighborhood of $y^{*}$, $U_{K_{n}}$ such that $\omega_{t}(1/n)<1/n$ for all $t\in K_{n}$ for each choice of $y\in U_{K_{n}}$, and this numerical integration approximates the value of $T(y)(t)$ with arbitrarily high accuracy.

We indicate our numerical integration scheme using the formula

where $s_{i}$ indicates the $i^{\rm th}$ grid point of $\{t_{i}\}\subset[0,1]$. We can therefore obtain the evaluation of $T(y)$ at the grid points $t_{j}$ as