Solution Interpolation Method for Highly Oscillating Hyperbolic Equations

Abstract

This paper deals with a novel numerical scheme for hyperbolic equations with rapidly changing terms. We are especially interested in the quasilinear equation 𝑢_𝑡+𝑎𝑢_𝑥=𝑓(𝑥)𝑢+𝑔(𝑥)𝑢^𝑛 and the wave equation 𝑢_𝑡𝑡=𝑓(𝑥)𝑢_𝑥𝑥 that have a highly oscillating term like 𝑓(𝑥)=sin(𝑥/𝜀),𝜀≪1. It also applies to the equations involving rapidly changing or even discontinuous coefficients. The method is based on the solution interpolation and the underlying idea is to establish a numerical scheme by interpolating numerical data with a parameterized solution of the equation. While the constructed numerical schemes retain the same stability condition, they carry both quantitatively and qualitatively better performances than the standard method.

1. Introduction

In the modern theory of partial differential equations (PDEs), the role of exact solutions is becoming more crucial, especially for a number of nonlinear equations or systems of PDEs. Since those problems do not allow many of the classical techniques in principle, exact solutions serve as windows revealing correct classes of existence, regularity, uniqueness, and specific asymptotics [1–4]. From the view point of numerical analysis, exact solutions are a basis for developing and testing numerical schemes and computer algebra software [5–10].

The research on solution interpolation started with the motivation to construct finite difference schemes, using a set of exact solutions to given differential equations [11]. The idea of solution interpolation is that once we find an exact solution that coincides with the given initial or boundary conditions on a grid, it can be used to make the numerical approximation which locally inherits all geometric properties of the original system, including oscillation embedded in the equation. This is not feasible over the whole domain in practical situations, but in the case that there exists a family of known solutions to the differential equation, we can still construct a finite difference scheme by interpolating them locally on a grid. That is, a numerical value at a specific point can be derived by interpolating the numerical data at around its neighborhood at the previous time step with some special solutions. Refer to Figure 1 for illustration of the idea.

Figure 1 

Once a specific solution that embeds three points 𝑝1, 𝑝2, and 𝑝3 is found, an approximated point 𝑝4 can be obtained from the solution.

In this paper, we first focus on the applications of the solution interpolation method to the quasilinear equation, 𝑢_𝑡+𝑎𝑢_𝑥=𝑓(𝑥)𝑢+𝑔(𝑥)𝑢^𝑘,𝑘∈𝐙−{1}. (1) The emphasis is placed on the cases where the equation involves highly oscillating or fast varying terms like 𝑓(𝑥)=sin(𝑥/𝜀), 𝜀≪1. The conventional numerical schemes have difficulty in dealing with such terms accurately.

We also develop a numerical scheme for the wave equation: 𝑢_𝑡𝑡=𝑐(𝑥)𝑢_𝑥𝑥,𝑐(𝑥)>0, (2) which describes wave propagation and reflection in inhomogeneous medium. If the local speed of propagation 𝑐(𝑥) slowly varies, asymptotic solutions can be found using the WKB method [12–14]. In the case of strongly inhomogeneous medium, 𝑐(𝑥) rapidly changes making both exact solution and numerical approximation hard to obtain.

In what follows, we will have a brief overview on the solution interpolation method. Sections 3 and 4 deal with applications of the method to various cases of (1), (2), and their numerical results. We will discuss error analysis of those schemes in the last section.

2. Solution Interpolation Method

For a given differential equation, suppose there is a family of exact solutions depending on 𝑛 parameters. If interpolation on 𝑛 points of a given stencil determines a unique function in this family, this naturally defines an (𝑛+1)-point numerical scheme. The scheme by solution interpolation shares geometrical properties with exact solutions such as Lie symmetry.

From here on we adopt common notations for numerical schemes; 𝑢^𝑛_𝑚 is an approximation of 𝑢(𝑥,𝑡) at a point (𝑥_𝑚,𝑡^𝑛) on a rectangular grid with the spacial and temporal grid size, ℎ=𝑥_𝑚+1−𝑥_𝑚, and 𝑘=𝑡^𝑛+1−𝑡^𝑛, respectively. A central difference operator 𝐷 is defined as 𝐷𝑢^𝑛_𝑚=𝑢^𝑛_𝑚+1−𝑢^𝑛_𝑚−1,𝐷2ℎ₂𝑢^𝑛_𝑚=𝑢^𝑛_𝑚+1−2𝑢^𝑛_𝑚+𝑢^𝑛_𝑚−1ℎ². (3) We also use left and right difference operators, 𝐷_𝑙𝑢^𝑛_𝑚=𝑢^𝑖_𝑗−𝑢^𝑛_𝑚−1ℎ,𝐷_𝑟𝑢^𝑛_𝑚=𝑢^𝑛_𝑚+1−𝑢^𝑖_𝑗ℎ. (4)

