Exact Nonlinear Spectral Analysis of Nonsmooth WENO-JS Solutions

. The ﬁ fth-order accurate Weighted Essentially Nonoscillatory space discretization developed by Jiang and Shu (WENO-JS) is studied theoretically. An exact Nonlinear Spectral Method (NSM) is developed based on an innovative yet simple methodology. The NSM explains the behaviour of nonsmooth solutions because it is valid for arbitrary modi ﬁ ed wave numbers (MWN). The NSM clari ﬁ es the e ﬀ ects of the time integration methods and the Courant number. The mode isolation assumption, extensively used to analyse WENO-JS, is elucidated, and several novel ﬁ ndings are presented. The improved performance of the combination of WENO-JS with the forward Euler time integration method, compared to the Linear Fifth-Order Upwind discretization, is illustrated. The overdamping of the combination of WENO-JS with the popular third-order total variation diminishing Runge-Kutta method is discovered. Thus, the NSM covers several gaps in the current literature.


Introduction
Generally, high-order accurate discretizations suffer from spurious oscillations near high gradients due to Gibbs oscillations [1].The fifth-order accurate Weighted Essentially Nonoscillatory method was developed by Jiang and Shu (WENO-JS) to achieve high-order accuracy while avoiding Gibbs oscillations [2].The design of WENO-JS is based on a clever nonlinear modification of the fifth-order Linear Upwind (LUW5) discretization.The relatively wide five-point stencil of LUW5 can be subdivided into three small substencils.In WENO-JS, each smooth substencil is assigned a higher weight and vice versa.Hence, the constant weights adopted by LUW5 are replaced by nonlinear rational weights in WENO-JS.
A wide variety of WENO-JS improvements were developed, and a sample of the recent developments was provided in [3].On the contrary, theoretical studies of WENO-JS are relatively rare.Spectral analysis of WENO-JS is a theoretical methodology that quantifies the solution in terms of dissipation and dispersion [4][5][6][7].The spectral analysis methodologies used to analyse WENO-JS can be classified into two main categories: Geometric Stability Analysis (GSA) and Approximate Dispersion Relation (ADR).
The pioneering study of [5] adopted GSA to study the Linear Stability of WENO-JS and concluded that the forward Euler time integration method is linearly unstable.The results of [5] were revised by [6] who illustrated how linear stability can be satisfied because the discretization spectrum is equal to a finite set.However, the GSA provided by [5,6] adopted linearized approximations of WENO-JS.For instance, many results were derived by [5] using Taylor expansions in terms of the modified wave number (MWN).Hence, these results were limited to small MWN that describe smooth solutions.This is a major disadvantage because the main goal of WENO-JS is modelling non-smooth solutions coexistent with high MWN.Moreover, the GSA of [6] was mainly based on LUW5.Hence, both [5,6] did not identify the discrepancy between WENO-JS and LUW5 at high MWN.
The ADR is an alternative methodology based on discrete Fourier transform (DFT) and was originally developed by [4].Within ADR, an initial condition defined on an arbitrary grid is advanced numerically in time.Next, DFT is applied to the solution to obtain the spectrum.The ADR algorithm accounts for nonlinearity and is valid for arbitrary MWN.The ADR was used to analyse the cost efficiency of WENO-JS [4], and a wide application was presented in [8].The ADR was further developed by [7] and applied to analyse synthetic turbulent energy spectra [9,10].The design of the ADR algorithm requires adopting an arbitrarily small time step to minimize the effects of time integration.There are no definite criteria to select the numerical grid parameters within the ADR.In addition, neither the effects of the Courant number (CN) nor those of the numerical integration method were provided in [4,[7][8][9][10].
A cornerstone assumption that was adopted in the past GSA and ADR studies was ignoring the higher-order modes induced due to WENO-JS nonlinearity.However, theoretical validation of this assumption was not provided yet.This assumption will be denoted by mode isolation in this work.An exact methodology is proposed in this study to address the challenges mentioned above.A novel Nonlinear Spectral Method (NSM) is introduced to analyse WENO-JS.The benefits of the application of spectral methods to nonlinear problems are well explained in [11].The development of the novel NSM is based on a continuous presentation of the discrete solution.This subtle and key idea does not compromise accuracy and proves to be effective within the analysis.Unlike [5,6], NSM includes the full nonlinear rational weights of WENO-JS, and the results are valid for the full range of MWN.The NSM does not adopt any arbitrary grids, and the effects of CN and various time integration methods are described.A quantitative analysis of the mode isolation assumption is provided.To the author's best knowledge, these contributions are provided for the first time.

Numerical Discretization
The numerical algorithms studied in this work will be detailed.Consider the scalar conservation law given by subject to the suitable boundary and initial conditions.Here, u = uðx, tÞ stands for a conserved quantity, and f ðuÞ denotes the corresponding flux.Space and time domains are defined as Ω x = ½0, L and Ω t = ½0,∞Þ, respectively.To obtain a numerical solution of Equation (1) using the finite difference method, the space domain Ω x is discretized using the grid X given by X = x i : Here, h = L/N stands for the grid spacing, and N ∈ ℕ stands for the total number of grid subdivisions.The numerical solution is defined as 2.1.Fifth-Order Space Discretization.The fifth-order accurate space discretization algorithm is explained in numerous sources and will be presented briefly.The reader may refer to [1] for more details.The main task is to obtain a high-order finite difference discretization of the flux space derivative ∂f ðuÞ/∂x .To simplify notations, let f stand for f ðuÞ.The space derivative of f at x i will be written as ∂f ∂x The attention is shifted towards the calculation of f i−1/2 .This quantity is commonly termed as the numerical flux at point calculated based on the five-point stencil shown in Figure 1.This five-point stencil includes three smaller substencils as also illustrated in Figure 1.The three substencils are given by ð5Þ Various choices exist to calculate f i−1/2 .Consider the choice of f j i−1/2 that yields a third-order approximation given by ∂f ∂x Here, the superscript j corresponds to the substencil S j selected to calculate f j i−1/2 .The linear combination of f i that is used to calculate f j i−1/2 is given by The values of the constants c j p are provided explicitly below for each S j .
The three choices for f j i−1/2 can be combined to provide a fifth-order accurate flux using the weighted sum given by The selection of the weights C j will be detailed in the next subsections.
2.1.1.Fifth-Order Linear Upwind Discretization.Generally, f i are interpolated using a high-order polynomial.The following constants are defined [2,12]: The same values of d j are provided by [1] in the opposite order.The reason for this slight difference is that [1] listed the formulas to calculate the flux at point x i+1/2 instead of x i−1/2 .For the fifth-order Linear Upwind (LUW5) method, C j are given by Hence, the fifth-order accurate linear numerical flux f i−1/2 is given by The algorithm of LUW5 is linear because constant values are assigned to C j , and accurate results may be obtained in smooth regions.However, the high-order polynomial used in LUW5 may produce spurious oscillations near high gradients due to the Gibbs phenomenon.

Fifth-Order Weighted Essentially Nonoscillatory
Discretization.The fifth-order WENO-JS method is designed to overcome the disadvantages of LUW5 [1], and C j are designed to be higher for smooth S j and vice versa.Similar to [5,6], the current analysis will focus on WENO-JS.Smoothness is quantified using the smoothness indicators (β j ) based on a sub-stencil S j .The L 2 norm squared was used to calculate β j in [2] using the following formula: Here, q ðmÞ j is the m th derivative of the polynomial used to interpolate f at S j .This definition of β j in Equation ( 13) yields the following expressions: The modified weights are calculated as Here, ϵ is an arbitrary small number introduced to avoid division by zero.The normalized weights are given by Hence, for WENO-JS, the weights C j are given by and WENO-JS flux is given by Substituting f i−1/2 by f i−1/2 (given by Equation ( 12)) or f i−1/2 (given by Equation ( 18)) in Equation ( 4) results in LUW5 or WENO-JS space discretizations, respectively.

Numerical Time Integration. The operator
The operator Lðu i Þ encapsulates the space discretization algorithm.The time integration algorithm is based on the method of lines [1,4,5].A semidiscrete approximation of Equation ( 1) is given by To perform the numerical integration of Equation ( 20), the time domain Ω t is discretized using T given by The time step Δ is defined as Δ = t n+1 − t n .Numerical integration can be done using a variety of standard methods.The single-stage, first-order Euler forward solution of Equation (20) is written as follows [1]: Here, Þ are the solution vectors at t n and t n+1 , respectively.The symbol λ = Δ/h is introduced.And Equation ( 22) is rewritten as follows: For a wave speed of unity, λ is commonly termed the Courant number (CN).The first-order Euler forward method can be regarded as a member of the set of Runge-Kutta (RK) methods [13].Hence, the method defined by Equation ( 23) will be denoted by RK1.Higher-order multistage time integration algorithms can be implemented using the intermediate solution vectors u i ð1Þ and u i ð2Þ .The second-order explicit total variation diminishing Runge-Kutta (TVDRK2) solution of Equation (20) is given by The third-order explicit total variation diminishing Runge-Kutta (TVDRK3) solution of Equation ( 20) is given by The combination of WENO-JS space discretization with Equations ( 23), (24), and (25) will be denoted by WENO-RK1, WENO-TVDRK2, and WENO-TVDRK3, respectively.Similarly the combinations of LUW5 with Equations ( 23), (24), and (25) will be denoted by LUW5-RK1, LUW5-TVDRK2, and LUW5-TVDRK3, respectively.

Analysis
Similar to [4][5][6]9], the analysis will focus on the linear flux f ðuÞ = u.Hence, Equation ( 1) is rewritten as The space periodic boundary condition uð0, tÞ = uðL, tÞ is adopted.Consider the initial condition uðx, t = 0Þ = FðxÞ.The exact solution u e ðx, tÞ is obtained by the method of characteristics as follows: Hence, the exact solution corresponds to a righttravelling wave whose form is identical to the initial condition and speed equal to 1.

