For some initial-boundary value problems (IBVP) in fluid mechanics, the solution of the partial differential equations (PDE) include discontinuous initial data or a discontinuity that forms in finite time, i.e. shockwaves. Numerical methods for solving these PDE must be specially tailored to properly resolve these discontinuities (LeVeque et al., 2002).

These shock-capturing methods are designed with the goal of sharply resolving a shock without inducing spurious oscillations, while also giving accurate solutions in smooth regions of the flow. One major breakthrough in this effort was the development of high-resolution methods (Harten, 1983) , as these methods were capable of achieving secondorder accuracy without introducing spurious oscillations around shocks. These methods gave rise to a class of high-resolution methods called essentially non-oscillatory (ENO) schemes (Harten et al., 1987) that measure the smoothness of the solution on several stencils, and then compute the flux based on the smoothest stencil to avoid interpolating through the discontinuity. These schemes are nonlinear (even when the PDE are linear) since the interpolation co-efficients depend on the solution. These ideas were then modified to create WENO-JS (weighted ENO-Jiang Shu) methods (Jiang & Shu, 1996), which again compute the smoothness on several stencils. However, instead of taking only the smoothest stencil, these methods take a weighted average of the fluxes predicted on each stencil to emphasize the smoother ones. When each stencil is equally smooth, the weights are designed to cause the method to converge to the constant coefficient scheme that maximizes the order of accuracy over the union of the sub-stencils, which gives these methods a high order of accuracy for smooth solutions.

Many efforts have built on the original WENO-JS schemes by modifying the smoothness indicators (Ha et al., 2013; Kim et al., 2016; Rathan & Raju, 2018b), modifying the nonlinear weights (Borges et al., 2008; Castro et al., 2011; Rathan & Raju, 2018a), and using WENO-JS as part of a hybrid scheme (Li & Qiu, 2010; Peer et al., 2009). While some of these references build off of each other rather than starting from WENO-JS, we will base our method on WENO-JS because our strategy for developing the method does not resemble other methods. However, our methodology could easily adopt various improvements that have been made to WENO-JS.

One commonality that has persisted since the original ENO scheme is a reliance on human intuition in shock-capturing method design, particularly in the nonlinear aspects of the schemes, i.e. smoothness indicators and weighting functions. While they have been well studied, there is no reason to believe that they are optimal. Efforts have been made to develop optimal spatial discretization methods by minimizing wave propagation errors (Kim & Lee, 1996; Lele, 1992; Liu, 2013; Tam & Webb, 1993) and minimizing error over certain frequency ranges (Zhang & Yao, 2013), and some of these techniques have even been combined with shock-capturing schemes (Fang et al., 2013; Wang & Chen, 2001). However, designing the optimization problem still requires human intuition with regards to balancing competing goals, rather than attempting to learn an optimal scheme from data in an unbiased way.

Over the past decades, machine learning has become ubiquitous in data analysis and is increasingly seen as having potential to improve (or reformulate) numerical methods for PDEs. Lagaris et al. (Lagaris et al., 1998) parameterized the solution to a PDE as a neural network and optimized the weights to minimize the residual of the solution. Yu et al. (Yu et al., 2018) trained a neural network to classify the local smoothness and apply artificial viscosity based on this classification. Bar-Sinai et al. (Bar-Sinai et al., 2019) used simulation data to embed coarse graining models into fi-nite difference schemes involving neural networks, allowing them to achieve low error on relatively coarse grids. Pfau et al. (Pfau et al., 2018) parameterized the eigenfunctions of eigenvalue problems as a neural network and cast the training as a bilevel optimization problem to reduce bias, resulting in significantly decreased memory requirements. Hsieh et al. (Hsieh et al., 2019) attempted to learn domain specific fast PDE solvers by learning how to iteratively improve the solution using a deep neural network, resulting in a 2-3 times speedup compared to state of the art solvers.

In the current work, we attempt to train a neural network to improve WENO5-JS. Our goal is to get closer to the optimal nonlinear finite-volume coefficients while introducing a minimal amount of bias. Unlike other references, we do not directly change the smoothness indicators or nonlinear weights of the method. Instead, we use a neural network to perturb the finite-volume coefficients determined using the original smoothness indicators and nonlinear weights of WENO5-JS. We attempt to learn an optimal function for this perturbation using data generated from waveforms that are representative of solutions of PDEs. These modifica-tions result in a finite-volume scheme that diffuses fine-scale flow features and discontinuities less severely than WENO5-JS. We start in the next section by giving a more detailed description of the proposed algorithm.