To illustrate how the idea of solution interpolation is implemented, let us take an example of a linear heat equation: 𝑢_𝑡=𝑎𝑢_𝑥𝑥,𝑎>0. (5) To construct an efficient four-point scheme by solution interpolation, it is natural to choose a set of solutions as simple as possible. In fact, the heat equation possesses a series of polynomial solutions, 𝜙₀=1,𝜙₁𝜙=𝑥,₂=𝑥²+2𝑎𝑡,𝜙₃=𝑥³𝜙+6𝑎𝑥𝑡,₄=𝑥⁴+12𝑎𝑥²𝑡+12𝑎²𝑡²,𝜙₅=⋯, (6) out of which we can construct a solution with three parameters as 𝑢(𝑥,𝑡)=𝐶₀𝜙₀+𝐶₁𝜙₁+𝐶₂𝜙₂=𝐶₀+𝐶₁𝑥+𝐶₂𝑥²+2𝑎𝑡. (7)

Interpolation of solution (7) by 𝑢^𝑛_𝑚+1, 𝑢^𝑛_𝑚, and 𝑢^𝑛_𝑚−1 at (𝑥_𝑚+1,𝑡^𝑛), (𝑥_𝑚,𝑡^𝑛), and (𝑥_𝑚−1,𝑡^𝑛), respectively, gives 𝐶₀=𝑢^𝑛_𝑚−𝑥_𝑚𝐷𝑢^𝑛_𝑚+𝑥²_𝑚2−𝑎𝑡^𝑛𝐷𝑢^𝑛_𝑚,𝐶₁=𝐷𝑢^𝑛_𝑚−𝑥_𝑚𝐷²𝑢^𝑛_𝑚,𝐶₂=12𝐷²𝑢^𝑛_𝑚. (8) Consequently, the numerical approximation for 𝑢_𝑚^𝑛+1 is determined as 𝑢_𝑚^𝑛+1=𝐶₀+𝐶₁𝑥_𝑚+𝐶₂𝑥_𝑚²+2𝑎𝑡^𝑛+1. (9) After simplification, the resulted scheme becomes 𝑢_𝑚^𝑛+1=𝑢^𝑛_𝑚+𝑎𝑘𝐷𝑢^𝑛_𝑚, (10) which is nothing but the traditional forward central scheme.

Now let us construct a six-point implicit scheme by the solution interpolation method. We need to interpolate five values of a function on a usual six-point stencil, using the first five polynomial solutions 𝜙₀,…,𝜙₄. After similar computation, the corresponding scheme turns out to be 1−6𝜇𝑢^𝑛+1_𝑚−1+10+12𝜇𝑢_𝑚^𝑛+1+1−6𝜇𝑢^𝑛+1_𝑚+1=1+6𝜇𝑢^𝑛_𝑚−1+10−12𝜇𝑢^𝑛_𝑚+1+6𝜇𝑢^𝑛_𝑚+1, (11) where 𝜇 is 𝑎(𝑘/ℎ²). This scheme is known as the Crandall method, which is unconditionally stable and has error 𝑂(ℎ⁴,𝑘²) [15].

We call a numerical scheme created from the above interpolation method a solution interpolation scheme (SI). For constant-coefficient differential equations like (5), this approach is no more than traditional discretization based the Taylor expansion. However, it is no more true if we deal with even little more complicated equations. One can refer to [11] to see many cases where simple solutions lead to interesting nonconventional numerical schemes.

3. Solution Interpolation for Quasilinear Equations

Equation (1) is known to have a general solution in the form of 𝑢^1−𝑘(𝑥,𝑡)=𝐻(𝑥)𝜙(𝑥−𝑎𝑡)+1−𝑘𝑎𝐻(𝑥)𝑔(𝑥)𝐻(𝑥)𝑑𝑥, (12) where 𝜙(𝑥) is an arbitrary differentiable function and 𝐻(𝑥)=exp(((1−𝑘)/𝑎)∫𝑓(𝑥)𝑑𝑥). Note that we can easily derive a family of special solutions of (1) by setting 𝜙 as a certain parameterized function.

In this section, we apply the solution interpolation method to several cases of convection equation (1) as nonhomogeneousequation:𝑢_𝑡+𝑎𝑢_𝑥=𝑓(𝑥),varying-coefficientequation:𝑢_𝑡+𝑎𝑢_𝑥(𝑥)𝑢=𝑓𝑢,quasilinearvarying-coefficientequation:_𝑡+𝑎𝑢_𝑥=𝑓(𝑥)𝑢², (13) to construct the corresponding SIs. The numerical performances of the schemes will be compared to those of conventional schemes, the upwind scheme, and the Lax-Wendroff scheme.