Spectral Analysis of Numerical Solutions. Any periodic intermediate solution u n
i is exactly presented as a summation of Fourier modes.
The total number of modes is equal to the grid size N.The modified wave number (MWN) termed as θ j is given by The operation defined by Equation ( 28) is termed the discrete Fourier transform (DFT).An extended explanation of DFT is provided in Appendix A. Generally, the range of θ j is θ j ∈ ½0, 2π.However, for any real-valued periodic function, the total number of unique modes is N/2 and the rest of the modes are obtained by Hermitian symmetry [14].Hence, the effective range of θ j is reduced to be θ j ∈ ½0, π as adopted by [4,8].The Hermitian symmetric property of DFT of real functions is derived in Appendix B.
The numerical solution u n+1 i of any numerical algorithm can be written as The spectral parameters G j and Φ j stand for the magnification factor and the phase shift, respectively.Generally, G j and Φ j are expressed as The effect of the numerical algorithm on each mode is fully described in terms of G j and Φ j .A value of G j > 1 or G j < 1 indicates whether mode j is magnified or damped, respectively.The significance of Φ j will be clarified as follows.Rewrite Equation (30) in terms of x i , and factor out θ j .
4 Journal of Applied Mathematics Each mode j in u n i can be regarded as a wave, that travels at a speed V j .Hence, u n+1 i is obtained by shifting u n i by a distance V j Δ.Hence, Equation (33) is rewritten in the form of the travelling wave solution: Comparing the arguments of the sin function in Equations ( 33) and (34), we get Hence, Φ j is proportional to V j .
3.1.1.Mode Isolation.Consider a monochromatic intermediate solution u n i that is described by a single mode [4] and defined as For this monochromatic case, let us denote all wave parameters with subscript j (e.g., θ j ) as the principal mode parameters.It will be useful to write explicitly the exact solution for this case using Equation (27).First, Equation (36) is rewritten as Next, the exact solution is obtained by replacing x i with x i − Δ and using the definition of Equation (2).
Finally, Equation (38) is rewritten in terms of the spectral parameters, introduced in Equation (30).
Upon comparison of Equations ( 38) and (39), the spectral parameters of the exact solution are given by Let us consider the numerical solution u n+1 i .Based on Equation (36), the expression given by Equation (30) is unaltered.However, the expressions given by Equations (31) and (32) are simplified to be As the accuracy of the numerical discretization is improved, the expressions provided by Equations ( 42) and (43) approach those provided by Equations ( 40) and (41), respectively.
So far, no limitations have been imposed on the presented analysis.However, it will be useful to introduce an assumption that was a cornerstone in [4][5][6].This assumption will be denoted by the Mode Isolation or MI for abbreviation.This assumption is stated as follows.
Assumption 1.Let u n i be a monochromatic intermediate discrete solution described by Equation (36).Then, u n+1 i resulting from Equation ( 23) is well approximated as In other words, when a monochromatic u n i is considered, the numerical solution u n+1 i can be equated to a single mode.Hence, all the summation terms in Equation ( 30) are equated to zero, except that corresponding to the principal mode.The validity of Assumption 1 will be thoroughly clarified in the next sections.