2.1. Description of WENO-NN

Although we focus on WENO5-JS in this paper, our approach could generally be used to enhance any shock capturing method (or perhaps any numerical method). The proposed algorithm involves pre-processing the flow variables on a stencil using a conventional shock capturing method and feeding those results into a neural network. The neural network then perturbs the results of the shock capturing method. Post-processing is then applied to the output of the neural network to guarantee consistency (Bar-Sinai et al., 2019) (or, more generally, could be used to enforce

other desirable properties). Hence, the augmented numerical scheme takes on many properties of the original. For example, applying the method to WENO5-JS results in an upwind-biased finite volume method with coefficients that depend on the local solution. The steps of the algorithm for enhancing WENO5-JS can be seen in algorithm ??.

We use WENO5-JS to pre-process the input data, so that the input to the neural network is the set of finite-volume coefficients found by WENO5-JS. We found that including this pre-processing step significantly improved performance. Once the nonlinear weights are determined according to the WENO5-JS algorithm, the coefficients for each cell average are computed as

These five coefficients are the inputs to the neural network, which outputs a change in each coefficient, . Our neural network uses 3 hidden layers, each with 3 neurons. We deliberately make the network as small as possible to reduce the computational cost of evaluating it. We are able to use such a small network because assuming that the WENO5-JS coefficients are a useful model input is a strong prior, so WENO5-JS performs a significant amount of the required processing. regularization is applied to the output of the neural network to penalize deviations from WENO5-JS, which encourages the network to only change the answer supplied by WENO5-JS when an improved result is expected. The new coefficients are computed by subtracting the change in coefficients from the old coefficients.

Additionally, the size of the input space is reduced by scaling cell averages within the stencil as

If all the cell averages have the same value, the scaling equation fails so the value at the cell edge is simply set to the cell average value.

To guarantee that WENO-NN is consistent, we apply an affine transformation to these coefficients that guarantees that they sum to one (Bar-Sinai et al., 2019). We derive this transformation by solving the optimization problem

which can be reformulated with the substitution to pose the problem as finding the minimum norm solution to an under-constrained linear system

which has the analytical solution

One can use the same approach to enforce arbitrarily high orders of accuracy since the optimization problem has an analytical solution for any constraint matrix of sufficiently high rank

This optimization problem has analytical solution is invertible.

We verify that our constraint is satisfied by looking at the convergence rate of WENO-NN for a smooth solution. For this test case, we will simply use WENO-NN and WENO5-JS to take the derivative of , and compare the results to the analytical solution using the error metric

In Figure 1, we can see that WENO-NN achieves first order accuracy, which confirms that the constraint is satisfied.

Figure 1. Convergence of WENO-NN and WENO5-JS for smooth solutions

We also see that, as expected, WENO5-JS converges at fifth order as . However, when discontinuities are present it is not possible to achieve better than first order accuracy with any finite volume method (LeVeque et al., 2002). Despite this fact, it is advantageous to use WENO5-JS over WENO3-JS in such situations, as WENO5-JS still tends to give lower error in discontinuous problems (Shu, 1998), which is why we chose to use WENO5-JS for processing the cell average values despite the fact that WENO-NN ends up being first-order accurate. Similarly, we see that for some discontinuous problems, WENO-NN gives lower error than WENO5-JS. If a higher order of accuracy is desired in smooth regions of the flow, one could develop a hybrid method using WENO-NN and any high-order method.

2.2. Other Numerical Methods Used

For all simulations shown, we use a third-order TVD RungeKutta scheme (Gottlieb & Shu, 1998) as our time-stepping method

For flux-splitting, we use a Lax-Friedrichs flux splitting procedure (Shu, 2003)

Figure 2. Comparing error trends between (A) exact generated data and (B) simulation results

In this expression, f(u) is defined as the flux of a 1-D hyperbolic conservation law . When solving the 1-D Euler equations, we apply the flux splitting to the characteristic decomposition of the system. For our numerical Riemann solver, we use the Lax-Friedrichs method (Chu, 1979).