The solution interpolation scheme often seemingly involves a function 𝐹(𝑥), an antiderivative of 𝑓(𝑥) that is in a coefficient or a nonhomogeneous term of the original differential equation. But it is important to note that the schemes do not necessarily require 𝐹(𝑥) in an explicit form. Since 𝐹(𝑥) only appears with the difference operator 𝐷 in the schemes, one can manage with numerical integration of 𝑓(𝑥) instead. For example, 𝐷𝐹(𝑥_𝑚) and 𝐷²𝐹(𝑥_𝑚) can be approximated as 𝐷𝐹𝑥_𝑚≈1𝑓𝑥12^𝑚+1+4𝑓𝑥_𝑚+1/2+2𝑓𝑥_𝑚+4𝑓𝑥_𝑚−1/2+𝑓𝑥_𝑚−1,𝐷²𝐹𝑥_𝑚≈1𝑓𝑥6ℎ^𝑚+1+4𝑓𝑥_𝑚+1/2−4𝑓𝑥_𝑚−1/2−𝑓𝑥_𝑚−1, (14) if Simpson’s rule is adopted. We call the scheme using substitution (14) an approximate solution interpolation scheme (ASI), to tell from SI using 𝐹(𝑥) explicitly.

Throughout the examples, we assume that the coefficient 𝑎 is fixed at 1 and the function 𝑓 is given as 𝑓(𝑥)𝑥=sin𝜀, (15) with 𝜀 varying between 0 and 1. We want to see how oscillation in the equations affects the performance of the scheme. In all examples, the spacial and temporal step sizes are fixed at ℎ=0.01 and 𝑘=0.005, respectively.

Example 1. Nonhomogeneous linear convection equation: 𝑢_𝑡+𝑎𝑢_𝑥=𝑓(𝑥). (16) This equation has a general solution, ()+1𝑢=𝜙𝑥−𝑎𝑡𝑎𝐹(𝑥), (17) where 𝜙 is an arbitrary differentiable function. Choosing 𝜙(𝑥)=𝐶₀+𝐶₁𝑥+𝐶₂𝑥², we solve the interpolation equations, 𝑢^𝑛_𝑙=𝜙𝑥_𝑙−𝑎𝑡^𝑛+1𝑎𝐹𝑥_𝑙,𝑙=𝑚−1,𝑚,𝑚+1. (18) Considering that (16) does not involve 𝑡, we can set 𝑡^𝑛=0 in (18) without losing generality. Then the solution is 𝐶₀=1−𝑥_𝑚𝑥𝐷+^𝑚22𝐷²𝑢^𝑛_𝑚−1𝑎𝐹𝑥_𝑚,𝐶₁=𝑥𝐷−^𝑚2𝐷²𝑢^𝑛_𝑚−1𝑎𝐹𝑥_𝑚,𝐶₂=𝐷²𝑢^𝑛_𝑚−1𝑎𝐹𝑥_𝑚. (19) Plugging this result into 𝑢_𝑚^𝑛+1=𝜙𝑥_𝑚+1−𝑎𝑘𝑎𝐹𝑥_𝑚 (20) gives the solution interpolation scheme, 𝑢_𝑚^𝑛+1=𝑢^𝑛_𝑚+𝑎−𝑎𝑘𝐷+²𝑘²2𝐷²𝑢^𝑛_𝑚+1𝑎𝐹𝑥_𝑚. (21) The scheme might be similar in appearance to the Lax-Wendroff scheme, 𝑢_𝑚^𝑛+1=𝑢^𝑛_𝑚+𝑎−𝑎𝑘𝐷+²𝑘²2𝐷²𝑢^𝑛_𝑚+𝑘−𝑎𝑘²2𝐷𝑓𝑥_𝑚. (22) However, there is significant difference between two schemes in numerical performance, which is illustrated in Figure 2. One can easily confirm that SI and ASI excel the upwind and the Lax-Wendroff schemes. More interestingly in Figure 2(b), SI and ASI maintain almost the same accuracy even with sever oscillation (𝜀=0.01), while two others fail.

(a) 𝜀=1

(b) 𝜀=0.01

Figure 2 

Numerical solutions for 𝑢𝑡+𝑎𝑢𝑥=sin(𝑥/𝜀): a solution 𝑢=(𝑥−𝑡)(sin(𝑥−𝑡)+2)−1−𝜀cos(𝑥/𝜀) is used to generate initial/boundary conditions and measure errors.