Linear von Neumann Analysis of LUW5.
In this subsection, the Linear von Neumann (VN) analysis will be presented by application to the LUW5-RK1 algorithm.The VN analysis is generally used to study finite domain problems with periodic boundary conditions [15].Adopting Equation ( 11) yields a linear scheme, and the operator Alternatively, Equation ( 45) is rewritten in the following compact form.
Here, D r are constants whose values are listed explicitly in Equation (45).The numerical solution of LUW5-RK1 based on Equation ( 23) can be written as The coefficients of u i+r on the right-hand side of Equation (47) are constants.Hence, Equation (47) is linear in terms of u n i+r , and the analysis can be simplified as follows.The VN analysis considers a discrete Fourier Mode [16] that corresponds to the intermediate solution given by Equation (36).Although one mode is considered out of the multiple modes in Equation ( 28), the analysis is still general [17].Also, adopting the complex notation is allowed, and Equation (36) can be rewritten in terms of i = ffiffiffiffiffi ffi −1 p as follows: The expression of Equation ( 48) is substituted into Equation (47), and we get To simplify Equation ( 49), the term a j e iðθ j i+ϕ j Þ is factored out to obtain Here, Hðθ j , λÞ = 1 + λ∑ 2 r=−3 D r e iθ j r , and the second equality is obtained by substitution by Equation (48).The spectral parameters are calculated as Hence, the numerical solution of LUW5-RK1 is written as The following points should be emphasized regarding the application of the VN method to LUW5-RK1 and other linear discretizations.
The expressions for G j and Φ j provided by Equations ( 51) and (52) are exactly valid for arbitrary values of h and Δ.
(i) As illustrated by Equation (50), Assumption 1 is exactly satisfied by LUW5-RK1 and all linear discretizations (ii) Referring to Equations ( 51) and ( 52), the spectral parameters G j and Φ j are independent of a j and ϕ j .Hence, Equations ( 42) and ( 43) are simplified to be 3.2.1.Extension to High-Order Time Integration.The analysis based on VN can be extended to high-order multistage time integration methods.It should be noted that the algorithms described in Equations ( 24) and (25) adopt Equation ( 23) as a building block.Hence, Hðθ j , λÞ will be used frequently.Since LUW5 is linear, the following simplification of expressions including L is justified: The symbol A in Equation ( 56) stands for any expression that does not depend explicitly on i.The analysis will be based on adopting u n i given by Equation (48).To analyse LUW5-TVDRK2, Equation (24) is rewritten as Hence, for LUW5-TVDRK2, the spectral parameters are defined as Similarly, LUW5-TVDRK3 is analysed using Equation (25).The derivation follows the same ideas listed in Equation (57).To avoid distraction, the derivation is provided in a relatively brief format.
Hence, for LUW5-TVDRK3, the spectral parameters are defined as 3.3.Nonlinear Spectral Analysis of WENO-JS.It should be noted that the numerical solution of Equation (26) using WENO-JS is nonlinear.The reason is that C j used in WENO-JS are not constants and are based on rational expressions in terms of u n i .Hence, the steps explained in Section 3.2 cannot be applied directly to WENO-JS.Still, the nonlinear features of WENO-JS could be analysed using a slight modification of the VN method.This modification is 6 Journal of Applied Mathematics based on a continuous presentation of WENO-JS.This presentation's main benefit is clarifying two main solution properties: periodicity and homogeneity.As explained in Section 3.2, the VN method is based on analysing Fourier modes for periodic finite domains.It will be shortly illustrated how the solution periodicity will be used to develop a nonlinear analysis that enjoys several advantages of the VN method.For instance, the grid parameters h and N will be eliminated and the analysis will depend only on λ and θ j .The homogeneity of WENO-JS weights and flux will eliminate the dependence of the numerical solution on a j if the MI is adopted.It could be shown that the periodicity and the homogeneity properties are satisfied by several other nonlinear numerical algorithms.However, only WENO-JS will be considered in the current work.