We construct our training data directly from known functions that we expect to represent the waveforms that WENONN will encounter in practice. For each datapoint, we start with some function u(x) and a discretized domain of n cells. The cell average is evaluated on each cell as

and because we chose the form of u(x) we can evaluate the cell average analytically. We also evaluate the function value on the cell boundary as analytically. We then move along the domain and form the dataset based on the stencil size. So for WENO-NN, one datapoint involves 5 cell averages as the input with the function value on the cell boundary as the output. The functions we use when creating the dataset are step functions, sawtooth waves, hyperbolictangent functions, sinusoids, polynomials, and sums of the above.

When adding a new entry to the dataset, we first check to see if it is close to other points already present in the sense. Sufficiently close points are not added to the dataset to prevent redundant data that will slow down the training process. The resulting dataset has 75241 entries.

When training the network, we use the Adam optimizer (Kingma & Ba, 2014), split the data into batches of 80 points to estimate the gradient, and optimize for 10 epochs using the Keras package in python (Chollet et al., 2015). We trained the network from many different randomly chosen initial guesses of the parameters, and chose the best one based on performance in simulating the linear advection of a step function. We apply regularization with a constant of to the neural network output, and find that when splitting the data into a training set of 80% of the data and a validation set of the other 20% of the data our in-sample error is 0.569 and the out-of-sample error is 0.571, averaged from 100 trials of training on the dataset, so overfitting within the generated dataset is not a concern. This difference is so small because the model we are training is of relatively low complexity, and is essentially underfitting the generated dataset. We use mean squared loss as our objective function to minimize.

Despite the fact that we do not see overfitting within the generated dataset, we still observe overfitting when we apply the method to an actual simulation. Figure 2 shows the average training error, average validation error, and average error when using the method to simulate a PDE for different regularization values of the neural network output. The training and validation error are computed using the mean square error,

while the simulation error is computed by using the learned numerical method to linearly advect a step function and computing the error at the end of the simulation,

One can see that adding regularization causes error to increase in both the training and validation datasets but decreases the error in the simulation results. Hence, we can see that we are overfitting to the training data, but because the validation data does not show this, we can conclude that the dataset does not exactly match the distribution we are trying to approximate.

4.1. Advection Equation

These results will focus on comparing WENO5-JS to WENO-NN. Every WENO-NN result we show in this paper was generated using the same neural network with the same weights. As such, our numerical method is broadly applicable to problems not discussed in this paper, in contrast with many machine learning solutions that are problem-specific. The first test case we look at is the linear advection of a step function on a periodic domain. Mathematically, this IBVP is posed as

For this simulation, we set c = 1 and L = 2. We split the domain into 100 cells, use a CFL number of 2/3, and run the simulation for 50 periods for a total time of T = 100. Figure 3 shows the solution of this PDE for WENO5-JS and WENO-NN at t = 0, 20, 50 and 100. The solution at t = 0 is also the exact solution at all the other times plotted.

One can see that the solution using WENO-NN provides a closer visual fit to the exact solution, as WENO5-JS diffuses the discontinuity more significantly than WENO-NN. WENO5-JS also introduces noticeable overshoot behind the discontinuity. The neural network has the interesting property that the waveform is nearly invariant to its propagation, while WENO5-JS continues to diffuse the solution. This behavior can be explained by examining the artificial fluid properties associated with the modified equation obtained by Taylor series expansion (assuming linearity of the scheme). The modified PDE is

The expansions give expressions for the artificial viscosity,

dispersion, and hyperviscosity, 0 after making the substitutions

, and are computed as

Figure 4 shows these quantities for WENO5-JS. In order to estimate the contribution of each term, we approximated the higher-order spatial derivatives using standard finite-volume methods, and scale each by the magnitude of that derivative. For example, the influence of artificial viscosity is computed as

Hence, we ignore regions of the flow where the coefficient may signify that artificial viscosity (or other properties) are being added when they would have a negligible effect because the derivative is small.

One can see that for WENO-JS there is no viscosity or dispersion, as the method is designed such that on each substencil and 0, so WENO5-JS applies only hyperviscosity. The method applies a small amount of negative hyperviscosity near the discontinuity. As time goes on and the discontinuity continues to diffuse, the influence of hyperviscosity decreases.