Example 2. Linear convection equation with a varying coefficient: 𝑢_𝑡+𝑎𝑢_𝑥=𝑓(𝑥)𝑢. (23) A general solution for (23) is 1𝑢=exp𝑎𝐹(𝑥)𝜙(𝑥−𝑎𝑡). (24) The choice of the same function 𝜙(𝑥)=𝐶₀+𝐶₁𝑥+𝐶₂𝑥² as in the previous example eventually leads to the scheme 𝑢_𝑚^𝑛+1=𝛼𝑢^𝑛_𝑚+1+𝛽𝑢^𝑛_𝑚+𝛾𝑢^𝑛_𝑚−1, (25) where the coefficients 𝛼, 𝛽, and 𝛾 are determined as 𝛼=𝑎𝜆(−1+𝑎𝜆)exp−ℎ𝐷_𝑟𝐹𝑥_𝑚,𝛽=1−𝑎²𝜆²,𝛾=𝑎𝜆(1+𝑎𝜆)expℎ𝐷_𝑙𝐹𝑥_𝑚. (26) Here 𝜆 denotes a Courant number 𝑘/ℎ. The performance of the scheme is shown in Figure 3. One can see that SI and ASI are just as good as the Lax-Wendroff scheme when 𝜀=1. But in contrast to the Lax-Wendroff that suffers from oscillation term in (b), two schemes keep their accuracy level regardless of oscillation intensity.

(a) 𝜀=1

(b) 𝜀=0.01

Figure 3 

Numerical solutions for 𝑢𝑡+𝑎𝑢𝑥=sin(𝑥/𝜀)𝑢: a solution 𝑢=exp(−𝜀cos(𝑥/𝜀)−𝑥+𝑡)(𝑥−𝑡)2+1 is used to generate initial/boundary conditions and measure errors.

Example 3. Quasilinear convection equation with a varying coefficient: 𝑢_𝑡+𝑎𝑢_𝑥=𝑓(𝑥)𝑢². (27) The equation has a general solution, 1𝑢=𝜙(𝑥−𝑎𝑡)−𝑎𝐹(𝑥)⁻¹. (28) We take the same parameterized function 𝜙 as in two previous examples and apply the solution interpolation. After some algebraic work, the scheme is derived as 𝑢_𝑚^𝑛+1=𝑢^𝑛_𝑚+1𝑢^𝑛_𝑚𝑢^𝑛_𝑚−1×𝑢^𝑛_𝑚+1𝑢^𝑛_𝑚−1+𝑎𝑎𝑘𝑆+²𝑘²2𝑆²𝑢^𝑛_𝑚+𝑢^𝑛_𝑚+1𝑢^𝑛_𝑚𝑢^𝑛_𝑚−1×𝑎−𝑎𝑘𝐷+²𝑘²2𝐷²𝐹(𝑥_𝑚)⁻¹. (29) Here, we use an additional difference operator 𝑆 as 𝑆𝑢^𝑛_𝑚=𝑢^𝑛_𝑚+1𝑢^𝑛_𝑚−𝑢^𝑛_𝑚𝑢^𝑛_𝑚−1,𝑆2ℎ₂𝑢^𝑛_𝑚=𝑢^𝑛_𝑚+1𝑢^𝑛_𝑚−2𝑢^𝑛_𝑚𝑢^𝑛_𝑚−1+𝑢^𝑛_𝑚𝑢^𝑛_𝑚−1ℎ², (30) for a simplified expression. Once again, it is shown in Figure 4 that SI and ASI overwhelm two others, especially when there is severe spatial oscillation generated by the coefficient sin(𝑥/𝜀), 𝜀=0.01.

(a) 𝜀=1

(b) 𝜀=0.01

Figure 4 

Numerical solutions for 𝑢𝑡+𝑎𝑢𝑥=sin(𝑥/𝜀)𝑢2: a solution 𝑢=(𝑥2−2𝑥𝑡+𝑡2+1)(1+𝜀(𝑥2−2𝑥𝑡+𝑡2+1)cos(𝑥/𝜀))−1 is used to generate initial/boundary conditions and measure errors.

4. Solution Interpolation for Wave Equations