Continuous Presentation of WENO-JS.
The algorithm of WENO-JS is designed originally to deal with discrete quantities.However, WENO-JS algorithm steps can be applied to continuous quantities as follows.The continuous version of β j is obtained by substituting with the linear flux f ðuÞ = u into Equation ( 14) and replacing f i−k by uðx − khÞ, where k ∈ f−1, 0, 1, 2, 3g.Hence, based on Equation ( 14) the continuous version of β j is given by The superscript c in β c j ðx ; hÞ stands for continuous.The reader should note the following: (i) Although t was skipped in Equation (61), time dependence is still included.In fact uðx − khÞ = uðx − kh, tÞ.This symbolic simplification is justified since Lðu i Þ operates only in space (ii) β c j ðx ; hÞ is defined for a continuous range of x ∈ ½0, L. This is in contrast to β j defined for the discrete points in S j defined by Equation (5) (iii) β c j ðx ; hÞ is equal to a second-order homogeneous polynomial [18] in terms of uðx − khÞ.Hence, β c j ðx ; hÞ can be regarded as a composite homogeneous function Based on the discussion above, the following equation is introduced: Here, ∘ stands for the composite function, and g β j stands for the polynomials defined in the right-hand side of Equation (61).Since g β j ðuÞ is a second-order homogeneous function, the following equation is satisfied: Lemma 1.Let uðxÞ be a space periodic function with period P ∈ ℝ.Then, β c j ðx ; hÞ is also space periodic with period P.
Proof.By definition, Substituting by x + P into Equation (62), we get The second equality is obtained based on Equation (64).Hence, β c j ðx ; hÞ is also periodic with period P. To proceed further, Equation ( 16) is rewritten as Following the same direction, a continuous version of ω j is defined as 7 Journal of Applied Mathematics A composite function form can be also used to describe ω c j ðx ; hÞ as follows.Substituting Equation (62) into Equation (67), we get Corollary 2. Let uðxÞ be a space periodic function with period P ∈ ℝ.Then, ω c j ðx ; hÞ is also space periodic with period P.
Corollary 2 can be deduced using the same steps used to prove Lemma 1.
The value of ϵ should be proportional to the machine precision [12].Hence, ϵ will be neglected and assigned to zero within the current exact analysis.Consequently, the expression for g ω j ða × uÞ is written based on Equation (68) as The second equality in Equation ( 69) is based on Equation (63).Hence, g ω j ðuÞ is a zero-order homogeneous function as illustrated by Equation (69) that.A continuous version of f j i−1/2 as defined by Equation ( 7) is given by Similarly, a continuous version of f i−1/2 defined by Equation ( 18) is given by The flux Fc ðx ; hÞ can be regarded as the composite function given by To clarify homogeneity, the flux Fc is applied to a × u as follows: The second equality is based on Equations ( 68) and (70).Hence, g F ðuÞ is a first-order homogeneous function.

Corollary 3.
Let uðxÞ be a space periodic function with period P ∈ ℝ.Then, Fc ðx ; hÞ is also space periodic with period P.
Corollary 3 can be deduced using the same steps used to prove Lemma 1 and Corollary 2.
Based on Fc ðx ; hÞ, a continuous operator L c ðuðxÞÞ is defined as Corollary 4. Let uðxÞ be a space periodic function with period P ∈ ℝ.Then, L c ðuðxÞÞ based on WENO-JS is also space periodic with period P.
The proof is straightforward and is based on Corollary 3. It can be illustrated that L c ðuðxÞÞ is a first-order homogeneous composite function as follows: The first equality is justified by Equation (73).
The main theorem will be presented through the analysis of the relatively simple WENO-RK1 algorithm.
Theorem 5. Consider the initial value problem of Equation (26), subject to space periodic boundary conditions.Let the initial condition uðx, 0Þ be a space periodic function with period P ∈ ℝ.Then, the numerical solution based on WENO-JS space discretization and RK1 time integration is also space periodic with period P.
Proof.The proof proceeds by induction.Let u n ðxÞ = uðx, t n Þ be a space periodic function with period P. Applying Equation (23), we get The periodicity of λL c ðu n ðxÞÞ is clarified by Corollary 4, and u n ðxÞ is periodic by assumption.Hence, the whole right-hand side of Equation ( 76) is periodic with period P. The proof is straightforward for the base case since u 0 ðxÞ = u ðx, 0Þ is periodic with period P by assumption.Corollary 6.Let u n ðxÞ = uðx, t n Þ be a space periodic function with period P. Consider the numerical solution u n+1 ðxÞ based on WENO-JS space discretization and RK1 time integration.If u n+1 ðxÞ is C 0 continuous, then it can be expressed using the following Fourier series: The coefficients ã0 , ãα , bα , and cα are the Fourier expansion coefficients defined as follows: Proof.By definition, Hence, from Theorem 5, Consequently, u n+1 ðxÞ can be expressed as the infinite trigonometric series provided by Equation (77).The series expression is provided by Equations (3.7, 3.10) in [19], and provided by Equation (77).The series coefficients ã0 , ãα , bα , and cα are provided by Equations (3.5, 3.6, 3.9) in [19] and given by Equations ( 78)-(81).The condition of C 0 continuity is included to assure the series convergence.For more details, the reader may refer to Theorem 4.3 in [19].