Figure 5 shows that unlike WENO5-JS, WENO-NN adds both artificial viscosity and dispersion to the solution. We see that near the discontinuity, negative viscosity is being added, which apparently provides the anti-diffusion that causes the discontinuity to retain its steepness, while hyperviscosity is applied to stabilize the solution.

We obtain a quantitative picture of the error in figure 6. We plot the error over time (measured to the exact solution), as well as the total variation, , to indicate when oscillations have been induced in the solution. We also measured the width over which the discontinuity is spread by counting the cells that have an error above a certain threshold (in this case chosen to be 0.01) and multiplying this number by since there are two discontinuities in the simulation.

Figure 3. Numerical solutions of the advection equation at t = 0, 20, 50 and 100 using (A) WENO-NN and (B) WENO5-JS. Note that the curves in (A) for t > 0 are indistinguishable.

Figure 4. Influence of (A) artificial viscosity, (B) dispersion, and (C) hyperviscosity of WENO5-JS

Figure 5. Influence of (A) artificial viscosity, (B) dispersion, and (C) hyperviscosity of WENO-NN

Figure 6. Comparing (A) error, (B) total variation, and (C) discontinuity width over time for WENO-NN and WENO-JS

Figure 7. error at the end of the simulation for (A) WENO-NN and (B) WENO5-JS

The figure shows that WENO-NN decreases the error by almost a factor of 2 1. Although the total variation spikes at the start of the WENO-NN simulation, it is damped out and returns back to approximately the true value of 2, while the WENO5-JS total variation steadily climbs to above 2.04. We see a similar behavior in the discontinuity width, where WENO-NN reaches its steady value relatively quickly, while WENO5-JS continues to spread.

In order to determine how WENO-NN performs in different settings, the spatial and temporal discretizations were varied, and the error at the end of the simulation was measured. We again use a domain of length 2 and simulate for 50 periods. These results can be seen in figure 7.

We can see that WENO-NN tends to outperform WENO5-JS in regions where the spatial discretization is fine, but results in a larger error for coarse discretizations. To further compare the methods, figure 8 shows the error against the runtime for the two methods within a range of CFL values. We will only look at moderate CFL numbers, between 0.25 and 0.75, as stability becomes a concern for both methods above this range, and it becomes inefficient to run the simulation with CFL numbers below this range. We will also restrict the cell width to be below 0.025, as coarser meshes cause the final waveform to be unrecognizable compared to the exact solution for both methods, so the error comparison becomes meaningless. We see that when the CFL number is of a moderate value and the grid is sufficiently refined, WENO-NN typically achieves lower errors with smaller run time than WENO5-JS.

Figure 8. Comparing the error and runtime of WENO-NN and WENO5-JS for

4.2. Inviscid Burgers’ Equation

We next consider the inviscid Burgers’ equation. Unlike the linear advection equation that included only contact (initial) discontinuities, the inviscid Burgers’ equation results in shocks from smooth initial data. The distinction here is important: for a shock, the dynamics of the PDE will drive the solution towards a discontinuity, countering any diffusive effects associated with the numerics. We will again consider periodic boundary conditions, though we will start the simulation with a Gaussian as the initial condition. Hence, the IBVP is posed as

We simulate the problem for a time of T = 4 on a domain of length L = 2, and a value of k = 20. We first approximate the exact solution by solving this simulation with and for a total of 6400 cells and 25601 timesteps. What we see is that the is roughly the same for WENO5-JS and WENO-NN, as shown by Figure 9. Hence, we should expect the method to perform similarly to WENO5 in the presence of a shock.

Figure 9. Comparing error vs. grid spacing of WENO-NN and WENO5-JS for the inviscid Burgers’ equation

4.3. 1-D Euler Equations

The last test case we will look at is the Shu-Osher problem, a test case involving the 1-D Euler equations. Note that the method was also tested on the Sod problem, but because this test case did not lead to any conclusions not drawn from either the advection equation or the inviscid Burgers’ equation, these results have been omitted. The Shu-Osher problem is a model problem for turbulence-shockwave interactions. It involves the following equations and initial conditions

The simulation takes place on a domain of length L = 10 and is run until a final time of is set to 0.2. We first obtain an approximately exact solution by discretizing the solution into 12800 cells and 10240 time-steps and use WENO5-JS for the simulation. This grid is fine enough to consider the solution exact. We then solve the problem

