{'fraction_non_alphanumeric': 0.05195038110895232, 'fraction_numerical': 0.014168285756986997, 'mean_word_length': 4.297862232779098, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Stability analysis of the numerical Method of characteristics applied to a class of energy-preserving systems. Part I: Periodic boundary conditions

We study numerical (in)stability of the Method of characteristics (MoC) applied to a system of non-dissipative hyperbolic partial differential equations (PDEs) with periodic boundary conditions. We consider three different solvers along the characteristics: simple Euler (SE), modified Euler (ME), and Leap-frog (LF). The two former solvers are well known to exhibit a mild, but unconditional, numerical instability for non-dissipative ordinary differential equations (ODEs). They are found to have a similar (or stronger, for the MoC-ME) instability when applied to non-dissipative PDEs. On the other hand, the LF solver is known to be stable when applied to non-dissipative ODEs. However, when applied to non-dissipative PDEs within the MoC framework, it was found to have by far the strongest instability among all three solvers. We also comment on the use of the fourth-order Runge--Kutta solver within the MoC framework.


Introduction
In this series of two papers we address numerical stability of the Method of characteristics (MoC) applied to a class of non-dissipative hyperbolic partial differential equations (PDEs) in 1 time + 1 space dimension. For reference purposes, we will now present the idea of the MoC using the following system as an example: where f 1,2 are some functions. In the MoC, each of the equations in (1) is transformed to an ordinary differential equation (ODE) in its respective variable, η ± = x ∓ t, and then is solved by an ODE numerical solver.
The MoC is widely used and is described in most textbooks on numerical solution of PDEs.
It is, therefore, quite surprising that its stability has not been investigated in as much detail as that of most other numerical methods. For example, in [1], it was considered by the von Neumann analysis (which implies periodic boundary conditions (BC)) for a simple scalar model with the ODE solver being the simple Euler method. Clearly, that analysis could be straightforwardly generalized for a vector model where A is any constant diagonalizable matrix with real eigenvalues; an analysis for a similar model in three spatial dimensions was presented in [2]. However, we have been unable to find a systematic stability analysis -even for periodic BC -for a system of the form where A, B are constant matrices. A partial exception is a conference paper [3], where a model with a somewhat more complicated right-hand side (rhs), describing a specific engineering application, was studied by the von Neumann analysis. The focus of that paper was on the examination of the effect of various parameters of that particular model rather than on the impact on stability of the MoC by the ODE solver used. (In fact, the ODE solver used in [3] was the simple Euler method.) In this paper we will examine the effect of the ODE solver on the stability of the MoC applied to the model problem (3) with periodic BC. We will further limit the scope of the problem as follows.
First, and most importantly, we will consider only non-dissipative and stable systems. This implies that matrix B possesses only imaginary eigenvalues. Second, we will consider only three ODE solvers: simple Euler (SE), modified Euler (ME), and the Leap-frog (LF). Let us note that for a non-dissipative stable ODE, the numerical solutions obtained by the SE and ME are known to be mildly, but unconditionally, unstable, with the growth rates of the numerical instability being O(∆t) and O(∆t 3 ), respectively; see, e.g., Secs. 3 and 4 below. On the contrary, the LF solver is known to quasi-preserve (i.e., not shift systematically) at least some of the conserved quantities of non-dissipative ODEs.
Our third assumption is mostly cosmetic and does not affect the methodology of our analysis.
Namely, we will assume that the constant matrix A in (3) can be diagonalized into the form where c 1,2 are scalars, I N is the N ×N identity matrix, and N is the dimension of vector u. Without loss of generality (i.e., by a simple change of variables) one can set We will also let N 1 = N 2 = N/2. Particular conclusions of the forthcoming analysis may change if the diagonalization of matrix A is different from (4a). However, the methodology of the analysis will not be affected.
Our analysis reveals two facts which, to our knowledge, have not been previously pointed out.
First, in contrast to the situation with most numerical methods, the most numerically unstable Fourier harmonics may occur not at the edges of the spectrum but in the "middle" (see the footnote in Section 4) of it. In fact, the (in)stability of the highest Fourier harmonics is the same as that of the lowest ones. This fact easily follows from the foregoing analysis and holds for the MoC employing any ODE solver as long as matrix A has the form (4).
Second, and quite unexpectedly, the LF, which outperforms the SE and ME when applied to non-dissipative ODEs, performs much worse than those methods when applied to non-dissipative stable PDEs within the MoC framework.
The main part of this work is organized as follows. In Section 2 we specify the physical model to be considered below. Our analysis will be developed for the linearized version of that model. In Section 3 we present the von Neumann analysis and its verification by direct numerical simulations of the PDE for the case when the MoC employs the SE solver. In Sections 4 and 5 we repeat steps of Section 3 for the cases where the ODE solvers are the ME and LF, respectively.
It should be noted that in all three cases, the von Neumann analysis leads one to 4 × 4 generalized eigenvalue problems which involve both matrices B and A (or, as per (4a), Σ). The largest-magnitude eigenvalue, which determines the stability of the numerical method, of those problems cannot be analytically found or even related to the individual eigenvalues of the "physical" matrices B and A. However, for given B and A, that eigenvalue can be easily and quickly found numerically by standard built-in commands in software like Matlab, and thus we use this semi-analytical approach.
In Section 6 we comment on the case when the ODE solver is the fourth-order classical Runge-Kutta (cRK) method. We will not, however, analyze it in any detail because, as we will demonstrate, there is an uncertainty about using the cRK solver within the MoC framework. A detailed investigation of this issue is outside the scope of this paper. In Section 7 we present our conclusions and demonstrate the validity of our analysis for a broader class of systems than the particular system considered in Sections 2-5. First, in Appendix A we apply our analysis of the MoC-LF scheme to a class of models of which the model considered in the main text is a special case. These models still result in linearized equations with constant coefficients. Then, in Appendix B, we present a different model with a spatially localized solution and demonstrate by direct numerical simulations that essential conclusions of our von Neumann analysis remain valid for (in)stability of the MoC applied to that non-constant-coefficient system.

Physical model
While our study will focus on a linear problem of a rather general form (see Eq. (10a) below), we will begin by stating a specific nonlinear problem which had originally motivated our study and whose linearization leads to (10a). We consider the system where , and superscript 'T' denotes the transposition. This system is a representative of a class of models that arise in studying propagation of light in birefringent optical fibers with Kerr nonlinearity [4]- [7]; see Appendix A for more detail. The nonlinear system (5) has a rather special form and arises in a specific application. However, its linearized form (see Eq. (10a) below), for which we will analyze the stability of the MoC, is quite general. A wide and diverse range of physical problems lead to the same equation, as we will explain after Eqs. (11).
In the component form, system (5) is: It has four families of soliton/kink solutions, a special case of one of which is (see, e.g., [6]): (The other three families differ from (7) by combinations of signs of the solution's components.) However, we will focus on the analysis of the (in)stability of the MoC applied to a simpler, constant solution: This solution corresponds to the asymptotic value of (7) as x → −∞. It can be shown (see, e.g., [8]) that this solution is stable in the sense to be defined below. However, we have found that when simulated by certain "flavors" of the MoC, it can be numerically unstable. If that asymptotic, constant solution is numerically unstable, then so will be any "physically interesting" non-constant solutions, like (7), possessing the asymptotics (8). Thus, numerical stability of the constant solution (8) is necessary for the successful performance of the MoC on the soliton/kink solution (7) and related ones. Moreover, considering the numerical stability of the MoC simulating the simpler solution (8) will allow us to focus on the behavior of the numerical method without being distracted by the complexity of the physical model. The methodology of our analysis, as well as at least some of our general conclusions, will hold for the MoC applied to other non-dissipative hyperbolic PDEs; we will discuss this in Section 7 and provide details in the Appendices.
For future reference, we linearize Eqs. (6) on the background of solution (8): here S ± j 0 are the components of the exact solution (8) and s ± j are small perturbations. The linearized system (6) reduces to: where s = [s + 1 , s + 3 , s − 1 , s − 3 ] T , the 4 × 4 matrices in (10a) are: and the 2 × 2 matrices in (11a) are: Non-dissipative linearized equations of the form (10a), with the same or similar Σ but a more general form of P, arise in a wide range of physical applications involving the interaction of waves  [12]) and, more generally, parametric interaction among two or more electromagnetic waves [13] in semi-classical optics [14,15], quantum optics [16,17], photonics [18], and plasma physics applications [19,20]; interaction of crossed beams in photorefractive materials ( [9], Sec. 10.6); interaction of counter-and co-propagating light waves in a medium with a periodic refractive index ( [21,22,23]); and the relativistic field theory [24,25] (also see [26] for a general form of coupledwave-like field models with a cubic nonlinearity and [27] for more recent references on these models' solutions). Interestingly, models similar to those in the relativistic field theory also arise in nonlinear fiber optics [28,29]. In Section 7 we will demonstrate that our results for system (10), (11) are also applicable to the one-dimensional Gross-Neveu model [25] of the relativistic field theory.
In Eqs. (10) and (11) and in what follows we will adopt the following notations. Boldfaced quantities with an underline, S ± , and with a hat,Ĵ, will continue to denote 3 × 1 vectors and 3 × 3 matrices, respectively, as in (5). Boldfaced quantities without an underline or a hat will denote 4 × 4 matrices or 4 × 1 vectors, as in (10a); the ambiguity of the same notations for matrices and vectors here will not cause any confusion. Finally, underlined letters in regular (not boldfaced) font will denote 2 × 1 vectors; e.g.: Seeking the solution of (10a) to be proportional to e ikx−iωt , one can show that ω ∈ R for all k ∈ R. This means that solution (8) is stable on the infinite line, as mentioned two paragraphs above. In particular, ω(k = 0) ∈ R, which is equivalent to the statement that all nonzero eigenvalues of P are purely imaginary. Indeed, from (11) one finds: Returning to the full set of equations (5) and scalar-multiplying each of them by its respective S + or S − , we notice that they admit the following "conservation" relations: where | . . . | stands for the length of the vector. Two other conservation relations can be obtained in a similar fashion. Note that for periodic BC, considered in this paper, these relations become conservation laws. E.g., (13a) yield: where L is the length of the spatial domain. Moreover, a Hamiltonian of system (5) can be constructed in certain action-angle variables [30].
Therefore, it is desirable that the numerical method also conserve or quasi-conserve (i.e., not lead to a systematic shift) at least some of these quantities. A simple, explicit such a method for ODEs is LF, whereas explicit Euler methods are known to lead to mild numerical instability for conservative ODEs. For that reason we had expected that using the LF as the ODE solver in the MoC would produce better results than explicit Euler solvers. However, we found that, on the contrary, it produces the worst results. We will now turn to the description and analysis of the SE, ME, and LF solvers for the MoC with periodic BC. We will refer to the corresponding "flavors" of the MoC as MoC-SE, MoC-ME, and MoC-LF, respectively.

MoC with the SE solver
The SE is a first-order method and, moreover, is well-known to lead to a mild yet conspicuous numerical instability when applied to conservative ODEs (see, e.g., [31]). The reason that we consider this method is that it is simple enough to illustrate the approach. Thus, describing sufficient details in this Section will allow us to skip them in subsequent sections, devoted to second-order methods, where such details are more involved.

Analysis
The form of the MoC-SE equations for system (6) is: where f ± j are the nonlinear functions on the rhs of (6), and m = 0, . . . M , with (M + 1) being the number of grid points. Note that in this paper we have set the temporal and spatial steps equal, ∆t = ∆x = h, to ensure having a regular grid. To impose periodic BC, we make the following identifications in (14): Linearizing Eqs. (14), one arrives at: where O is the 2 × 2 zero matrix and Seeking the solution of (16) with periodic BC in the form (s) n m = (s) n e ikmh reduces (16) to: where z = kh, and Σ is defined in (11).
Below we will show that eigenvalues of N satisfy |λ N | > 1, which implies that the MoC-SE with periodic BC is unstable. One can easily see this for the longest-and shortest-period Fourier harmonics (z = 0 and z = π, respectively), where one has Q(z = 0) = I and Q(z = π) = −I and invokes relation (12). Then one has: However, for other Fourier harmonics of the numerical error, i.e., for z = 0 or π, an analytical relation between the eigenvalues λ N of the numerical method and the eigenvalues λ P of the "physical" matrix P cannot be established. In particular, other than by accident. Therefore, to determine |λ N |, one has to find these eigenvalues numerically.
Let us stress that this limits only the result, but not the methodology, of our von Neumann analysis to a specific matrix P. Indeed, for any P, the relation between the amplification matrix N of the numerical scheme and the "physical" matrix P is given by (18b). The eigenvalues of N are found within a second by modern software, whereas direct numerical simulations of scheme (14) (and other MoC schemes considered in Sections 4 and 5) take much longer.
The largest (in magnitude) eigenvalue of N, found by Matlab, is shown as a function of the normalized wavenumber z in Fig. 1(a). The two curves illustrate the general trend as h is varied.
They also show that the strongest instability of the MoC-SE applied to (5) occurs in the ODEand "anti-ODE" limits, i.e. for z = 0 and z = π. In particular, since Q(z = π) = −Q(z = 0), one concludes that |λ N (z = π)| = |λ N (z = 0)|. In other words, the instability of the MoC-SE for the highest and lowest Fourier harmonics is the same, in contrast to that of other numerical methods.

Numerical verification
First, in Fig. 1(b) we show the logarithm of the Fourier spectrum of the numerical error obtained by applying the MoC-SE to (6). Namely, we simulated scheme (14) with the initial condition being (8) plus white (in space) noise of magnitude 10 −12 and subtracted from the numerical solution the exact solution (8). The curves in Fig. 1(b) have the same qualitative shapes as the corresponding curves in Fig. 1(a), which confirms the validity of the analysis in Section 3.1. Indeed, it follows from (18a) that (individual) Fourier harmonics of the numerical error satisfy where . . . denotes the Euclidean norm of the vector.
Second, to verify our analytical results from yet another perspective, we will compare the growth rate measured for a "total" numerical error (see below) with the growth rate predicted by the analysis of Section 3.1. The "total" numerical error is computed as Note that the summation over m smooths out the noisy spatial profile of the error. The perturbations s ± 2 were indeed found to be much smaller than s ± 1,3 (as long as the latter are themselves sufficiently small), as predicted by (10b), and hence they are not included in the "total" error (21). Now, according to Section 3.1, the highest growth rate occurs for harmonics with z = 0, π. It follows from (19a) that Thus, the theoretical growth rate is where we have used (12).
On the other hand, the growth rate can be measured from the numerical solution: where t 1,2 are some times when the dependence of ln |s tot | on t appears to be linear (as, e.g., in Fig. 1(c) for t > 40). Using (24), we have found for the parameters listed in the caption to Fig. 1 that γ meas agrees with its theoretical value (23) to two significant figures.

MoC with the ME solver
The MoC-ME applied to system (6) yields the following scheme: where the notations are the same as in Section 3. Note that the rhs of (25a) is the same as that of (14). Following the analysis of Section 3.1, one obtains a counterpart of (18) where now As we discussed in Section 3.1, the eigenvalues of N cannot be analytically related to the eigenvalues of P and hence need to be found numerically. (As for the MoC-SE before, this takes fractions of a second for modern software, unlike the much longer direct simulations of scheme (25).) The magnitude of the largest-in-magnitude such eigenvalue is shown in Fig. 2(a). The logarithm of the Fourier spectrum of s is plotted in Fig. 2(b) and is seen to agree qualitatively with the result in Fig. 2(a); see also (20).  (25), (15) with that established by formula (27b). Note that the value 1 of the ratio would correspond to the perfect agreement between the numerics and analysis.
Let us note that the numerical instability of the MoC-ME scheme in the ODE and anti-ODE limits is much weaker than that of the MoC-SE, because from (26), , and hence similarly to (22), one has However, Fig. 2(a) shows that unlike the MoC-SE, the strongest instability of the MoC-ME occurs in the "middle" 1 of the Fourier spectrum. Moreover, this growth rate has the same order of magnitude as the growth rate of the MoC-SE's instability. Therefore, it is the growth rate of the harmonics with k ∼ π/(2h) that one would measure in the numerical experiment. Using the information from In Section 7 we will demonstrate that these conclusions about the numerical instability of the MoC-ME remain qualitatively valid for a different system of PDEs, which has spatially localized (as opposed to constant) coefficients; see Appendix B.
Unlike the Euler ODE solvers considered in the previous two sections, the LF involves the solution at three time levels: Given the initial condition, one can find the solution at time level n = 1 with, e.g., the MoC-SE.
To enforce periodic BC, convention (15) needs to be supplemented by The counterpart of (18) is now The (in)stability of the MoC-LF is determined by whether the magnitude of the largest eigenvalue of the quadratic eigenvalue problem obtained from (30) and leading to det λ 2 I − 2λhQ P + Q 2 = 0 (31) exceeds unity. As previously for the MoC-SE and MoC-ME, this 4 × 4 eigenvalue problem has to be solved numerically for a given P. The largest magnitude of its eigenvalue, computed with Matlab's command polyeig for P given by (11), is plotted in Fig. 3(a).  It should be noted that the quantity (|λ| − 1) for the MoC-LF scales as h, in contrast to the cases of MoC-SE and MoC-ME, where it scales as h 2 . Therefore, the instability growth rate of the MoC-LF is: γ = O(1), which is much greater than those rates of the MoC-SE and MoC-ME. The reason for this will come out as a by-product when we qualitatively explain, in the next paragraph, why the strongest instability of the MoC-LF for the PDE system in question occurs near z = π/2. First note that in the ODE limit, z = 0, one has Q(0) = I, and hence from (31) one recovers the well-known relation between the amplification factor λ of the ODE-LF and eigenvalues of P: Given λ P ∈ iR, as would be the case for any stable non-dissipative system, this yields |λ(z = 0)| = 1 (here and below we assume that h max |λ P | < 1). A similar situation holds in the anti-ODE limit, where Q = −I. However, for z = π/2, Q(π/2) = −i diag(I, −I) can be said to be the most dissimilar from I ≡ diag(I, I), and the counterpart of (32) is: where P mod = Q(π/2) P. For λ P mod ∈ iR, this yields max |λ(z = π/2)| > 1, i.e., a numerical instability. In the case of matrix P defined in (11), λ P mod = 0, 0, ±i √ 2, and thus the above statements explain our result shown in Fig. 3(a). In Section 7 we will discuss how this generalizes to some other energy-preserving hyperbolic PDEs.
Moreover, it also follows from (33) that the growth rate of this instability is about max |λ P mod | = O(1). (For the specific P in (11), we will show in the second paper of this work that the growth rate is 3/2 instead of √ 2 = max |λ P mod |; this minor difference is responsible for the subtly double-peaked profile of the curves in Fig. 3(a).) We confirmed this by direct numerical simulations: representative spectra of the numerical error are shown in Fig. 3(b) and are seen to agree with the shapes of the curves in Fig. 3(a); see (20). The solution obtained by the MoC-LF becomes completely destroyed

Comments about the MoC with the cRK solver
The fourth-order cRK method may appear as a natural choice of a higher-order ODE solver to be combined with the MoC. Indeed, not only is its accuracy for ODEs O(∆t 4 ), but also the systematic error that it introduces into the numerically computed conserved quantities of energy-preserving ODEs is usually negligible for sufficiently small time steps.
However, a straightforward use of the cRK solver in the MoC framework produces an unsatisfactory numerical method. The corresponding scheme is: A von Neumann analysis analogous to that presented in Section 4 yields a dependence max |λ N (z)| similar to the one shown in Fig. 2(a), and direct numerical simulations confirm that result for system (6). Thus, scheme (34) has the instability growth rate γ = O(h), just as the MoC-ME (25). In fact, for the specific PDE system (6), this γ is twice that of the MoC-ME's, and hence scheme (34) offers no advantage over the MoC-ME.
Perhaps, a proper MoC-cRK needs to ensure that the arguments of f ± in (34c) and (34d) be taken on the same characteristics. E.g., (34d) would then need to be replaced with which mimics the last term in (25). The node numbers of some of the arguments of f ± in (34c) would then also need to be modified in some way. This may result in a more stable MoC-cRK.
However, a detailed exploration of this topic will have to begin with an analysis of the accuracy of the resulting method away from the ODE limit. This is outside the scope of this paper, and therefore we will not pursue it here. Let us mention that we have been unable to find any published research which would study the problem of setting up the proper MoC-cRK for the case of crossing characteristics.

Summary and discussion
We have presented results of the von Neumann analysis of the MoC combined with three ODE solvers: SE, ME, and LF, as applied to an energy-preserving hyperbolic PDE system, (5) or (6).
Our main findings are as follows.
First, we have found that the numerical instability of the highest and lowest Fourier harmonics in the MoC schemes has the same growth rate. This is in contrast with the situation of other numerical methods, where it is the highest Fourier harmonics that become numerically unstable before the lower ones. Moreover, we have found that the harmonics in the "middle" of the spectrum, i.e., with |k| ≈ π/(2h), can become the most numerically unstable. This occurs, for the PDE system in question, for the MoC-ME and MoC-LF schemes. While such harmonics for the MoC-LF occupy a narrow spectral band and hence can, in principle, be filtered out, for the MoC-ME they occupy a substantial portion of the spectrum and hence cannot be filtered out without destroying the solution. It should also be pointed out that at least for some energy-preserving PDEs, as the one considered here, the growth rate of the most unstable harmonics has the same order of magnitude, Second, we have found that the MoC based on the LF solver, which (the solver) is known to considerably outperform the Euler schemes for energy-preserving ODEs, performs the worst for the energy-preserving PDE that we considered. Namely, the Fourier harmonics with |k| ≈ π/(2h) grow at the rate γ = O(1), which by far exceeds the growth rate γ = O(h) of the numerically most unstable harmonics in both the MoC-SE and MoC-ME.
In Section 6 we have pointed out that, for hyperbolic PDE systems with crossing characteristics, the formulation of a MoC based on the cRK ODE solver and having numerical stability (or, at worst, a milder instability than the MoC-ME) for |kh| ∼ π/2, remains an open problem.
There remains a question as to whether the above conclusions of our analysis, which necessarily had to be obtained for a specific system (as we explained in the main text), apply to other energypreserving hyperbolic PDEs. Again, for the same reasons, this cannot be answered in general, as no analytical relation between the eigenvalues of the "physical" matrix P and the numerical amplification factor |λ N | can be obtained for an arbitrary wavenumber k. Therefore, in Appendices A and B we consider several specific PDEs with the same linearized form (10a), where matrix P has a different form than (11). In Appendix A we limit ourselves to systems with spatially constant coefficients and focus on how qualitatively different eigenvalues (see below) of P affect the amplification factor for the MoC-LF. In Appendix B we will demonstrate that our conclusions about numerical instability of both MoC-LF and MoC-ME remain qualitatively valid for a different system whose linearization has spatially-localized entries in matrix P.
Thus, we will now explore what "patterns" of numerical instability can be expected in simulations of energy-preserving PDEs by the MoC-LF. To that end, we consider a number of systems 2 representing a broader class of such PDEs, which include our system (6), (8) as a special case.
The total of 18 systems, including it, were considered, and a summary of results is presented in Before we proceed, let us remind the reader that the stability of the physical problem (10a), as opposed to the stability of the numerical method, was defined before Eq. (12). For example, system (6), (8) is stable, which can be seen from Fig. 3(a). Indeed, the MoC-LF accurately approximates the solution for k = O(1) (i.e., for z ≡ kh = O(h) in Fig. 3(a)), and the fact that |λ| = 1 there corresponds to ω ∈ R in the text before Eq. (12). Let us stress that our focus is on the presence of numerical instability of the MoC-LF for physically stable systems. Indeed, even if the numerical method is found to be free of numerical instability for a physically unstable system, then the corresponding solution will eventually (typically within t = O(1)) evolve towards the physically stable solution.
We now continue with the summary of stability results for the MoC-LF. These were obtained by solving the quadratic eigenvalue problem (31) with Matlab's command polyeig for 18 specific matrices P corresponding to each of the systems listed in Appendix A. For all the seven of those systems which are physically stable, the dependence max |λ(z)| was found to be similar to that shown in Fig. 3(a). Thus, for all of the stable systems considered, the MoC-LF has numerically unstable modes with |kh| ≈ π/2 and growth rate γ = O(1). Thus, this method should be deemed unsuitable for simulation of these, and possibly other, stable energy-preserving PDEs with crossing characteristics.
Let us mention, for completeness, that for the physically unstable systems considered in Appendix A, we found two possibilities with respect to the numerical (in)stability of the MoC-LF.
In the first group, there are unstable systems such that: λ P = ±iα, ±iβ for α, β ∈ R (i.e., the mode with k = 0 is stable), with α = β and at most one of α or β could be zero. We have found that such systems have ω✚ ∈ R in the immediate vicinity of k = 0. For systems in this group, the MoC-LF has also numerically unstable modes with |k|h near (but not exactly at) π/2. This is illustrated in Fig. 4(a). In the second group there are unstable systems with all other possibilities of λ P compatible with the energy-preserving nature of the problem: λ P = ±iα (two pairs of repeated imaginary eigenvalues), including the case α = 0; and λ P = ±α, ±β (in this case, even the k = 0 mode is unstable). This is illustrated in Figs. 4(b,c), which show no numerical instability in the MoC-LF. However, as we have mentioned above, once the initial, physically unstable, solution evolves sufficiently near a stable one, the MoC-LF will be invalidated by the numerically unstable harmonics, which are seen in Figs. 3(a) and 4(a).
Let us now demonstrate the validity of qualitative conclusions of our von Neumann analysis for a system whose linearization (10a) has a spatially dependent matrix P. This system is the soliton (see Fig. 5(a)) of the Gross-Neveu model [25] in the relativistic field theory; its details are presented in Appendix B. It has received considerable attention in the past decade both from the analytical and numerical perspectives: see [27], [32]- [38], and references therein. The plots of the spectra of the numerical errors, obtained for this soliton with the MoC-LF (Fig. 5(b)) and MoC-ME ( Fig. 5(c)), are qualitatively similar to such plots in Sections 5 and 4, respectively. In particular, the result in panel (b) agrees with our conclusion that the MoC-LF is unsuitable for simulations of a physically stable hyperbolic PDE system with crossing characteristics: the numerical instability at k ≈ π/(2h) will destroy the solution soon after t = 20. In contrast, panel (c) shows that the MoC-ME would be a suitable method for such systems, even for relatively long times. Moreover, we also verified our conclusion (see (27b) and Fig. 2(c)) that the growth rate of the most unstable Fourier harmonics, with k ≈ π/(2h), of the MoC-ME scales as O(h) (which is considerably greater than the growth rate of that method's error for energy-preserving ODEs). To that end, we repeated the simulations for several values of h = L/M , where M = 2 11 , 2 12 , 2 13 , 2 14 and computed γ meas using (24) and (41). The results fit the dependence γ meas = 6.4 · 10 −3 · (2 11 /M ) · ln(10), which is equivalent to γ meas = O(h), to two significant figures. and that of an isotropic (i.e., non-birefringent) fiber [5]: Let us stress that only the above forms of matricesĴ c, s correspond to physically meaningful situations in the context of birefringent fibers, and thus it would not make sense to consider model (35) with arbitraryĴ c, s .
Models (35), (36) each have many stationary constant solutions. We will consider only six of them for each of those three models: for one of j = 1, 2, or 3, with the other two components of S ± being 0.
For brevity, we will refer to these solutions as (j±), where j and ± correspond to the particular choice of the component and the sign in (37). For example, solution (8) of the main text corresponds to (2−) in (37).
However, for Ω sufficiently away from 0, the soliton has been found to be stable by both analysis and numerics (see references in [37]). In Fig. 5(a) we show this solution for Ω = 0.7.
The linearized system (39) on the background of soliton (40) has the form (10a), where the entries of matrix P are localized functions of x. Therefore, a von Neumann analysis cannot be rigorously applied to it. However, one can expect that its predictions could be valid qualitatively, based on the principle of frozen coefficients. To demonstrate that this is indeed the case, we simulated system (39) with the initial condition shown in Fig. 5(a), to which a small white noise was added. The length of the computational domain was L = 64. The spectra of the numerical errors obtained by schemes (28) Note that the peaks near z = 0, seen in Figs. 5(b,c), do not correspond to any numerical instability.
In panel (b) that peak is caused directly by the discretization error of the scheme. In panel (c), a much higher such peak is also caused by the discretization error, but indirectly. Indeed, a slight error in the propagation constant Ω, which inevitably occurs due to a limited accuracy of the numerical scheme, causes {u n m , v n m } in (41) to differ from their respective {U (x m ), V (x m )} exp[−iΩt n ] by O(1) for t n ≫ 1.