In this section, we apply the solution interpolation method to the wave equation (2) with variable coefficient 𝑐(𝑥). It is known that (2) has two families of particular solutions as follows [16]. Solutions with even powers of 𝑡 are 𝑢_𝑙(𝑥,𝑡)=^𝑙_𝑘=0𝜙_𝑘(𝑥)𝑡^{2(𝑙−𝑘)},𝑙=0,1,2,…, (31) where the function 𝜙=𝜙_𝑘(𝑥) is defined by the recurrence relations 𝜙₀(𝑥)=𝐴₀𝑥+𝐵₀,𝜙_𝑘+1(𝑥)=𝐴_𝑘𝑥+𝐵_𝑘()×+2𝑙−𝑘)(2𝑙−2𝑘−1_𝑥𝑎(𝑥−𝑠)𝜙_𝑘(𝑠)𝑐(𝑠)𝑑𝑠,𝑘=1,2,…,𝑙−1. (32) Here 𝐴_𝑘, 𝐵_𝑘, and 𝑎 are arbitrary constants. Similarly, there is another family of particular solutions with odd powers of 𝑡 as 𝑢_𝑙(𝑥,𝑡)=^𝑙_𝑘=0𝜑_𝑘(𝑥)𝑡^{2𝑙−2𝑘+1}, (33) where the function 𝜑=𝜑_𝑘(𝑥) is defined by the recurrence relations: 𝜑₀(𝑥)=𝐴₀𝑥+𝐵₀,𝜑_𝑘+1(𝑥)=𝐴_𝑘𝑥+𝐵_𝑘()×+2𝑙−𝑘)(2𝑙−2𝑘+1_𝑥𝑎(𝑥−𝑠)𝜑_𝑛(𝑠)𝑐(𝑠)𝑑𝑠,𝑘=1,2,…,𝑙−1. (34) Again, 𝐴_𝑘, 𝐵_𝑘, and 𝑎 are arbitrary constants.

Among the above solutions, let us choose a four-parameter solution, 𝑢(𝑥,𝑡)=𝐶₀+𝐶₁𝑥+𝐶₂𝑡+𝐶₃𝑡²+2𝐶₃^𝑥₀(𝑥−𝑠)𝑐(𝑠)𝑑𝑠, (35) out of which we construct a five-point scheme for the equation. Since (2) does not explicitly depend on time 𝑡, we can assume 𝑡^𝑛=0 without loss of generality. Interpolating the numerical data 𝑢^𝑛_𝑚−1, 𝑢^𝑛_𝑚, 𝑢^𝑛_𝑚+1, and 𝑢_𝑚^𝑛−1 on the grid with solution (35) determines the parameters as 𝐶₀=𝐹𝑥_𝑚𝐷²𝑢^𝑛_𝑚+𝑥_𝑚𝐷²𝐹𝑥_𝑚𝐷𝑢^𝑛_𝑚−𝑢^𝑛_𝑚𝐷²𝐹𝑥_𝑚−𝑥_𝑚𝐷²𝑢^𝑛_𝑚𝐷𝐹𝑥_𝑚×𝐷²𝐹(𝑥_𝑚)⁻¹,𝐶₁=𝐷²𝑢^𝑛_𝑚𝐷𝐹𝑥_𝑚−𝐷²𝐹𝑥_𝑚𝐷𝑢^𝑛_𝑚𝐷²𝐹𝑥_𝑚,𝐶₂=𝑘²𝐷²𝑢^𝑛_𝑚+2𝑢^𝑛_𝑚−𝑢_𝑚^𝑛−1𝐷²𝐹𝑥_𝑚2𝑘𝐷²𝐹𝑥_𝑚,𝐶₃=𝐷²𝑢^𝑛_𝑚2𝐷²𝐹𝑥_𝑚, (36) where 𝐹(𝑥)=∫^𝑥₀((𝑥−𝑠)/𝑐(𝑠))𝑑𝑠. We plug this result to 𝑢_𝑚^𝑛+1=𝑢(𝑥_𝑚,𝑡^𝑛) in (35) and then eventually obtain the scheme: 𝑢_𝑚^𝑛+1=2𝑢^𝑛_𝑚−𝑢_𝑚^𝑛−1+𝑘²𝐷²𝑢^𝑛_𝑚𝐷²𝐹𝑥_𝑚, (37) after simplification. Note that we do not have to integrate (𝑥−𝑠)/𝑐(𝑠) and find 𝐹(𝑥) explicitly. Instead, by using numerical integration method, we can replace the factor involving 𝐹(𝑥_𝑚) with values of 𝑐(𝑥) on some points. For example, adopting Simpson’s rule gives 𝐷²𝐹𝑥_𝑚1≈3𝑐𝑥_𝑚−1/2+1𝑐𝑥_𝑚+1𝑐𝑥_𝑚+1/2. (38)

One can easily see from (37) similarity and difference between the approximate solution interpolation scheme and the standard central forward scheme (CF). Both schemes share the same form: 𝑢_𝑚^𝑛+1=2𝑢^𝑛_𝑚−𝑢_𝑚^𝑛−1+𝑘²𝛼_𝑚𝐷²𝑢^𝑛_𝑚, (39) where 𝛼_𝑚=𝑐𝑥_𝑚𝛼∶CF,_𝑚=3𝑐𝑥_𝑚−1/2𝑐𝑥_𝑚𝑐𝑥_𝑚+1/2×𝑐𝑥_𝑚−1/2𝑐𝑥_𝑚+𝑐𝑥_𝑚−1/2𝑐𝑥_𝑚+1/2+𝑐𝑥_𝑚𝑐𝑥_𝑚+1/2⁻¹∶ASI. (40)

Now let us show the performances of ASI in comparison to CF, focusing especially on the cases where the propagation speed 𝑐(𝑥) changes rapidly or discontinuously.

Example 4. Wave equation with exponentially decreasing speed: 𝑢_𝑡𝑡−𝑥=exp𝜀𝑢_𝑥𝑥. (41) We apply two schemes ASI and CF to this equation with two different values of 𝜀. Figure 5(a) shows that there are no qualitative differences between two schemes for 𝜀=1. However, when 𝜀 is 0.05 and the local propagation speed rapidly diminishes along the 𝑥-axis, one can see that ASI yields much better result while CF loses its accuracy quickly as time goes.

(a) 𝜀=1

(b) 𝜀=0.05

Figure 5 

Time evolution of the errors in numerical solutions for 𝑢𝑡𝑡=exp(−𝑥/𝜀)𝑢𝑥𝑥,0≤𝑥≤1. The exact solution 𝑢=(1+𝑥)𝑡2 + 𝑥𝑡+𝑥+2(−𝜀2+2𝜀3+(−𝜀+𝜀2)𝑥+(𝜀2−2𝜀3+𝜀2𝑥)exp(𝑥/𝜀)) is used to generate the initial/boundary conditions and measure the errors.

Example 5. Wave equation with oscillating speed: 𝑢_𝑡𝑡=1()𝑢2+sin−𝑥/𝜀^𝑥𝑥. (42) This equation involves oscillation in the wave speed because of the sine function in the denominator. Once again, if 𝜀=1 and therefore the speed changes slowly, two schemes show little qualitative differences. The outstanding difference occurs when 𝜀 drops to 0.01 and the propagation speed undergoes severe oscillation. It is obvious from Figure 6(b) that CF suffers from rapid change, but ASI successfully captures the solution in numerical simulation. Even more interestingly, the error of ASI with oscillation in (b) slightly decreases, compared to that in (a).

(a) 𝜀=1

(b) 𝜀=0.01

Figure 6 

Time evolution of the errors in numerical solutions for 𝑢𝑡𝑡=(2+sin(−𝑥/𝜀))−1𝑢𝑥𝑥,0≤𝑥≤1. The exact solution 𝑢=(1+𝑥)𝑡2+𝑥𝑡+𝑥+2(2𝜀3+𝜀𝑥+𝑥2+𝑥3/3−2𝜀3cos(𝑥/𝜀)−𝜀2(1+𝑥)sin(𝑥/𝜀)) is used to generate the initial/boundary conditions and measure the errors.

Example 6. Wave equation with discontinuous speed: 𝑢_𝑡𝑡=𝑐(𝑥)𝑢_𝑥𝑥, (43) where 𝑐(𝑥) is a step function as 𝑐(𝑥)=1if𝑥<0,2if𝑥≥0. (44) Equation (43) describes the propagation of waves across two different types of medium and its solution may be discontinuous. Applications of ASI and CF to the equation with two different initial/boundary conditions are illustrated in Figure 7. In both cases, CF fails to reproduce this discontinuity and loses sharp transition layer, even generating spurious oscillation. On the contrary, ASI well captures the qualitative aspect of the exact solutions, handling discontinuity successfully.

(a)

(b)

Figure 7 

Snapshots of numerical solutions for (43) at 𝑡=0.3. In (a), the exact solution is 𝑢=exp(𝑥−𝑡)−exp(−𝑥+𝑡)/2 if 𝑥<0 and 𝑢=exp(−𝑥/√2−𝑡) − exp(𝑥/√2+𝑡)/2 otherwise. In (b), the exact solution is 𝑢=exp(𝑥+𝑡) − exp(−𝑥−𝑡) if 𝑥<0 and 𝑢=exp(𝑥/√2−𝑡)+exp(−𝑥/√2+𝑡) otherwise.

5. Analysis of Solution Interpolation Schemes

The examples from the last two sections confirm that the schemes based on the solution interpolation methods achieve generally better, or compatible performances at the least, in comparison to the conventional schemes. Consistency of those schemes is guaranteed by the fact that they are created from solutions. In fact, Taylor expansion of schemes (21), (25), and (29) shows that they are of the second order, just as the Lax-Wendroff scheme. Likewise, solution interpolation scheme (37) has the same order of error as the central forward method.

Now, let us investigate the stability conditions for the solution interpolation schemes. In order to handle a nonhomogeneous term, it is helpful to modify usual notations for the von Neumann stability analysis. We replace 𝑢^𝑛_𝑚 and 𝑢_𝑚^𝑛+1 in a scheme with 𝑔(𝑡^𝑛)𝑒^𝑖𝑚𝜃 and 𝑔(𝑡^𝑛+𝑘)𝑒^𝑖𝑚𝜃, respectively. The scheme is stable as long as the amplification factor 𝑔(𝑡) satisfies 𝑔𝑡^𝑛+𝑘≤𝑔𝑡^𝑛for−𝜋≤𝜃≤𝜋,𝑘≥0. (45) In fact, this stability condition can be relaxed as 𝑔𝑡^𝑛+𝑘≤𝑔𝑡^𝑛+𝑂(𝑘)for−𝜋≤𝜃≤𝜋,𝑘≥0. (46) We use the above substitution and rewrite scheme (21) as 𝑔𝑡^𝑛+𝑘𝑒^𝑖𝑚𝜃=𝑔𝑡^𝑛𝑒^𝑖𝑚𝜃−𝑎𝜆2𝑔𝑡^𝑛𝑒^{𝑖(𝑚+1)𝜃}−𝑔𝑡^𝑛𝑒^{𝑖(𝑚−1)𝜃}+𝑎²𝜆²2𝑔𝑡^𝑛𝑒^{𝑖(𝑚+1)𝜃}−2𝑔𝑡^𝑛𝑒^𝑖𝑚𝜃+𝑔𝑡^𝑛𝑒^{𝑖(𝑚−1)𝜃}+𝑎−𝑎𝑘𝐷+²𝑘²2𝐷²𝐹𝑥_𝑚, (47) where 𝜆=𝑘/ℎ. This can be simplified and reorganized as 𝑔𝑡^𝑛+𝑘=𝛼𝑔𝑡^𝑛+𝛽, (48) where 𝛼=1−𝑎²𝜆²(1−cos𝜃)−𝑖𝑎𝜆sin𝜃,𝛽=𝑒^{−𝑖𝑚𝜃}𝑎−𝑎𝑘𝐷+²𝑘²2𝐷²𝐹𝑥_𝑚. (49) If we assume |𝛼|≤1, then 𝑔𝑡^𝑛+𝑘=𝛼𝑔𝑡^𝑛≤|𝛼𝑔𝑡+𝛽|^𝑛+𝛽≤𝑔𝑡^𝑛+𝛽, (50) and since 𝛽=𝑂(𝑘), the scheme would be stable. By the way, |𝛼|²=1−2𝑎²𝜆²sin²𝜃2²+(𝑎𝜆sin𝜃)²=1−4𝑎²𝜆²1−𝑎²𝜆²sin⁴𝜃2. (51) Therefore we see |𝛼|≤1⇔(1−𝑎²𝜆²)≥0, and this implies that the stability condition for scheme (21) is |𝑎𝜆|≤1. (52) In the same way, we can confirm that the Lax-Wendroff scheme (22) shares the same condition for stability.

Scheme (29) also has the identical stability region because it can be viewed as a numerical scheme for the same equation 𝑣_𝑡+𝑎𝑣_𝑥=𝑓(𝑥) after substitution 𝑣=1/𝑢. Similarly, scheme (25) can be converted to another solution interpolation scheme for the homogeneous advection equation: 𝑣_𝑡+𝑎𝑣_𝑥=0, (53) by substitution of 𝑣=exp(−1/𝑎𝐹(𝑥))𝑢. Note that this substitution does not affect the stability region. The solution interpolation scheme for (53) is 𝑢_𝑚^𝑛+1=𝑢^𝑛_𝑚+𝑎−𝑎𝑘𝐷+²𝑘²2𝐷²𝑢^𝑛_𝑚, (54) which is nothing but the Lax-Wendroff scheme and is stable if |𝑎𝜆|≤1. Therefore, scheme (25) has the same numerical stability too.

For the variable coefficient wave equation (2), the general procedure is to regard 𝑐(𝑥) as a frozen coefficient for each 𝑥, 𝑡 values in question. If the scheme is stable for each frozen coefficient, then the scheme is guaranteed to be stable. The corresponding CFL condition for the central forward scheme is |𝑐(𝑥1)|𝜆≤1or𝜆≤|𝑐(𝑥.)| (55) In the case of the solution interpolation scheme (37), one can show that the condition becomes 1𝜆≤|𝑐(𝑥ℎ)|+𝑂², (56) which is almost the same as the original stability condition. This is derived by expanding the coefficient functions in (40) in a Taylor series. The condition is based on assumption that 𝑐(𝑥) is differentiable. However, the analysis is valid regardless of its differentiability since the scheme approximates 𝑐(𝑥) with polynomials locally.

Since the schemes based on solution interpolation share the same consistency and stability with the conventional schemes, the reason why those are more robust against oscillation or rapid change embedded in the equation should be sought elsewhere. One of the possible explanations is that we utilize analytic information on the equations in the scheme construction. Solution interpolation makes numerical schemes reflect the additional structures of differential equations that are likely lost in usual discretization techniques.

Let us discuss characteristic lines of the equations, as an example of the structure which is successfully inherited to numerical schemes through the solution interpolation method. The following differential equations of increasing order can be derived one by one from (16) as 𝑢_𝑡𝑡−𝑎²𝑢_𝑥𝑥=−𝑎𝑓^′(𝑥),𝑢_𝑡𝑡𝑡+𝑎³𝑢_𝑥𝑥𝑥=𝑎²𝑓^′′(𝑥),𝑢_{𝑡𝑡𝑡𝑡}−𝑎⁴𝑢_{𝑥𝑥𝑥𝑥}=−𝑎³𝑓^′′′(𝑥),𝑢_{𝑡𝑡𝑡𝑡𝑡}+𝑎⁵𝑢_{𝑥𝑥𝑥𝑥𝑥}=−𝑎⁴𝑓^′′′′(𝑥)⋮ (57) Suppose 𝑢(𝑥_𝑚,𝑡^𝑛+1) is the exact value of a solution and 𝑢_𝑚^𝑛+1 is the value obtained from scheme (21). Application of the Taylor expansion to 𝑢(𝑥_𝑚,𝑡^𝑛+1)−𝑢_𝑚^𝑛+1 at (𝑥_𝑚,𝑡^𝑛) gives 𝑢𝑥_𝑚,𝑡^𝑛+1−𝑢_𝑚^𝑛+1=𝑢_𝑡+𝑎𝑢_𝑥(𝑥)𝑘+1−𝑓2𝑢_𝑡𝑡−𝑎²𝑢_𝑥𝑥+𝑎𝑓^′(𝑥)𝑘²+16𝑢_𝑡𝑡𝑡+𝑎𝜆²𝑢_𝑥𝑥𝑥−1𝜆²𝑓^′′(𝑥)𝑘³+1𝑢24^{𝑡𝑡𝑡𝑡}−𝑎²𝜆²𝑢_{𝑥𝑥𝑥𝑥}+𝑎𝜆²𝑓^′′′(𝑥)𝑘⁴+1𝑢120^{𝑡𝑡𝑡𝑡𝑡}+𝑎𝜆⁴𝑢_{𝑥𝑥𝑥𝑥𝑥}−1𝜆⁴𝑓^′′′′(𝑥)𝑘⁵+⋯. (58) Note that the first two terms on the right hand side are always zero. This indicates that the error of the scheme is 𝑂(𝑘³,ℎ²𝑘). Now, also suppose 𝜆=1/𝑎 and compare the differential equations in the list in (57) and the terms in the series (58). We can see that the third, fourth, and fifth terms become zero as well. In fact, a little bit of algebraic work shows that all terms in the series vanish simultaneously with 𝜆=1/𝑎, which implies that 𝑢_𝑚^𝑛+1 which we found is exact.

It is no wonder this happens since 𝑢_𝑚^𝑛+1 is on the same characteristic line with 𝑢^𝑛_𝑚 when 𝜆=1/𝑎. In this case, change of 𝑢 along this line can be traced by solving a simple ordinary differential equation and especially at (𝑥_𝑚,𝑡^𝑛+1), we obtain 𝑢𝑥_𝑚,𝑡^𝑛+1=𝑢𝑥_𝑚−1,𝑡^𝑛+^𝑥𝑚_𝑥𝑚−1𝑓(𝑥)𝑑𝑥. (59) This is in fact exactly what we get from scheme (21) with ℎ=𝑎𝑘, 𝑢_𝑚^𝑛+1=𝑢^𝑛_𝑚−1+𝐷_𝑙𝐹𝑥_𝑚. (60) The similar situation occurs with schemes (25) and (29). Note that this observation holds as long as a characteristic line connects the nodes of the grid. However, for general 𝜆, we can view the solution interpolation method as a natural way to extend characteristic lines and find an approximate value at the target point.

Acknowledgments

The corresponding author (P. Kim) acknowledges that this work was supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education, Science and Technology (2011-0023486). C.H. Lee's work was also supported by NRF funded by the Ministry of Education, Science and Technology (2010-0024849).