Figure 10. Comparing (A) density, (B) pressure, and (C) velocity of WENO-NN and WENO5-JS to the exact solution for the Shu-Osher problem

Figure 11. Zoomed in view of the turbulent section for different grid resolutions of (A) 250, (B) 300, (C) 800, and (D) 3200 cells for the Shu-Osher problem

using 300 cells and 240 time-steps using both WENO5-JS and WENO-NN, and compare the numerical results to the exact solution. Figure 10 shows the density, pressure, and velocity at the end of the simulation.

The most interesting aspect of the solution is the highly oscillatory section of the density profile, which is considered to behave similarly to turbulence. Figure 11 shows a zoomed in view of this section at different grid resolutions.

One can see that the neural network diffuses the oscillations significantly less than WENO5 for coarse grids, which is an encouraging result in terms of simulating actual turbulence. As the mesh is further refined, the WENO-NN appears to overcompensate, though on the very fine grid both WENO-J5 and WENO-NN are similar (provided WENO-NN is stable, then it is constrained to converge as at least first-order).

By training a neural network to process the outputs of the WENO5-JS algorithm, we were able to improve its accuracy, particularly in problems where the artificial diffusion introduced in WENO5-JS was excessive. While WENONN is more expensive per evaluation than WENO5-JS, it achieved lower errors on coarser grids, which indicates some potential to be useful more generally. We trace these performance improvements to increased flexibility in the neural network compared to WENO5-JS, as it can add artificial viscosity and dispersion while WENO5-JS coefficients are constrained to make these quantities zero. By analyzing the advection of a step function, we found that WENO-NN applies negative artificial viscosity near the discontinuity, which allows it to maintain its sharp profile (this takes place sometime into the simulation after the initial profile has been slightly smoothed due to artificial viscosity that prevents spurious oscillations). We then observe similar behavior in the Shu-Osher problem, where we see that WENO5-JS diffuses the fine features of the solution more than WENO-NN. However, we also found that at certain resolutions WENONN applies too much negative artificial viscosity, resulting in too much amplification of these fine scale features, though this amplification does not develop into an instability. For true shocks, as opposed to contact discontinuities, we found that our method performs very similarly to WENO5-JS.

There is still room for improvement for the method. We believe that the development of training data could be improved. Perhaps better results could be achieved by formulating the training data in a more systematic way. Some papers generate the data directly from simulations (Bar-Sinai et al., 2019). While this approach results in neural networks that are specific to the equation that they were trained on, it could lead to methods that outperform methods designed to work well for general problems. Another drawback of WENO-NN is that it does not inherit the high-order convergence of WENO5-JS. It would be an improvement to the method to be able to structure the network such that its coefficients more quickly converge to those of either WENO5-JS or the constant coefficient scheme that maximizes order of accuracy in the presence of smooth solutions. However, this must be done in a way that does not interfere with predictions in non-smooth regimes that benefit from low-order behavior, which is a non-trivial task. Until such a method is developed, one would need to use WENO-NN as part of a hybrid scheme if higher order convergence is desired in smooth regions.

Another outstanding issue with machine-learned schemes is stability. The WENO-NN scheme used here seemed to inherit the stability of the underlying WENO5-JS scheme that it was based on, but this needn’t have been the case, and we cannot offer proof or an estimate for the maximal CFL.

In future work, we aim to test the method on large-scale, multidimensional problems, since we would expect the ben-efits seen in 1-D problems to be more significant when multiple spatial dimensions are present.

The authors declare that they have no conflict of interest.

This material is based upon work supported by the National Science Foundation Graduate Research Fellowship under Grant No. 1745301

Bar-Sinai, Y., Hoyer, S., Hickey, J., and Brenner, M. Learn- ing data-driven discretizations for partial differential equations. Proceedings of the National Academy of Sciences, 116(31):15344–15349, 2019.

Borges, R., Carmona, M., Costa, B., and Don, W. An improved weighted essentially non-oscillatory scheme for hyperbolic conservation laws. Journal of Computational Physics, 227(6):3191–3211, 2008.

Castro, M., Costa, B., and Don, W. High order weighted es- sentially non-oscillatory weno-z schemes for hyperbolic conservation laws. Journal of Computational Physics, 230(5):1766–1792, 2011.

Chollet, F. et al. Keras. https://keras.io, 2015.