Fourier Analysis of a Monochromatic Solution.
Consider a monochromatic u n i with a single period P j and MWN θ j .As illustrated in Theorem 5, the subsequent solutions will be also exactly periodic with the same period P j .The variable i c ∈ R is defined as a continuous version of index i to be Both variables x and i c are interchangeable.Hence, ũn ði c Þ is defined as Also, the numerical solution u n+1 ðxÞ defined in Equation (76) can be written as Using the continuous presentation, Equation (37) is rewritten in terms of x instead of x i as Also, Equation (37) is rewritten in terms of i c as The periodicity of u n ðxÞ and ũn ði c Þ are defined as The periods are defined as The main analysis step is substituting Equation (88) into Equation (76).Based on Corollary 6, u n+1 ðxÞ can be expressed as a series expansion of Equation (77).In principle, the coefficients ã0 , ãα , bα , and cα can be obtained from Equations (78), (79), (80), and (81).However, equivalent formulas can be written in terms of i c as follows: 9 Journal of Applied Mathematics Although u n ðxÞ and ũn ði c Þ refer to the same quantity, the latter is in Fourier analysis.The specific advantage of ũn ði c Þ as defined in Equation ( 88) is that h is eliminated.The solution behaviour is fully clarified considering only i c ∈ f0, P j g, instead of considering the whole domain i ∈ f0, Ng.
Although u n ðxÞ in Equation ( 87) is monochromatic, the numerical solution u n+1 ðxÞ includes multiple modes as described in Equation (77).In other words, a monochromatic initial condition induces high-order modes.Assumption 1 will be revisited within the continuous presentation provided.As implied by MI, zero and high order modes can be neglected and Equation (77) can be rewritten as follows: The spectral parameters are defined based on Fourier expansion coefficients as follows: Using G j and Φ j , Equation ( 96) is rewritten as The parameter ϕ j describes a shift in the waveform of u n ðxÞ.This shift does not affect neither c1 nor c−1 [19].In addition, c1 and c−1 depend linearly on a j as illustrated by Equation (75).Consequently, the values a j = 1 and ϕ j = 0 can be adopted without loss of generality.Hence, G j and Φ j can be expressed using Equations ( 54) and (55), respectively.
The analysis described in Equations ( 97) and (98) will be denoted by the Nonlinear Spectral Method, abbreviated as NSM.Upon introduction of MI, the focus is limited to calculating only the principal mode spectral parameters, namely, ã1 and b1 .Still, it should be emphasized that the error introduced by MI does not affect the accuracy of NSM.In fact the principal mode parameters ã1 and b1 are still recovered exactly.The following features of the NSM should be noted: (i) The results of the NSM are valid for any h and Δ similar to VN, and the computational effects are eliminated.Hence, the limitations of the ADR developed by [4] are avoided (ii) The expressions resulting from plugging Equation (88) into Equation (76) are mightily lengthy.Hence, this work uses the computer algebra system MAXIMA [20] (iii) The integrands of Equations ( 93) and (94) include rational expressions in terms of trigonometric functions.This is attributed to the form of the nonlinear weights of WENO-JS.Still, the definite integrals can be calculated up to arbitrary accuracy using numerical quadrature.The reader could refer to [21] or other standard textbooks for the details of numerical integration using the trapezoidal approximation (iv) The nonlinear nature of the WENO-JS algorithm should be well captured by the present analysis.
The rational expressions of ω j are adopted, unlike [6] where ω j were approximated as constants (v) Past theoretical studies [5,22] adopted the assumption of small MWN (θ j ) because Taylor expansion was used.Such limitations are avoided, and the NSM is valid for arbitrary values of θ j 3.3.3.Mode Isolation Error.The error introduced by Assumption 1 will be quantified.Consider the well-known Parseval identity: It is straightforward to verify that Hence, the percentage error normalized with respect to kũ n+1 ði c Þk 2 is defined as The error Eðθ j , λÞ is proportional to the magnitude of zero and higher-order modes truncated by ũa n+1 ði c Þ.The reader is 10 Journal of Applied Mathematics reminded that MI is exactly satisfied by LUW5.Hence, a considerable discrepancy may exist between the results of WENO-JS and LUW5 due to Eðθ j , λÞ, even though the of both methods may be similar in some cases.
11 Journal of Applied Mathematics The algorithm described in Section 3.3.2 is modified by replacing Equation (76) with Equations ( 103) and ( 104) to analyse WENO-TVDRK2 and WENO-TVDRK3, respectively.

