# Solution Interpolation Method for Highly Oscillating Hyperbolic Equations

Authors: Pilwon Kim; Chang Hyeong Lee
Journal: Journal of Applied Mathematics
 (2013)

---

## Abstract

This paper deals with a novel numerical scheme for hyperbolic equations with
rapidly changing terms. We are especially interested in the quasilinear equationut+aux=f(x)u+g(x)un and the wave equation utt=f(x)uxx that have a highly
oscillating term like f(x)=sin(x/ε),ε≪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.


---

## Body

## 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 pointsp1, p2, and p3 is found, an approximated point p4 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)ut+aux=f(x)u+g(x)uk,k∈Z-{1}.
The emphasis is placed on the cases where the equation involves highly oscillating or fast varying terms like f(x)=sin(x/ε), ε≪1. The conventional numerical schemes have difficulty in dealing with such terms accurately.We also develop a numerical scheme for the wave equation:(2)utt=c(x)uxx,c(x)>0,
which describes wave propagation and reflection in inhomogeneous medium. If the local speed of propagation c(x) slowly varies, asymptotic solutions can be found using the WKB method [12–14]. In the case of strongly inhomogeneous medium, c(x) 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. Sections3 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 onn parameters. If interpolation on n points of a given stencil determines a unique function in this family, this naturally defines an (n+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;umn is an approximation of u(x,t) at a point (xm,tn) on a rectangular grid with the spacial and temporal grid size, h=xm+1-xm, and k=tn+1-tn, respectively. A central difference operator D is defined as
(3)Dumn=um+1n-um-1n2h,D2umn=um+1n-2umn+um-1nh2.
We also use left and right difference operators,
(4)Dlumn=uji-um-1nh,Drumn=um+1n-ujih.To illustrate how the idea of solution interpolation is implemented, let us take an example of a linear heat equation:(5)ut=auxx,a>0.
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,
(6)ϕ0=1,ϕ1=x,ϕ2=x2+2at,ϕ3=x3+6axt,ϕ4=x4+12ax2t+12a2t2,ϕ5=⋯,
out of which we can construct a solution with three parameters as
(7)u(x,t)=C0ϕ0+C1ϕ1+C2ϕ2=C0+C1x+C2(x2+2at).Interpolation of solution (7) by um+1n, umn, and um-1n at (xm+1,tn), (xm,tn), and (xm-1,tn), respectively, gives
(8)C0=umn-xmDumn+(xm22-atn)Dumn,C1=Dumn-xmD2umn,C2=12D2umn.
Consequently, the numerical approximation for umn+1 is determined as
(9)umn+1=C0+C1xm+C2((xm)2+2atn+1).
After simplification, the resulted scheme becomes
(10)umn+1=umn+akDumn,
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ϕ0,…,ϕ4. After similar computation, the corresponding scheme turns out to be
(11)(1-6μ)um-1n+1+(10+12μ)umn+1+(1-6μ)um+1n+1=(1+6μ)um-1n+(10-12μ)umn+(1+6μ)um+1n,
where μ is a(k/h2). This scheme is known as the Crandall method, which is unconditionally stable and has error O(h4,k2) [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
(12)u1-k(x,t)=H(x)ϕ(x-at)+(1-ka)H(x)∫‍g(x)H(x)dx,
where ϕ(x) is an arbitrary differentiable function and H(x)=exp(((1-k)/a)∫f(x)dx‍). 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
(13)nonhomogeneousequation:ut+aux=f(x),varying-coefficientequation:ut+aux=f(x)u,quasilinearvarying-coefficientequation:ddddddddddut+aux=f(x)u2,
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 functionF(x), an antiderivative of f(x) 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 F(x) in an explicit form. Since F(x) only appears with the difference operator D in the schemes, one can manage with numerical integration of f(x) instead. For example, DF(xm) and D2F(xm) can be approximated as
(14)DF[xm]≈112(f(xm+1)+4f(xm+1/2)ddd+2f(xm)+4f(xm-1/2)+f(xm-1)),D2F[xm]≈16h(f(xm+1)+4f(xm+1/2)ddd-4f(xm-1/2)-f(xm-1)),
if Simpson’s rule is adopted. We call the scheme using substitution (14) an approximate solution interpolation scheme (ASI), to tell from SI using F(x) explicitly.Throughout the examples, we assume that the coefficienta is fixed at 1 and the function f is given as
(15)f(x)=sin(xε),
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 h=0.01 and k=0.005, respectively.Example 1.
Nonhomogeneous linear convection equation:(16)ut+aux=f(x).
This equation has a general solution,
(17)u=ϕ(x-at)+1aF(x),
where ϕ is an arbitrary differentiable function. Choosing ϕ(x)=C0+C1x+C2x2, we solve the interpolation equations,
(18)uln=ϕ(xl-atn)+1aF(xl),l=m-1,m,m+1.
Considering that (16) does not involve t, we can set tn=0 in (18) without losing generality. Then the solution is
(19)C0=(1-xmD+(xm)22D2)(umn-1aF(xm)),C1=(D-xm2D2)(umn-1aF(xm)),C2=D2(umn-1aF(xm)).
Plugging this result into
(20)umn+1=ϕ(xm-ak)+1aF(xm)
gives the solution interpolation scheme,
(21)umn+1=umn+(-akD+a2k22D2)(umn+1aF(xm)).
The scheme might be similar in appearance to the Lax-Wendroff scheme,
(22)umn+1=umn+(-akD+a2k22D2)umn+(k-ak22D)f(xm).
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.Numerical solutions forut+aux=sin(x/ε): a solution u=(x-t)(sin(x-t)+2)-1-εcos(x/ε) is used to generate initial/boundary conditions and measure errors.
(a)
ε
=
1
(b)
ε
=
0.01Example 2.
Linear convection equation with a varying coefficient:(23)ut+aux=f(x)u.
A general solution for (23) is
(24)u=exp(1aF(x))ϕ(x-at).
The choice of the same function ϕ(x)=C0+C1x+C2x2 as in the previous example eventually leads to the scheme
(25)umn+1=αum+1n+βumn+γum-1n,
where the coefficients α, β, and γ are determined as
(26)α=aλ(-1+aλ)exp(-hDrF[xm]),β=1-a2λ2,γ=aλ(1+aλ)exp(hDlF[xm]).
Here λ denotes a Courant number k/h. 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.Numerical solutions forut+aux=sin(x/ε)u: a solution u=exp(-εcos(x/ε)-x+t)(x-t)2+1 is used to generate initial/boundary conditions and measure errors.
(a)
ε
=
1
(b)
ε
=
0.01Example 3.
Quasilinear convection equation with a varying coefficient:(27)ut+aux=f(x)u2.
The equation has a general solution,
(28)u=(ϕ(x-at)-1aF(x))-1.
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
(29)umn+1=(um+1numnum-1n)×(um+1num-1n+(akS+a2k22S2)umn+um+1numnum-1nddd×(-akD+a2k22D2)F(xm))-1.
Here, we use an additional difference operator S as
(30)Sumn=um+1numn-umnum-1n2h,S2umn=um+1numn-2umnum-1n+umnum-1nh2,
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(x/ε), ε=0.01.Numerical solutions forut+aux=sin(x/ε)u2: a solution u=(x2-2xt+t2+1)(1+ε(x2-2xt+t2+1)cos(x/ε))-1 is used to generate initial/boundary conditions and measure errors.
(a)
ε
=
1
(b)
ε
=
0.01


## 4. Solution Interpolation for Wave Equations

In this section, we apply the solution interpolation method to the wave equation (2) with variable coefficient c(x). It is known that (2) has two families of particular solutions as follows [16]. Solutions with even powers of t are
(31)ul(x,t)=∑k=0l‍ϕk(x)t2(l-k),l=0,1,2,…,
where the function ϕ=ϕk(x) is defined by the recurrence relations
(32)ϕ0(x)=A0x+B0,ϕk+1(x)=Akx+Bk+2(l-k)(2l-2k-1)×∫ax‍(x-s)ϕk(s)c(s)ds,k=1,2,…,l-1.
Here Ak, Bk, and a are arbitrary constants. Similarly, there is another family of particular solutions with odd powers of t as
(33)ul(x,t)=∑k=0l‍φk(x)t2l-2k+1,
where the function φ=φk(x) is defined by the recurrence relations:
(34)φ0(x)=A0x+B0,φk+1(x)=Akx+Bk+2(l-k)(2l-2k+1)×∫ax‍(x-s)φn(s)c(s)ds,k=1,2,…,l-1.
Again, Ak, Bk, and a are arbitrary constants.Among the above solutions, let us choose a four-parameter solution,(35)u(x,t)=C0+C1x+C2t+C3t2+2C3∫0x‍(x-s)c(s)ds,
out of which we construct a five-point scheme for the equation. Since (2) does not explicitly depend on time t, we can assume tn=0 without loss of generality. Interpolating the numerical data um-1n, umn, um+1n, and umn-1 on the grid with solution (35) determines the parameters as
(36)C0=(F(xm)D2umn+xmD2F(xm)Dumnd-umnD2F(xm)-xmD2umnDF(xm))×(D2F(xm))-1,C1=D2umnDF(xm)-D2F(xm)DumnD2F(xm),C2=k2D2umn+2(umn-umn-1)D2F(xm)2kD2F(xm),C3=D2umn2D2F(xm),
where F(x)=∫0x‍((x-s)/c(s))ds. We plug this result to umn+1=u(xm,tn) in (35) and then eventually obtain the scheme:
(37)umn+1=2umn-umn-1+k2D2umnD2F(xm),
after simplification. Note that we do not have to integrate (x-s)/c(s) and find F(x) explicitly. Instead, by using numerical integration method, we can replace the factor involving F(xm) with values of c(x) on some points. For example, adopting Simpson’s rule gives
(38)D2F(xm)≈3(1c(xm-1/2)+1c(xm)+1c(xm+1/2)).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:
(39)umn+1=2umn-umn-1+k2αmD2umn,
where
(40)αm=c(xm):CF,αm=(3c(xm-1/2)c(xm)c(xm+1/2))×((xm+1/2)c(xm-1/2)c(xm)+c(xm-1/2)c(xm+1/2)dd.+c(xm)c(xm+1/2))-1:ASI.Now let us show the performances of ASI in comparison to CF, focusing especially on the cases where the propagation speedc(x) changes rapidly or discontinuously.Example 4.
Wave equation with exponentially decreasing speed:(41)utt=exp(-xε)uxx.
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 x-axis, one can see that ASI yields much better result while CF loses its accuracy quickly as time goes.Time evolution of the errors in numerical solutions forutt=exp(-x/ε)uxx,0≤x≤1. The exact solution u=(1+x)t2 + xt+x+2(-ε2+2ε3+(-ε+ε2)x+(ε2-2ε3+ε2x)exp(x/ε)) is used to generate the initial/boundary conditions and measure the errors.
(a)
ε
=
1
(b)
ε
=
0.05Example 5.
Wave equation with oscillating speed:(42)utt=12+sin(-x/ε)uxx.
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).Time evolution of the errors in numerical solutions forutt=(2+sin(-x/ε))-1uxx,0≤x≤1. The exact solution u=(1+x)t2+xt+x+2(2ε3+εx+x2+x3/3-2ε3cos(x/ε)-ε2(1+x)sin(x/ε)) is used to generate the initial/boundary conditions and measure the errors.
(a)
ε
=
1
(b)
ε
=
0.01Example 6.
Wave equation with discontinuous speed:(43)utt=c(x)uxx,
where c(x) is a step function as
(44)c(x)={1ifx<0,2ifx≥0.
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.Snapshots of numerical solutions for (43) at t=0.3. In (a), the exact solution is u=exp(x-t)-exp(-x+t)/2 if x<0 and u=exp(-x/2-t) − exp(x/2+t)/2 otherwise. In (b), the exact solution is u=exp(x+t) − exp(-x-t) if x<0 and u=exp(x/2-t)+exp(-x/2+t) otherwise.
(a)
(b)


## 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 replaceumn and umn+1 in a scheme with g(tn)eimθ and g(tn+k)eimθ, respectively. The scheme is stable as long as the amplification factor g(t) satisfies
(45)|g(tn+k)|≤|g(tn)|for-π≤θ≤π,k≥0.
In fact, this stability condition can be relaxed as
(46)|g(tn+k)|≤|g(tn)|+O(k)for-π≤θ≤π,k≥0.
We use the above substitution and rewrite scheme (21) as
(47)g(tn+k)eimθ=g(tn)eimθ-aλ2(g(tn)ei(m+1)θ-g(tn)ei(m-1)θ)+a2λ22(g(tn)ei(m+1)θ-2g(tn)eimθdddddd+g(tn)ei(m-1)θ)+(-akD+a2k22D2)F(xm),
where λ=k/h. This can be simplified and reorganized as
(48)g(tn+k)=αg(tn)+β,
where
(49)α=1-a2λ2(1-cosθ)-iaλsinθ,β=e-imθ(-akD+a2k22D2)F(xm).
If we assume |α|≤1, then
(50)|g(tn+k)|=|αg(tn)+β|≤|α||g(tn)|+|β|≤|g(tn)|+|β|,
and since β=O(k), the scheme would be stable. By the way,
(51)|α|2=(1-2a2λ2sin2(θ2))2+(aλsinθ)2=1-4a2λ2(1-a2λ2)sin4(θ2).
Therefore we see |α|≤1⇔(1-a2λ2)≥0, and this implies that the stability condition for scheme (21) is
(52)|aλ|≤1.
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 vt+avx=f(x) after substitution v=1/u. Similarly, scheme (25) can be converted to another solution interpolation scheme for the homogeneous advection equation:
(53)vt+avx=0,
by substitution of v=exp(-1/aF(x))u. Note that this substitution does not affect the stability region. The solution interpolation scheme for (53) is
(54)umn+1=umn+(-akD+a2k22D2)umn,
which is nothing but the Lax-Wendroff scheme and is stable if |aλ|≤1. Therefore, scheme (25) has the same numerical stability too.For the variable coefficient wave equation (2), the general procedure is to regard c(x) as a frozen coefficient for each x, t 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
(55)|c(x)|λ≤1orλ≤1|c(x)|.
In the case of the solution interpolation scheme (37), one can show that the condition becomes
(56)λ≤1|c(x)|+O(h2),
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 c(x) is differentiable. However, the analysis is valid regardless of its differentiability since the scheme approximates c(x) 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)utt-a2uxx=-af′(x),uttt+a3uxxx=a2f′′(x),utttt-a4uxxxx=-a3f′′′(x),uttttt+a5uxxxxx=-a4f′′′′(x)⋮
Suppose u(xm,tn+1) is the exact value of a solution and umn+1 is the value obtained from scheme (21). Application of the Taylor expansion to u(xm,tn+1)-umn+1 at (xm,tn) gives
(58)u(xm,tn+1)-umn+1=(ut+aux-f(x))k+12(utt-a2uxx+af′(x))k2+16(uttt+aλ2uxxx-1λ2f′′(x))k3+124(utttt-a2λ2uxxxx+aλ2f′′′(x))k4+1120(uttttt+aλ4uxxxxx-1λ4f′′′′(x))k5+⋯.
Note that the first two terms on the right hand side are always zero. This indicates that the error of the scheme is O(k3,h2k). Now, also suppose λ=1/a 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/a, which implies that umn+1 which we found is exact.It is no wonder this happens sinceumn+1 is on the same characteristic line with umn when λ=1/a. In this case, change of u along this line can be traced by solving a simple ordinary differential equation and especially at (xm,tn+1), we obtain
(59)u(xm,tn+1)=u(xm-1,tn)+∫xm-1xm‍f(x)dx.
This is in fact exactly what we get from scheme (21) with h=ak,
(60)umn+1=um-1n+DlF(xm).
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.



---

*Source: 546031-2013-12-08.xml*