Chu, C. Numerical methods in fluid dynamics. In Advances in Applied Mechanics, volume 18, pp. 285–331. Elsevier, 1979.

Fang, J., Li, Z., and Lu, L. An optimized low-dissipation monotonicity-preserving scheme for numerical simulations of high-speed turbulent flows. Journal of Scientific Computing, 56(1):67–95, 2013.

Gottlieb, S. and Shu, C. Total variation diminishing runge- kutta schemes. Mathematics of computation of the American Mathematical Society, 67(221):73–85, 1998.

Ha, Y., Kim, C., Lee, Y., and Yoon, J. An improved weighted essentially non-oscillatory scheme with a new smoothness indicator. Journal of Computational Physics, 232 (1):68–86, 2013.

Harten, A. High resolution schemes for hyperbolic con- servation laws. Journal of computational physics, 49(3): 357–393, 1983.

Harten, A., Engquist, B., Osher, S., and Chakravarthy, S. Uniformly high order accurate essentially non-oscillatory schemes, iii. In Upwind and high-resolution schemes, pp. 218–290. Springer, 1987.

Hsieh, J., Zhao, S., Eismann, S., Mirabella, L., and Ermon, S. Learning neural pde solvers with convergence guarantees. arXiv preprint arXiv:1906.01200, 2019.

Jiang, G. and Shu, C. Efficient implementation of weighted eno schemes. Journal of computational physics, 126(1): 202–228, 1996.

Kim, C., Ha, Y., and Yoon, J. Modified non-linear weights for fifth-order weighted essentially non-oscillatory schemes. Journal of Scientific Computing, 67(1):299– 323, 2016.

Kim, J. and Lee, D. Optimized compact finite difference schemes with maximum resolution. AIAA journal, 34(5): 887–893, 1996.

Kingma, D. and Ba, J. Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980, 2014.

Lagaris, I., Likas, A., and Fotiadis, D. Artificial neural net- works for solving ordinary and partial differential equations. IEEE transactions on neural networks, 9(5):987– 1000, 1998.

Lele, S. Compact finite difference schemes with spectral- like resolution. Journal of computational physics, 103(1): 16–42, 1992.

LeVeque, R. J. et al. Finite volume methods for hyperbolic problems, volume 31. Cambridge university press, 2002.

Li, G. and Qiu, J. Hybrid weighted essentially non-oscillatory schemes with different indicators. Journal of Computational Physics, 229(21):8105–8129, 2010.

Liu, Y. Globally optimal finite-difference schemes based on least squares. Geophysics, 78(4):T113–T132, 2013.

Peer, A., Dauhoo, M., and Bhuruth, M. A method for improving the performance of the weno5 scheme near discontinuities. Applied Mathematics Letters, 22(11): 1730–1733, 2009.

Pfau, D., Petersen, S., Agarwal, A., Barrett, D., and Stachen- feld, K. Spectral inference networks: Unifying deep and spectral learning. arXiv preprint arXiv:1806.02215, 2018.

Rathan, S. and Raju, G. Improved weighted eno scheme based on parameters involved in nonlinear weights. Applied Mathematics and Computation, 331:120–129, 2018a.

Rathan, S. and Raju, G. A modified fifth-order weno scheme for hyperbolic conservation laws. Computers & Mathematics with Applications, 75(5):1531–1549, 2018b.

Shu, C. Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws. In Advanced numerical approximation of nonlinear hyperbolic equations, pp. 325–432. Springer, 1998.

Shu, C. High-order finite difference and finite volume weno schemes and discontinuous galerkin methods for cfd. International Journal of Computational Fluid Dynamics, 17(2):107–118, 2003.

Tam, C. and Webb, J. Dispersion-relation-preserving finite difference schemes for computational acoustics. Journal of computational physics, 107(2):262–281, 1993.

Wang, Z. and Chen, R. Optimized weighted essentially nonoscillatory schemes for linear waves with discontinuity. Journal of Computational Physics, 174(1):381–404, 2001.

Yu, J., Hesthaven, J., and Yan, C. A data-driven shock capturing approach for discontinuous galerkin methods. Technical report, 2018.

Zhang, J. and Yao, Z. Optimized explicit finite-difference schemes for spatial derivatives using maximum norm. Journal of Computational Physics, 250:511–526, 2013.