Geometric Stability Analysis.
The spectral methods developed by [5,6] termed the Geometric Stability Analysis (GSA) will be briefly presented.The spectral properties of the spatial discretization operator are defined in terms of z ðθ j Þ given by Here, u n i is given by Equation (48) because MI is a cornerstone in GSA [5,6].The continuous spectrum S is given by Generally, numerical time integration methods are characterized in terms of ẑ given by The amplification factor gðẑÞ for a numerical scheme is given by The relation between gðẑÞ and the spectral parameters in the current work is given by For any numerical time integration method, the Linear Stability Domain and its boundary are given by The Geometric Linear Stability Limit (GLSL) is defined as the largest λ such that λS ⊆ S t : ð111Þ It was illustrated by [5] that the GLSL exists for neither WENO-RK1 nor WENO-RK2, and both were concluded to be linearly unstable.However, the grid X given by Equation (2) was used in [6] to define a discrete spectrum S that consists of the set of eigenvalues given by Hence, Equation (111) should be replaced by It was illustrated by [6] that Equation ( 113) can be satisfied by both WENO-RK1 and WENO-RK2.Given a specific value of λ, define the intersection point between λS and ∂S t to be ẑ * ðλÞ.The nonzero eigenvalue of S that is closest to the origin −zðθ j = 2πhÞ is equated to ẑ * ðλÞ/λ = −zðθ * j ðλÞÞ.Hence, the linear stability condition can be given by For LUW5-RK1, the following formula of θ * j ðλÞ that is valid for small MWN was derived by [6]  The function zðθ j Þ plays a key role in the GSA.However, linearised procedures were used to calculate zðθ j Þ in [5,6].For instance, the GSA of [6] was mainly based on Equation (50) which describes LUW5 instead of WENO-JS.The GSA provided by [5] was mainly based on Taylor expansion and was limited to small θ j .The proposed NSM can be used to calculate zðθ j Þ exactly as follows.Define ãL 1 and bL 1 as Hence, zðθ j Þ is given by 3.5.Analysis Based on Numerical Experiments.The results provided above can be further verified using numerical experiments.The methodology follows closely the procedure outlined by [4].Also, complete derivations are provided in Appendix A. For given values of N, L, and λ, a grid x i is adopted based on Equation (2).A value of θ j is selected, where j ∈ f0, 1, ⋯, N/2g.A monochromatic intermediate condition u n i = u n ðx i Þ is calculated based on Equation (87).The solution u n+1 i is obtained by advancing one time step using one of Equations ( 23), (24), and (25).Finally, Fast Fourier transform (FFT) procedure is applied to u n+1 i to get the spectrum U n+1 w = Fðu n+1 i Þ, where w ∈ f0, 1, ⋯, N/2g.According to the MI assumption, only U n+1 j is considered and other modes are neglected.Hence, the magnification factor and the phase shift are calculated as The spectral parameters are calculated using the FFT of the solution obtained using a specific grid size N. Hence, the procedure will be sensitive to the machine's precision.A unity domain size L = 1:0 was adopted for this paper's numerical experiments.

Significance of the Nonlinear Analysis
4.1.1.Linear and Nonlinear Weights.To clarify the significance of the NSM, the weights of both LUW5 (d j ) and WENO-JS (ω j ) were plotted.The curves are provided in Figure 2 for four values of θ j = 0:1π, 0:3π, 0:5π, and 0:9π, covering small and large MWN.A low MWN corresponds to smooth solutions at which the performance of LUW5 and WENO-JS should be close.Hence, for θ j = 0:1π, the difference between d j and ω j is small.As the value of MWN increases the solution smoothness decreases and the influence of nonlinear WENO-JS weights should be stronger.For instance, at θ j = 0:5π, we get P j = 4 based on Equation (91), and a complete wave is located within 4 grid spacings h.Consequently, for large θ j , the difference between d j and ω j cannot be ignored.This can be easily observed in Figure 2 for θ j = 0:3π, 0:5π, 0:9π.The value of Eðθ j , λÞ vanishes as θ j ⟶ 0. This corresponds to the trivial case of a constant u n i , where Lðu n i Þ = 0.
(i) On the other hand, Eðθ j , λÞ also vanishes as θ j ⟶ π.As already illustrated, this corresponds to the maximum θ j sampled by the grid.Hence, the solution is described exactly using a single mode and the MI is exactly satisfied (ii) The value of Eðθ j , λÞ is almost negligible for the interval θ j < π/2.For this range of θ j , modes extend over more than 5 grid points and the solution will be relatively smooth.Hence, WENO-JS accuracy should be relatively higher and MI is well satisfied.This is the opposite for θ j > π/2, where modes extend over less than 5 grid points and u n i exhibits strong variation -0.5 0 0 0.5 1 1.5 2 2.5 3 (iii) As λ increases, Eðθ j , λÞ increases and MI becomes less satisfied.This is expected since the numerical accuracy is generally inversely proportional to λ (iv) The value of Eðθ j , λÞ is generally similar and less than 1% for the three methods WENO-RK1, WENO-TVDRK2, and WENO-TVDRK3 at λ = 0:025,0:1.However, at λ = 0:4, the value of Eðθ j , λÞ is considerably higher for WENO-RK1 than both WENO-TVDRK2 and WENO-TVDRK3 Based on the above discussion, it is expected that the discrepancy between WENO-JS and LUW5 results should be high for solutions including modes of θ j > π/2 and using λ ≥ 0:4.

