# Exploratory Study of a Green Function Based Solver for Nonlinear Partial Differential Equations
## 3 Numerical Implementation
### Spatial Discretization  
The chosen spatial discretization and its corresponding interpolation of the unknown function will drive the main error of the numerical method generating a the time evolving accumulated error, the evolution of which will be dependent of the timestep, \(\tau\). The original works of [PERSON] and [PERSON] (2001); [PERSON] and [PERSON] (2001); [PERSON] and [PERSON] (2001) used a simple zero-th-order discretization were both the Green function and the unknown function were evaluated over each of the points when computing the time marching integral scheme. Furthermore, while it is possible to keep this methodology, or even improve it by means of a Chebyshev-like grid distribution, in order to define a general numerical methodology that can be extended to a high performance numerical solver with an important number of grid points, \(N\), we will explore the application of the method by means of piecewise interpolation, trying to recover the usual methodology for numerical computation such as _stencils_ definitions and parallelization capabilities. Moreover, we choose a piecewise interpolation that, by optimizing the position of the grid points reduces the possible interpolation errors at the boundaries of the domain. Therefore, for the whole set of \(2N+1\) points we are using, the interpolation will be of order \(2q+1\), being able to increase or reduce this value in case additional precision or frequency representation are needed. The expression of the unknown function approximated by piecewise interpolants, \(u_{pw}\), will be defined as,  
\[ \begin{split} u_{pw}(x,t)&=\sum_{i=-M+1}^{M-1}H(x, \Omega_{i})I_{i}(x,t)+H(x,\Omega_{-M})I_{-M}(x,t)+H(x,\Omega_{M})I_{M}(x,t),\\ I_{i}(x,t)&=\sum_{j=-q}^{q}u\left(x_{i+j,t}\right) l_{p,j}(x-x_{i}),\quad M=N-q,\\ \Omega&=[a,b)\longrightarrow H(x,\Omega)=\text{Heaviside }(x-a)-\text{Heaviside}(x-b),\end{split}  \tag{34}\]  
where \(H\) stands for the rectangular function that defines the domain of validity for each of the piecewise interpolations, \(\Omega_{i}\), and \(l_{p,j}\) is the Lagrange polynomial distributed over the \(2q+1\) points of the stencil and centered at each of the points \(x_{i}\). Each of the domains of validity will be delimited by a set of \(N\), intermediate points, \(\{y_{j}\}\) in the following manner,  
\[\Omega_{-N}=[-1,y_{-N+1}),\quad\Omega_{i}=[y_{i},y_{i+1}),\quad\Omega_{N}=[y_{ N-1},+1]. \tag{35}\]  
where the last domains of validity, i.e., \(\Omega_{\pm M}\) have been sliced too in \(q\) subdomains to apply systematically the following methodology, yet these new sliced subdomains will not be domains of validity by definition, as the last stencil is valid all over \(\Omega_{M}\). We have normalized the ending points of the numerical domain at \(\pm 1\), yet the general \([x_{-N},x_{N}]\) can be easily recovered by a renormalization of the grid. The reason for this particular interpolation method is rather simple: it was shown by [PERSON] (1976) that the combination of collocation techniques and piecewise interpolation produce, for relative small ratios \(q/N\), a set of grid points that distributes homogeneously at the center of the domain and then adapts at the boundary to avoid the interpolation errors in the likes of the Chebyshev-Gauss-Lobatto grid distribution. To achieve such a grid distribution an iterative process is required as the mathematical problem has no analytic solution. The condition that is imposed is the usual equioscillation property (1976) over a lower order, \(2q\), piecewise interpolation defined using the \(\{y_{j}\}\) points to then recover the actual grid point distribution, \(\{x_{i}\}\) from their extrema. We attempt to summarize the whole process in a mathematical description, only for half of the domain as the whole \(2N+1\) grid distribution is symmetrical. \[