# Exploratory Study of a Green Function Based Solver for Nonlinear Partial Differential Equations
## 3 Numerical Implementation
### Spatial Discretization
\begin{split}\{y_{j}\}:\big{|}G_{2q,i+1}\big{|}(x_{i+1},\vec{y})& =\big{|}G_{2q,i}\big{|}(x_{i},\vec{y}),\quad\quad\{x_{i}\}:\begin{cases} \dfrac{d}{dx}G_{2q,i}(x,\vec{y})=0\\ \text{subject to }x_{i}\in\Omega_{i}=[y_{i},y_{i+1})\end{cases},\\ G_{2q,i}(x,\vec{y})&=\prod_{j=1}^{2q}\Big{[}x-\zeta_{s(i)+j} \Big{]},\;\;\{\zeta_{j}\}=\underbrace{\{-y_{q},\ldots,-y_{1},y_{1},\ldots,y_{ N}\}}_{q\text{ times}},\\ s(i)&=\{0,\ldots,N-q,\underbrace{N-q,\ldots N-q}_{q\text{ times}}\},\end{split}  \tag{36}\]  
where \(s(i)\) serves as a re-indexing for the half-domain problem, adapted from what was introduced by (1976). The equioscillation property has to be imposed over the canonical function (1976) of the interpolation as it serves as a measure of its error,  
\[ \begin{split} R(x,q,\Omega_{i})&=u(x,t)-u_{pw}(x,t) =\frac{1}{(2q+2)!}G(x)u^{(2q+2)}(\zeta,t),\;\;x\in\Omega_{i}=[y_{i},y_{i+1}) \\ G(x)&=\prod_{j=-N}^{N}\big{(}x-x_{j}).\end{split}  \tag{37}\]Recall that, as this expression is a consequence of the application of [PERSON]'s theorem, the derivative must be evaluated in some point \(\zeta(x)\) in the interval we are doing the reconstruction, hence the only parameters we have to improve the error are the nodes defining the canonical function, \(\{x_{i}\}\). Therefore, optimizing the function towards the reduction of its maximal absolute value will produce an optimized grid in what is called the _minimax_ problem [59]. Regarding the specific solver we use for the iterative process, we apply the standard Levenberg-Marquardt algorithm as described in [60] that will have the error vector, \(\vec{\epsilon}\), with \(N\) components as the cost function to be reduced to zero. \[\epsilon_{i}=\left\|\frac{\big{|}G_{2q,i+1}\big{|}(x_{i+1},\vec{y})}{\big{|}G_ {2q,i}\big{|}(x_{i},\vec{y})}-1\right\|. \tag{38}\]  
We present in Figure 3 an example on how the iterative procedure produces an optimized grid and the shape of the canonical functions for both interpolations, the one defined over the \(2q\) points where equioscillation must be achieved and the optimal one over which the unknown function will be discretized and interpolated. In this figure we can see how the optimal grid \(\{x_{i}\}\) already presents a two-zones distribution, with an almost homogeneous grid distribution close to the center of the domain. We can actually force this kind of two-zones distribution by fine tuning the Hermanns Navarro algorithm and pre-imposing that this region will exist. As a result, we define a set of boundary points, \(b\), that will be the ones optimized after the \(N-b\) homogeneous ones, redefining the error vector, now of dimension \(b\) for the iterative process as,  
\[\epsilon_{b,i}=\left\|\frac{\big{|}G_{2q,i}\big{|}(x_{i},\vec{y})}{\big{|}G_ {2q,N_{b}-1}\big{|}\big{(}x_{N_{b}-1},\vec{y})}-1\right\|,\quad i=N_{b},...,N \tag{39}\]  
Hence, instead of imposing an equioscillation over the whole extrema of the \(2q\) polynomials, we will impose the value of the already equioscillating homogeneous inner \(N-b\) point distribution and obtain the points that adapt to it. The algorithmic logic for the forced two-zones problem can be described in the following manner:  
* Select the number of points \(b\) to be optimized and obtain the limit of the homogeneous zone, \(N_{b}=N-b\), defining the grid step at this region as \(\Delta x_{h}=N_{b}/N^{2}\). * Obtain the value of \(\big{\|}G_{2q,i}(x_{i},\vec{y})\big{\|}\) for an \(x_{i}\) that does not require a stencil, \(\vec{y}\), composed of points from the region to optimize. See the initial condition of Figure 3 as a visual key. * Compute the error vector \(\epsilon_{b,i}\) and perform the optimization process over the \(2q\) interpolation with the \(b\) points as the only unknowns. Therefore, we can define two separated regions, an inner, homogeneous, and a boundary one, for the whole \(2N+1\) grid points.