Numerical Solution Spectra.
The spectra of WENO-JS are plotted for λ = 0:025,0:1, and 0:4, and for the three time integration methods RK1, TVDRK2, and TVDRK3 in Figure 4. Hence, the effects of both λ and the time integration methods are clarified.For each value of λ, three spectra are provided based on the following: the NSM explained in Section 3.3.2,FFT numerical experiments explained in Section 3.5, and the exact solution spectrum given by Equations ( 40) and (41).The results of both FFT and NSM are almost identical, except for some minor discrepancies.Hence, the validity of NSM at the full range of θ j ∈ ð0, π is illustrated.The deviation between the spectra of WENO-JS and the exact solution is proportional to λ.This inverse proportionality between accuracy and λ is expected.Most of the spectra of WENO-JS are close to the spectra of the exact solution at θ j ∈ ð0, 1Þ.However, G j of WENO-RK1 is considerably high for λ = 0:4 at θ j ∈ ð0, 1Þ, as shown in Figure 4(a).
Due to the current limited scope, the GSA will be applied to RK1 only.The results of the GSA of WENO-RK1 and LUW5-RK1 are shown in Figure 5.In Figure 5(a), the red unit circle represents ∂S t of RK1 [13], and the blue and black curves represent S of LUW5 and WENO-JS, respectively.Both spectra are calculated exactly using Equations (50) and (117) for LUW5 and WENO-JS, respectively.The full range of the MWN θ j ∈ ð0, π is included in Figure 5(a).For small θ j in the neighbourhood of the imaginary axis, WENO-JS and LUW5 have very close S.However, S of LUW5 is closer to the imaginary axis than that  115) to be ≈0:135π and 0:132π, respectively.Hence, the relation with the GSA of [6] is illustrated.Also, it is shown in Figure 5(b) that θ * j ð0:001Þ ≈ 0:075π for WENO-RK1.Hence, the intersection between λS and S t is higher for WENO-RK1 than that for LUW5-RK1 at λ = 0:001.Also, the grid spacing is restricted to h ≥ 0:075/2 and h ≥ 0:135/2 for WENO-RK1 and LUW5-RK1, respectively.Hence, WENO-RK1 admits a wider range of h compared to LUW5-RK1.
The spectra are plotted for WENO-TVDRK3 and LUW5-TVDRK3 using moderate and high values of λ in Figure 6.Generally, the magnitudes of G j and Φ j of WENO-TVDRK3 are considerably lower than those of LUW5-TVDRK3 for θ j ∈ ðπ/2, πÞ.The low values of G j and jΦ j j correspond to more dissipation and deceleration of high MWN modes.This is expected because WENO-JS is designed to filter out high MWN modes.Specifically, the weighting algorithm of WENO-JS is focused on the solution behaviour within the five-point stencil.The width of this stencil is 4h, and it is covered by the MWN modes of θ j ≥ π/2.Hence, the effect of WENO-JS should be stronger on the MWN modes of θ j ≥ π/2.The spectra of WENO-TVDRK3 and LUW5-TVDRK3 show somehow similar trends for λ = 0:3,0:7.However, the spectra differ considerably at higher λ.Specifically, two small peaks can be observed for WENO-TVDRK3 at λ = 0:9 in the range θ j ∈ ð1:5,2:5Þ as shown in Figure 6(e).These two peaks merge into a single spike at λ = 1:1 located at the high MWN θ j ≈ 2:1 as shown in Figure 6(g).This spike contradicts the main requirement of WENO-JS, which is damping high MWN modes.In addition, the effect of this spike will extend to other modes due to the high Eðθ j , λÞ at θ j ≈ 2:1 (Figure 3).Hence, the convergence of WENO-TVDRK3 solutions at high λ will deteriorate.Figure 8: Spectra of the exact (initial) solution --and numerical solutionsusing TVDRK3 time integration based on Equation (120).The numerical parameters are provided in the text.The corresponding curves for λ = 0:3,0:7,0:9, and 1:1 are plotted in magenta, black, red, and blue colours, respectively.to elucidate any long-term effects.Three experiments are performed adopting the initial conditions uðx, t = 0Þ = FðxÞ given by The numerical and exact solutions based on Equations (119), (120), and (121) are provided in Figures 7-9, respectively.The results of LUW5-TVDRK3 and WENO-TVDRK3 are provided in the upper and lower parts of Figures 7-9, respectively.The solutions uðx, t = 20Þ and the magnification factors are provided in the left and right parts of Figures 7-9, respectively.Based on Equation (27), the solutions plotted at t = 20 correspond to traversing outside and inside the domain and returning to the original position 20 times.Hence, the exact and the initial conditions are coincident and plotted using the same dashed line.
The initial condition of Equation (119) includes a sharp discontinuity and is selected to clarify the performance of WENO-JS for nonsmooth solutions.For the exact solution, G j approximates a Sinc function with multiple lobes.This behaviour is expected and is explained in [14] and other standard textbooks.
Considering LUW5-TVDRK3, G j approximates the first three lobes in the range θ j ∈ ½0, 1 as shown in Figure 7(b).However, G j almost vanishes for high MWN in the range θ j ∈ ð1, π and the spectrum is almost truncated.This truncation leads to an oscillatory behaviour and Gibbs oscillations can be observed in the solution plotted in Figure 7(a).This truncated spectrum is expected since G j was close to 1 in the range θ j ∈ ½0, 1 for LUW5-TVDRK3 spectra presented in Figure 6.The behaviours of LUW5-TVDRK3 solutions are almost similar for λ = 0:3,0:7,0:9,1:1.
For WENO-TVDRK3, the solution behaviour is strongly affected by λ.Considering low and moderate CN (λ = 0:3,0:7), multiple lobes are included in G j and cover the whole MWN range as shown in Figure 7(d).This distribution of G j approximates the product of two Sinc functions that describes a relatively smeared but non-oscillatory solution.Hence, WENO-TVDRK3 solutions at λ = 0:3,0:7 approximate the exact solution avoiding Gibbs oscillations as shown in Figure 7(c).Inverse proportionality between G j and θ j is observed in Figure 7(c).However, G j does not vanish at high MWN in the range θ j ∈ ð1, π.This behaviour is attributed to the nonlinear nature of WENO-JS, where low MWN modes contribute to high MWN modes as illustrated in Section 3.3.2and quantified by Eðθ j , λÞ. Figure 9: Spectra of the exact (initial) solution --and numerical solutionsusing TVDRK3 time integration based on Equation (121).The numerical parameters are provided in the text.The corresponding curves for λ = 0:3,0:7,0:9, and 1:1 are plotted in magenta, black, red, and blue colours, respectively.18 Journal of Applied Mathematics The distribution of G j for WENO-TVDRK3 considerably deteriorates at the high CN (λ = 1:1).The Sinc distribution disappears and a strong peak is observed at θ j ≈ 2:1 as shown in Figure 7(d).This peak is expected and corresponds to the spike that is clarified in Figure 6(g).As already mentioned in Section 4.2, this spike has an adverse effect on the solution.The spike tends to magnify high MWN modes leading to nonsmooth and irregular solutions.On the other hand, the WENO-JS mechanism tends to smooth the solution and remove the irregularity.As a result, the solution of WENO-TVDRK3 for λ = 1:1 shown in Figure 7(c) is overdamped and infested by small yet distinguishable spurious oscillations.
The detection of this overdamping phenomenon is generally difficult using numerical experiments.For instance, consider the solutions based on Equation (120) plotted in Figure 8.This initial condition, given by Equation ( 120), includes only one low MWN mode (θ j = π/64).Both WENO-TVDRK3 and UW5-TVDRK3 solutions are almost indistinguishable, and the effect of λ disappears.The reason for this is overdamping occurs due to the influence of high MWN modes.A numerical experiment was done by [5] at λ = 1:4 using a low MWN initial condition.Hence, overdamping was hidden, and it was erroneously concluded that the performance of WENO-TVDRK3 is generally satisfactory at this high CN.
Finally, the sensitivity of WENO-TVDRK3 to high wave numbers is illustrated by considering the initial condition of Equation ( 121).This condition consists of a high amplitude low MWN mode (θ j = π/64), and a low amplitude high MWN mode (θ j ≈ 2:2).The amplitude of the high MWN mode (1/100) is almost negligible compared to that of the low MWN mode.The influence of θ j ≈ 2:2 mode is negligible for UW5-TVDRK3 solution plotted in Figure 9(a), and it is effectively filtered in UW5-TVDRK3 spectra plotted in Figure 9(b).However, the solution of WENO-TVDRK3 suffers overdamping at λ = 1:1 as illustrated in Figure 9(c).A small yet distinguishable peak located at θ j ≈ 2 can be observed in the WENO-TVDRK3 spectrum for λ = 1:1 in Figure 9(d).This peak corresponds to the spike that is clarified in Figure 6(g).

Conclusions
The nonlinear Weighted Essentially Nonoscillatory method developed by Jiang and Shu (WENO-JS) is analysed to clarify its spectral properties.An exact Nonlinear Spectral Method (NSM) is developed based on the Fourier expansion of a continuous presentation of WENO-JS solutions.The nonlinear rational formula of WENO-JS weights is fully included.Unlike other theoretical methods, the NSM is valid for arbitrary modified wave numbers (MWN).The mode isolation (MI) assumption, which was used extensively in the literature, is quantified for the first time, and its impact on the performance of WENO-JS is explained.The MI error increases for high MWN and moderate Courant numbers (CN), and the discrepancy between the results of WENO-JS and those of the fifth-order Linear Upwind (LUW5) discretization cannot be neglected.
The NSM can reproduce the spectra of numerical experiments while being independent of any grid parameters.The effects of both the CN and the numerical time integration methods are elucidated.At low CN, the performance of the combination of WENO-JS with the forward Euler time integration method is better than that of LUW5.At high CN, the combination of WENO-JS with the third-order variation diminishing Runge-Kutta method suffers extra damping compared to LUW5.This peculiar phenomenon was overlooked in past studies because it occurs due to a spike located at the high MWN spectra, and its detection via numerical experiments is difficult.The results of the current study should be extended to include other time integration algorithms and nonlinear discretizations.

Appendix A Discrete Fourier Transform
Consider the periodic function uðxÞ of period L ∈ ℝ.Samples of uðxÞ are taken at N ∈ ℕ points x i , where i ∈ f0, 1, ⋯, N − 1g.The samples are obtained at equally spaced points x i .
x i = i L N = ih: ðA:1Þ Hence, the sequence u i = uðx i Þ is obtained.The sequence can be written as the following sum [14]: U w e i 2wπ/L ð Þ x i : ðA:2Þ Here, i = ffiffiffiffiffi ffi −1 p .The coefficients U w are termed as the discrete Fourier transform or DFT as an abbreviation.The operator F is used to define DFT.This is defined as u i e −i 2wπ/L ð Þ x i : ðA:3Þ Similarly, the inverse operation of DFT can be defined using the operator notation: For a better presentation, the wave number is defined as u i e −iθ w i : ðA:8Þ

B DFT of a Monochromatic Wave
Consider the real monochromatic data series, with an amplitude A, a wave number k j , and an arbitrary phase shift ϕ.
The discrete Fourier transform is obtained as e iϕ e i k j −k w ð Þ x i − e −iϕ e −i k j +k w ð Þ x i , e iϕ e i j−w Following the same steps, it can be easily deduced that that U n = 0 for n∈fj, N − jg.These results can be verified by substitution by Equations (A.12) and (A.Finally, referring to Equations (A.12) and (A.14), it is deduced that ðA:16Þ The conclusions are summarized in the following points: (i) The amplitude of mode j is recovered from the DFT coefficient U j as A = jð2i/NÞU j j (ii) The phase shift of mode j is recovered from the DFT coefficient U j as ϕ = arg ðð2i/NÞU j Þ (iii) The range of θ w as defined in Equation (A.6) is θ w ∈ ½0, 2π.However, the effective range of θ w is reduced to θ w ∈ ½0, π due to the property of Hermitian Symmetry, as illustrated by Equations (A. 16) and (A.17) (iv) The results above were derived for a monochromatic signal.However, due to the linearity of DFT [14], the results still hold for any periodic signal.Consider the general case: A j sin θ j i + ϕ j : ðA:17Þ The DFT can be obtained as shown in Equation (A.3).Hence, the amplitude and phase shift of mode j can be recovered back from the DFT as follows: ðA:18Þ As mentioned in [6], N can be assumed even without loss of generality.Based on the above discussion, a discrete mesh yields only a finite number of distinct modes.These 20 Journal of Applied Mathematics modes can be represented using k j = 2jπ/L, or θ j = 2jπ/N where j ∈ f0, 1, ⋯, N/2g.

Figure 1 :
Figure 1: The full five-point stencil used to calculate f i−1/2 and the substencils S j where j ∈ f0, 1, 2g.