2021 UNC Greensboro PDE Conference, Electronic Journal of Differential Equations, Conference 26 (2022), pp. 97–113. ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu SEMI-LAGRANGIAN FORWARD METHODS FOR SOME TIME-DEPENDENT NONLINEAR PARTIAL DIFFERENTIAL EQUATIONS DANIEL X. GUO Abstract. In this article, we study one-step Semi-Lagrangian forward method for computing the numerical solutions of time-dependent nonlinear partial dif- ferential equations with initial and boundary conditions in one space dimen- sion. Comparing with classic Semi-Lagrangian method, this method is more straight forward to analyze and implement. This method is based on La- grangian trajectory from the departure points (regular nodes) to the arrival points by Runge-Kutta methods. The arrival points are traced forward from the departure points along the trajectory of the path. Most likely the ar- rival points are not on the regular grid nodes. However, it is convenient to approximate the high order derivative terms in spatial dimension on regular nodes. The convergence and stability are studied for the explicit methods. The numerical examples show that those methods work very efficient for the time-dependent nonlinear partial differential equations. 1. Introduction Differential equations in mathematical models are very important part, and de- pending on how many variables are involved, the formulations of ordinary differ- ential equations (ODEs), and of partial differential equations (PDEs). For some equation, we can use the traditional techniques to obtain an explicit or implicit solutions. However, most realistic mathematical models cannot be solved in this way; instead, they must be dealt with by computational methods that provide approximate solutions in some sense, [6, 12, 13]. A lot of numerical methods have been investigated to produce the approxima- tions of the desired solutions. Semi-Lagrangian methods have been introduced at the beginning of the eighties [21] and since then have been extensively used in the numerical simulation of models for weather forecast, oceanography, and for more general problems in fluid dynamics [4, 18, 20, 24]. The basic idea of the semi- Lagrangian method is to construct the numerical solution over a set of grid points by advancing the characteristics, or the trajectory path. In the article [17], a Semi-Lagrangian backward method (traditional method) was investigated for its stability and convergence. However there are two disadvantages. 2020 Mathematics Subject Classification. 35Q35, 65M12, 76M20. Key words and phrases. Semi-Lagrangian forward methods; trajectory; Runge-Kutta method; time-dependent nonlinear partial differential equations. ©2022 This work is licensed under a CC BY 4.0 license. Published August 25, 2022. 97 98 D. GUO EJDE-2022/CONF/26 One is to compute the departure points by an implicit method. In each iteration at the approximated departure point, a same order interpolation method comparing to the iteration method is needed on no-regular grids. The computation cost is increasing according to the order. The another one is to approximate the terms with spatial derivatives at the irregular grid points. The approximations are not based on the regular nodes since the departure points most likely is not on the regular node. Recently, more applications of semi-Lagrangian method are reported. Typically they are for the Navier-Stokes equations [7], shallow-water equations [9, 1, 19], Advection-diffusion equations [11, 2], and Burgers equation [3]. Some even tried to use semi-Lagrangian forward method [10, 19]. The advantages of semi-Lagrangian forward method are the overall order of approximation, less computations from interpolation, and standard discretization of terms involving spatial derivatives. In this article, we study the stability and convergence of semi-Lagrangian forward method. We consider the time-dependent nonlinear differential equation ut + uux = f(t, x, u, ux, ...), t > 0, x ∈ Ω ⊂ Rn, with initial-value and boundary conditions. Where n ≥ 1 could be any positive integer. On the Lagrangian trajectory path, those equations could be considered as ODEs. To obtain a numerical approximation of the solution, the arrival and depar- ture points at tn+1 and t = tn are needed. One could be specified on regular nodes, but the other one needs to be computed. For example, if the arrival points are on the grid, then the departure points must be calculated [17]. In this article, the depar- ture points are fixed on regular grids, but the arrival points must be computed and most likely they are not on the regular grid. Also Runge-Kutta methods could be applied to different types of partial differential equations [5, 8, 22, 27]. In the paper [16], the truncation error of the first-order Semi-Lagrangian method was studied. The complete numerical analysis including stability and convergence is presented in this paper for the first-order and the second-order Semi-Lagrangian methods. The numerical analysis and numerical example showed that it is possible to construct overall order-two/three one-step difference method for the time-dependent partial differential equations. This article is organized as follows: Section 2 shows the specified partial dif- ferential equations with given initial-value and boundary conditions to be studied; In Section 3, the arrival point is computed; An one-step difference method is con- structed and analyzed in Section 4; In Section 5, the Semi-Lagrangian methods are investigated; Numerical examples are presented in Section 6; The final section is the summary and the consideration of future work. 2. Time-dependent partial differential equations In this article, we will focus on the partial differential equation du(t, x) dt = f(t, x, u), t > 0, x ∈ Ω, with initial-value condition u(0, x) = u0(x), x ∈ Ω, where Ω is a domain in one, two, or three dimensions. The boundary condition is provided for u(t, x) on ∂Ω if it is needed. EJDE-2018/CONF/26 SEMI-LAGRANGIAN FORWARD METHODS FOR SOME PDES 99 The total derivative reads d dt = ∂ ∂t + u · ∂ ∂x , where u is the velocity field and x is the spatial variable. In this article, the dimension one is considered for demonstration. It is more complicated for higher dimensions and it will be reported later. For examples, the Kortweg-de-Vries (KdV) equation ut + uux + δ2uxxx = 0, with boundary and initial conditions. Where δ is a given constant. Also the Navier- Stokes equation reads ut + (u · ∇)u = f + ν∆u−∇p, div u = 0, with boundary and the initial conditions. Even the reformulated Shallow-Water equations [14, 15] have the similar form. To focus on the proposed method, we will study the general time-dependent differential equation du dt = ∂u ∂t + u ∂u ∂x = g(t, x, u), 0 < t < T, a ≤ x ≤ b u(0, x) = u0(x), a ≤ x ≤ b (2.1) with boundary conditions about x and g(t, x, u) = f(t, x) +Lu, where L is a linear derivative operator in spatial x, for example Lu = cuxx, or cuxxx, where c is given constant. The one-step forward explicit methods for the problem (2.1) are proposed. Also we study the stability and convergence of these methods for problem (2.1). Theorem 2.1. If g(t, x, u) is a continuous function of t and x, and satisfies the Lipschitz condition in u in the region 0 ≤ t ≤ T , a ≤ x ≤ b, −∞ < u < ∞, and the first order derivative of u0(x) is continuous on [a, b], then there exists at least a differentiable solution u(t, x) of the initial-value boundary problem (2.1). This result can be found in many books on partial differential equations, such as [26, page 282]. 3. Computing the arrival points It is well-known that the exact solution of the problem (2.1) for most cases is hard to find. The approximations to u will be computed first at the given points, called grid points, in the domain [0, T ] × [a, b]. At other points in the domain, the approximate solutions are calculated from the interpolation of the approximate solutions on the grid points. Now assume that the grid points are equally distributed throughout the domain [0, T ] × [a, b]. Let N and M be two positive integer numbers. The grid points are formed by calculating tk = kτ, for k = 0, 1, 2, . . . , N ; xi = a+ ih, for i = 0, 1, 2, . . . ,M ; with the step sizes τ = T/N and h = (b− a)/M . Assume that all solutions u on (tk, x k i ) for i = 0, 1, . . . ,M are known, and we need to find the solutions u on (tk+1, x k+1 i ) for i = 0, 1, 2, . . . ,M . From the view of 100 D. GUO EJDE-2022/CONF/26 Lagrangian trajectory, the particles on arrival points at t = tk+1 shall come from some points, called departure points at t = tk along the trajectory of path. To compute the arrival point xk+1 A at t = tk+1 for the departure point xkD, we consider the equation dx dt = u. (3.1) T X x D1 k x D2 k x D1 k+1 x D2 k+1 x A1 k+1 x A2 k+1 Figure 1. Departure and arrival points on grid. Note that the arrival point xk+1 A may not be any of xi, i = 0, 1, . . . ,M at t = tk+1 as shown in Figure 1. Then the arrival point xk+1 A is calculated by integrating the equation along the trajectory from the departure point xkD as xk+1 A = xkD + ∫ tk+1 tk u dt, xkD ∈ {xi : i = 0, 1, . . . ,M}. It is obvious that Euler’s method will provide the first order approximation as xk+1 A ≈ xkD + τu(tk, x k D). This is an explicit equation about xk+1 A . It is possible to obtain a higher order approximation by combining the computation of solution u(tk+1, x k+1 A ) at the time t = tk+1, like, the modified Euler’s method. 4. One-step forward difference method For the problem (2.1), along the trajectory of the path, we propose a general one-step forward difference method in the form w0 D = u0(x0D), x0D ∈ {xi : i = 0, 1, . . . ,M}, xk+1 A = xkD + τwkD, wk+1 A = wkD + τφ(tk, x k D, w k D, τ, h), wk+1 D Calculated from linear interpolation of wk+1 A (4.1) for k = 0, 1, . . . , N−1, where φ(t, x, u, τ, h) is a continuous function of t, x, τ , and h and satisfies the Lipschitz condition in u. xkD is the departure point on the regular grids; xk+1 A is the arrival point and most likely is not any of the regular grids like EJDE-2018/CONF/26 SEMI-LAGRANGIAN FORWARD METHODS FOR SOME PDES 101 showed in the Figure 1, like xk+1 A1 or xk+1 A2 ; wkD and wk+1 A are approximations of u(tk, x k D) and u(tk+1, x k+1 A ) respectively for k = 0, 1, . . . , N − 1. Note that wk+1 D is calculated from linear interpolation of wk+1 A . If we consider a higher order difference method, a higher order interpolation is needed. Assuming that τ = O(h) (or τ = αh, α > 0 a constant), and let u(t, x) be the unique solution of the initial-value problem (2.1). We define the local truncation error for each time step as follows. Definition 4.1. The local truncation error of the difference method (4.1) to the problem (2.1) is defined as εk+1(τ) = uk+1 A − [ukD + τφ(tk, x k D, u k D, τ, h)] τ = uk+1 A − ukD τ − φ(tk, x k D, u k D, τ, h), where ukD = u(tk, x k D) and uk+1 A = u(tk+1, x k+1 A ) on the departure and arrival points respectively. We need two lemmas whose proofs can be found in any numerical analysis book, for example [6]. Lemma 4.2. For all t ≥ −1 and any positive m, we have 0 ≤ (1 + t)m ≤ emt. Lemma 4.3. If s and t are positive real numbers, {ak}nk=0 is a sequence satisfying a0 ≥ −t/s, and ak+1 ≤ (1 + s)ak + t, for each k = 0, 1, 2, . . . , n− 1, then ak+1 ≤ e(k+1)s ( a0 + t s ) − t s . T X x D1 k x D2 k x D1 k+1 x D2 k+1 x Au1 k+1 x Av1 k+1 x Au2 k+1 x Av2 k+1 Figure 2. Different arrival points on grid. Theorem 4.4. The initial-value boundary problem (2.1) is approximated by a one- step difference method (4.1). Suppose that there exists a number τ0 > 0 such that φ(t, x, w, τ, h) is continuous and satisfies a Lipschitz condition in all variables with Lipschitz constant L on the set D = {(t, x, w, τ, h)| 0 ≤ t ≤ T, a ≤ x ≤ b,−∞ < w <∞, 0 ≤ τ ≤ τ0, 0 ≤ h ≤ ατ0}. Then 102 D. GUO EJDE-2022/CONF/26 (i) the difference method (4.1) is stable; (ii) the difference method (4.1) is convergent if and only if it is consistent, which is equivalent to φ(t, x, w, 0, 0) = g(t, x, w), for all 0 ≤ t ≤ T, a ≤ x ≤ b,−∞ < w <∞; (iii) if a function ε exists and, for each k = 0, 1, . . . , N , the local truncation error εk(τ) satisfies |εk(τ)| ≤ ε(τ) whenever 0 ≤ τ ≤ τ0, then there exists a constant C such that max D |u(tk, x k D)− wkD| ≤ eTL max D |u0(x0D)− w0 D|+ ε(τ) + Ch L (eTL − 1) ≤ ε(τ) + Ch L (eTL − 1), for k = 1, 2, . . . , N. Proof. Let u and v be two solutions of the initial-value problem (1) with the initial values u0 and v0 respectively. There exist three constants M , K1, and K2 such that |u|, |v| ≤M, |∂u ∂x |, |∂u ∂t |, |∂v ∂x |, |∂v ∂t | ≤ K1, |∂ 2u ∂x2 |, |∂ 2v ∂x2 | ≤ K2. (4.2) and τK1 ≤ 1/2. (4.3) (i) Suppose that {ukD}Nk=0 and {vkD}Nk=0 satisfy the difference equation (4.1) with the initial-values u0 and v0 respectively. Let EkD = ukD − vkD for the point (tk, x k D) and Ek = max{|EkD|, for all regular grid points (departure points) at t = tk}. Without lose generality, let the arrival points be as in the Figure 2. If xk+1 D2 is located in between xk+1 Au2 and xk+1 Av2 , the procedure is very similar. Then uk+1 D2 = ωuu k+1 Au2 + (1− ωu)uk+1 Au1 , vk+1 D2 = ωvv k+1 Av2 + (1− ωv)vk+1 Av1 , (4.4) where ωu = xk+1 D2 − x k+1 Au1 xk+1 Au2 − xk+1 Au1 , ωv = xk+1 D2 − x k+1 Av1 xk+1 Av2 − xk+1 Av1 , and 0 ≤ ωu, ωv ≤ 1. Then from (4.4), |uk+1 D2 − v k+1 D2 | = |ωuuk+1 Au2 + (1− ωu)uk+1 Au1 − ωvvk+1 Av2 − (1− ωv)vk+1 Av1 | ≤ |ωuuk+1 Au2 − ωuvk+1 Av2 + (1− ωu)uk+1 Au1 − (1− ωu)vk+1 Av1 | + |ωuvk+1 Av2 − ωvvk+1 Av2 + (1− ωu)vk+1 Av1 − (1− ωv)vk+1 Av1 | ≤ ωu|uk+1 Au2 − vk+1 Av2 |+ (1− ωu)|uk+1 Au1 − vk+1 Av1 |+ |ωu − ωv‖vk+1 Av2 − vk+1 Av1 |. (4.5) For each departure point, for example xkD1, it reads from (4.1), uk+1 Au1 = ukD1 + τφ(tk, x k D1, u k D1, τ, h) and vk+1 Av1 = vkD1 + τφ(tk, x k D1, v k D1, τ, h). EJDE-2018/CONF/26 SEMI-LAGRANGIAN FORWARD METHODS FOR SOME PDES 103 Then |uk+1 Au1 − vk+1 Av1 | ≤ |ukD1 − vkD1|+ τ |φ(tk, x k D1, u k D1, τ, h)− φ(tk, x k D1, v k D1, τ, h)| ≤ |ukD1 − vkD1|+ τL|ukD1 − vkD1| ≤ (1 + τL)|ukD1 − vkD1| ≤ (1 + τL)|EkD1| ≤ (1 + τL)Ek. (4.6) Similarly, |uk+1 Au2 − vk+1 Av2 | ≤ (1 + τL)Ek. (4.7) From the definition of ωu and ωv and (4.1), we have |ωu − ωv| = ∣∣xk+1 D2 − x k+1 Au1 xk+1 Au2 − xk+1 Au1 − xk+1 D2 − x k+1 Av1 xk+1 Av2 − xk+1 Av1 ∣∣ = ∣∣ h− τukD1 h+ τ(ukD2 − ukD1) − h− τvkD1 h+ τ(vkD2 − vkD1) ∣∣ = |hτ(ukD2 − vkD2) + τ2(ukD1v k D2 − vkD1u k D2)| |h+ τ(ukD2 − ukD1)‖h+ τ(vkD2 − vkD1)| . From (4.2), it follows that h(1− τK1) ≤ |h+ τ(ukD2 − ukD1)| ≤ h(1 + τK1). Similarly, h(1− τK1) ≤ |h+ τ(vkD2 − vkD1)| ≤ h(1 + τK1). Then with help of (4.2) and (4.3), |ωu − ωv| ≤ τ h(1− τK1)2 [ |ukD2 − vkD2|+ τ h (|ukD1 − vkD1‖vkD2|+ |ukD2 − vkD2‖vkD1|) ] ≤ τ h(1− τK1)2 (1 + 2αM)Ek, (where τ = αh) ≤ 4 τ h (1 + 2αM)Ek. (4.8) Putting (4.6), (4.7), and (4.8) in (4.5), with the help of (4.2), yields |uk+1 D2 − v k+1 D2 | ≤ ωu(1 + τL)Ek + (1− ωu)(1 + τL)Ek + 4 τ h (1 + 2αM)EkhK1 ≤ [1 + τL+ 4τK1(1 + 2αM)]Ek (4.9) i.e. |Ek+1 D | ≤ [1 + τL+ 4τK1(1 + 2αM)]Ek, Ek+1 ≤ (1 + τL1)Ek where L1 = L+ 4K1(1 + 2αM). Therefore, from Lemmas 4.2 and 4.3, Ek ≤ (1 + τL1)kE0 ≤ ekτL1E0 ≤ eTL1E0, where k ≤ N = T/τ . So, the difference method (4.1) is stable. (ii) Let φ(t, x, w, 0, 0) = g(t, x, w). From the assumptions of φ, g satisfies the conditions of Theorem 2.1, then the differential equation dv dt = g(t, x, v), 0 < t < T, a ≤ x ≤ b, 104 D. GUO EJDE-2022/CONF/26 v(0, x) = u0(x), a ≤ x ≤ b, has a solution v(t, x). The numerical solution satisfies z0D = u0(x0D), x0D ∈ {xi : i = 0, 1, . . . ,M}, xk+1 A = xkD + τ zkD, zk+1 A = zkD + τφ(tk, x k D, z k D, τ, h), zk+1 D Calculated from linear interpolation of zk+1 A (4.10) for k = 0, 1, . . . , N − 1. By the mean value theorem and referring to the Figure 1, v(tk+1, x k+1 A2 ) = v(tk, x k D2) + τg(tk + ξτ, xkD2 + ηh, v(tk + ξτ, xkD2 + ηh)), for some constants ξ and η (0 ≤ ξ, η ≤ 1). Similarly, define EkD = zkD − v(tk, x k D) and Ek = max |EkD| on the regular grid points (departure points) at t = tk. Subtract this from the numerical solution, we have zk+1 A2 − v(tk+1, x k+1 A2 ) = zkD2 − v(tk, x k D2) + τ [φ(tk, x k D2, z k D2, τ, h)− g(tk + ξτ, xkD2 + ηh, v(tk + ξτ, xkD2 + ηh))] = EkD2 + τ [φ(tk, x k D2, z k D2, τ, h)− φ(tk, x k D2, v(tk, x k D2), τ, h) + φ(tk, x k D2, v(tk, x k D2), τ, h)− φ(tk, x k D2, v(tk, x k D2), τ, 0) + φ(tk, x k D2, v(tk, x k D2), τ, 0)− φ(tk, x k D2, v(tk, x k D2), 0, 0) + φ(tk, x k D2, v(tk, x k D2), 0, 0)− g(tk + ξτ, xkD2 + ηh, v(tk + ξτ, xkD2 + ηh))], Thanks to the Lipschitz condition of φ, |φ(tk, x k D2, z k D2, τ, h)− φ(tk, x k D2, v(tk, x k D2), τ, h)| ≤ L|zkD2 − v(tk, x k D2)| ≤ L|EkD2|, |φ(tk, x k D2, v(tk, x k D2), τ, h)− φ(tk, x k D2, v(tk, x k D2), τ, 0)| ≤ Lh, |φ(tk, x k D2, v(tk, x k D2), τ, 0)− φ(tk, x k D2, v(tk, x k D2), 0, 0)| ≤ Lτ, and |φ(tk, x k D2, v(tk, x k D2), 0, 0)− g(tk + ξτ, xkD2 + ηh, v(tk + ξτ, xkD2 + ηh))| ≤ |g(tk, x k D2, v(tk, x k D2))− g(tk + ξτ, xkD2, v(tk, x k D2))| + |g(tk + ξτ, xkD2, v(tk, x k D2))− g(tk + ξτ, xkD2 + ηh, v(tk, x k D2))| + |g(tk + ξτ, xkD2 + ηh, v(tk, x k D2)) − g(tk + ξτ, xkD2 + ηh, v(tk + ξτ, xkD2 + ηh))| ≤ Lξτ + Lηh+ L|v(tk, x k D2)− v(tk + ξτ, xkD2 + ηh))| ≤ Lξτ + Lηh+ LK1ξτ + LK1ηh ≤ Lτ + Lh+ LK1τ + LK1h ≤ L(1 +K1)(τ + h). Then |Ek+1 A2 | ≤ |E k D2|+ Lτ |EkD2|+ Lτh+ Lτ2 + L(1 +K1)τ(τ + h) ≤ (1 + τL)Ek + L(2 +K1)τ(τ + h). (4.11) EJDE-2018/CONF/26 SEMI-LAGRANGIAN FORWARD METHODS FOR SOME PDES 105 From the Figure 1, define a constant ω as ω = xk+1 D2 − x k+1 A1 xk+1 A2 − x k+1 A1 . Then zk+1 D2 = ωzk+1 A2 + (1− ω)zk+1 A1 and there exists a constant η such that v(tk+1, x k+1 D2 ) = ωv(tk+1, x k+1 A2 ) + (1− ω)v(tk+1, x k+1 A1 ) + 1 2! ∂2v(tk+1, η) ∂x2 (xk+1 D2 − x k+1 A1 )(xk+1 D2 − x k+1 A2 ), Then |zk+1 D2 − v(tk+1, x k+1 D2 )| ≤ ω|zk+1 A2 − v(tk+1, x k+1 A2 )| + (1− ω)|zk+1 A1 − v(tk+1, x k+1 A1 )|+ 1 2 K2|xk+1 D2 − x k+1 A1 ‖x k+1 D2 − x k+1 A2 | ≤ ω|Ek+1 A2 |+ (1− ω)|Ek+1 A1 |+ 1 2 K2|h+ τv(tk, x k D1)|τ |v(tk, x k D2)| ≤ ω|Ek+1 A2 |+ (1− ω)|Ek+1 A1 |+ 1 2 K2M(1 + αM)τh. (4.12) Thanks to (4.11) and the similar result for Ek+1 A1 , it follows that Ek+1 ≤ (1 + τL)Ek + L(2 +K1)τ(τ + h) + 1 2 K2M(1 + αM)τh ≤ (1 + τL)Ek + L1(τ2 + τh), where L1 = L(2 +K1) + 1 2K2M(1 + αM). Thanks to Lemma 4.3 we have Ek ≤ ekτLE0 + L1 L (τ + h)(ekτL − 1), Ek ≤ eTLE0 + L1 L (τ + h)(eTL − 1). This converges to zero as τ and h (E0 = 0) go to zero. So the numerical solution converges to the solution of initial-value boundary problem (2.1) if f = g. Therefore the numerical solution from the difference method (4.1) is convergent. On the other hand, if we have convergence, then u and v are the same. Suppose f and g differ at some point.Then we consider the initial-value boundary prob- lem starting from that point. Two solutions should be different. This leads to a contradiction. (iii) Let EkA, EkD, and Ek are defined as in (i). Then from the definition of the local truncation error, we have u(tk+1, x k+1 A ) = u(tk, x k D) + τφ(tk, x k D, u(tk, x k D), τ, h) + τεk+1(τ). Subtracting this from the difference method yields Ek+1 A = EkD + τ [φ(tk, x k D, w k D, τ, h)− φ(tk, x k D, u(tk, x k D), τ, h)]− τεk+1(τ). therefore, |Ek+1 A | ≤ |EkD|+ τL|wkD − u(tk, x k D)|+ τ |εk+1(τ)| ≤ (1 + τL)|EkD|+ τ |εk+1(τ)|. 106 D. GUO EJDE-2022/CONF/26 Using a similar technique as in (4.11) and (4.12), with the assumption for εk(τ) leads to Ek+1 ≤ (1 + τL)Ek + τε(τ) + 1 2 K2M(1 + αM)τh. Then thanks Lemma 4.3, we have Ek ≤ ekτL(E0 + ε(τ) + CH L )− ε(τ) + Ch L , i.e., Ek ≤ eTLE0 + ε(τ) + Ch L (eTL − 1), for k = 1, 2, . . . , N and C = 1 2K2M(1 + αM). � 5. Semi-Lagrangian forward methods Now we study the modified semi-Lagrangian forward Euler method for the initial- value boundary problem (2.1), w0 D = u0(x0D), x0D ∈ {xi : i = 0, 1, . . . ,M}, k1 = g(tk, x k D, w k D), w̄ = wkD + τk1, xk+1 A = xkD + 1 2 τ(wkD + w̄), k2 = g(tk+1, x k+1 A , w̄), wk+1 A = wkD + 1 2 τ(k1 + k2), wk+1 D Calculated from second-order interpolation of wk+1 A (5.1) for k = 0, 1, . . . , N − 1, for all of arrival points A and k = 0, 1, 2, . . . , N − 1. Theorem 5.1. If the function g(t, x, u) in (2.1) is continuous and satisfies a Lip- schitz condition in all variables with Lipschitz constant L on the set D = {(t, x, w) : 0 ≤ t ≤ T, a ≤ x ≤ b,−∞ < u <∞}. Then (i) The difference method (5.1) is stable; (ii) The difference method (5.1) is convergent. The stability and convergence of (5.1) could be proved following the similar idea with Lipschitz continuity in Theorem 4.4. Theorem 5.2. Suppose h = O(τ) and the function g(t, x, u) in (2.1) is continuous and differentiable. If εk is the local truncation error of (5.1), then ‖εk(τ)‖L∞ = O(τ2). Proof. Let u(t, x) be the unique solution with up to order three continuous bounded partial derivatives on [0, T ]× [a, b], so that for each k = 0, 1, . . . , N − 1, g(tk+1, x k+1 A , wkD + τg(tk, x k D, w k D)) = g(tk, x k D, w k D) + τ ∂g(tk, x k D, w k D) ∂t + τwkD ∂g(tk, x k D, w k D) ∂x EJDE-2018/CONF/26 SEMI-LAGRANGIAN FORWARD METHODS FOR SOME PDES 107 + τg(tk, x k D, w k D) ∂g(tk, x k D, w k D) ∂u +O(τ2), and u(tk+1, x k+1 A ) = u(tk, x k D) + τ du(tk, x k D) dt + 1 2 τ2 d2u(tk, x k D) dt2 +O(τ3) = u(tk, x k D) + τg(tk, x k D, u k D) + 1 2 τ2[ ∂g(tk, x k D, u k D) ∂t + ukD ∂g(tk, x k D, u k D) ∂x + g(tk, x k D, u k D) ∂g(tk, x k D, u k D) ∂u ] +O(τ3). From the definition of local truncation error and the method (15), the local trun- cation error reads εk+1(τ) = u(tk+1, x k+1 A )− u(tk, x k D) τ − 1 2 [ g(tk, x k D, u(tk, x k D)) + g(tk+1, x k+1 A , u(tk, x k D) + τg(tk, x k D, u(tk, x k D))] So it is easy to see that the local truncation error at the (k + 1)th step for the method (5.1) is ‖εk+1(τ)‖L∞ = O(τ2). � Note that the big difference in forward and backward methods is the computing of arrival points and departure point. In backward method, it needs the iteration to find better approximation. Using similar techniques, the semi-Lagrangian method of order-three is con- structed as follows w0 D = u0(x0D), x0D ∈ {xi : i = 0, 1, . . . ,M} k1 = g(tk, x k D, w k D) xk+1 A = xkD + τ wkD + 1 2 τ2k1 k2 = g(tk + 1 2 τ, xD + 1 2 τwkD, w k D + 1 2 τk1) k3 = g(tk+1, x k+1 A , wkD − τk1 + 2τk2) wk+1 A = wkD + 1 6 τ(k1 + 4k2 + k3) wk+1 D Calculated from the third-oder interpolation of wk+1 A (5.2) for k = 0, 1, . . . , N − 1. The numerical analysis is much more complicated be- cause of the need of higher order methods for computing of departure points and interpolation at the departure points. 6. Numerical experiments The main computational effort in applying the Semi-Lagrangian forward meth- ods (SLFW) is the evaluation of g and the interpolation of approximation on depar- ture points. In the second-order method, the local truncation error is O(τ2), and the cost is two functional evaluations per step and the computation of the arrival points. The Semi-Lagrangian method of order three requires three evaluations per step and the computation of the arrival points, and the local truncation error is O(τ3). 108 D. GUO EJDE-2022/CONF/26 -0.4 1 -0.2 0 0.5 0.2u 0.4 0.4 x 0.5 0.6 0.3 t 0.8 0.2 0.1 0 0 Figure 3. Approximated Solution of the initial-value problem. Example 6.1. Consider the problem ut + uux = t+ sin(2πx) + 2πu t cos(2πx), u(0, x) = 0, 0 < t ≤ 1, 0 < x < 1, (6.1) with periodic boundary condition for x. -4e-05 -2e-05 -2e-05 -2e-05 -2e-05 0 0 0 0 0 0 0 0 0 0 2e-05 2e-05 2e-05 2e-05 2e-05 4e-05 4e-05 0 0.1 0.2 0.3 0.4 0.5 t 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 x -4 -3 -2 -1 0 1 2 3 4 5 6 10 -5 Figure 4. Contour lines of the difference between the solution and its approximation. The exact solution for this problem is u(t, x) = 1 2 t 2 + t sin(2πx). Figure 3 shows the graph of the approximate solution of u(t, x) by using the Semi-Lagrangian method of order three with τ = 0.01 and h = 0.01. To focus on the truncation EJDE-2018/CONF/26 SEMI-LAGRANGIAN FORWARD METHODS FOR SOME PDES 109 Table 1. Numerical orders of the Semi-Lagrangian forward method (SLFW1), SLFW2, and SLFW3 as applied to the problem (6.1) τ 0.1 0.05 0.025 0.02 SLFW1 (‖u−u0‖∞ τ ) 1.2189 1.1656 1.1347 1.1270 SLFW2 (‖u−u0‖∞ τ2 ) 2.0234 2.0242 2.0263 2.0262 SLFW3 (‖u−u0‖∞ τ3 ) 2.7279 2.60209 2.5076 2.4823 error from the proposed Semi-Lagrangian forward methods, the periodic boundary condition is supplied. For other boundary conditions, the treatment near boundary is more complicated when the trajectories of the departure points hit the boundary. 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 t 10-6 10-5 10-4 10-3 10-2 10-1 || u -u 0 || SLFW1 SLFW2 SLFW3 Figure 5. Truncation errors: − ∗ − first-order Semi-Lagrangian forward method (SLFW1); −+− second-order method (SLFW2); − ◦ − third-order method (SLFW3). Figure 4 shows the contour lines of the difference between the exact solution and the approximation from the Semi-Lagrangian forward method of order three at the common grid points. The truncation error is increasing as t goes from 0 to 0.5. The contour lines are matching the changes of the function. The larger errors occur when the function has the bigger changes. The maximum error reached at about t = 0.5 for some points of x as we expected. The first-order Semi-Lagrangian forward method, the second-order Semi-Lagrangian forward method, and the Semi-Lagrangian forward method of order three are com- pared with τ = 0.01 and h = 0.005 in the Figure 5. The maximum truncation error presented at the each step of t for those three methods. In the table 1 and the Figure 6, we studied the numerical orders of convergence with h = (1/2)τ . It is clear that the first-order Semi-Lagrangian forward method 110 D. GUO EJDE-2022/CONF/26 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 t 10-6 10-5 10-4 10-3 10-2 10-1 E rr o r ( t) 2 ( t) 3 Error of SL1 Error of SL2 Error of SL3 Figure 6. Scaling of the error of Semi-Lagrangian methods of order one, two and three to the problem (6.1) and the second-order Semi-Lagrangian forward method have the order as we ex- pected. However, the Semi-Lagrangian forward method of order three is not exactly order three, but is close to 2.5. The main difficulty for Semi-Lagrangian method of order three are the calculation of the departure points and the third-order interpo- lation from arrival points(non-uniform grids) to departure points (uniform grids). Example 6.2. Consider the Kortweg-deVries (KdV) equation [25] ut + uux + δ2uxxx = 0 where δ = 0.022 with the initial conditions u(x, 0) = cos(πx), 0 ≤ x ≤ 2 and ux, uxx, and uxxx are periodic on [0, 2] for all t. The KdV equation describes the long time asymptotic behaviour of small but finite amplitude shallow water waves in one dimension. In the equation, u = u(x, t) measures the elevation (the height of water above the equilibrium level) at time t and position x. Two different mechanisms are involved: Non-linearity (uux) and Dispersion (uxxx). The delicate balance between these two effects leads to travelling waves of permanent form, which is called Solitary Waves, or Solitons. Various numerical methods are used to approximate the KdV equations. The comparison of efficiencies of different methods was studied in [23]. The aim here is to show the efficiency of SLFW as an explicit method applying to KdV equations. It is possible to apply higher order (higher than order two) to the KdV equations. By using SLFW, the nonlinear term was perfectly combined in the direction of moving waves. In Figure 7, the first-order SLFW is applied to the KdV equations with the initial condition of u(x, t) = cos(πx) on the interval [0, 2] with ∆x = 0.01 and ∆t = 0.00002. The wave profile at t = 1/π does not produce a shock, but close. EJDE-2018/CONF/26 SEMI-LAGRANGIAN FORWARD METHODS FOR SOME PDES 111 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 -1 -0.5 0 0.5 1 1.5 Figure 7. First-order SLFW method for δ = 0.022: initial values and solutions at t = 1/π and t = 3.6/π. 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 -1 -0.5 0 0.5 1 1.5 2 2.5 Figure 8. Second-order SLFW method for δ = 0.022: initial val- ues and solutions at t = 1/π and t = 3.6/π. When computed at time t = 3.6/π, another shock is also present. The trains of eight solitons are detected, but not well-defined. The results present in Figure 8 are from applying the second-order SLFW to the KdV equations with the initial condition of u(x, t) = cos(πx) on the interval [0, 2] with ∆x = 0.01 and ∆t = 0.00002. The dot line is the initial wave profile. The dash line presents the wave profile at t = 1/π. The wave almost produces a shock and has a noticeable oscillation for x < 1/2. At time t = 3.6/π the profile shows a train of eight well-defined waves which have developed since the early stage of evolution. 112 D. GUO EJDE-2022/CONF/26 7. Conclusions In this article, one-step Semi-Lagrangian forward methods were proposed and studied for the first-order time-dependent nonlinear partial differential equations. The second-order Semi-Lagrangian forward method was proved to be consistent and stable. When τ goes to zero, the numerical solution converges to the solution of the corresponding initial-boundary problem. Theorem 2 in the section 4 and the examples in the section 6 have shown that the Semi-Lagrangian forward methods can be applied to some time-dependent partial differential equations. They are very convenient for implementation. The only disadvantage is that at each time step, spatial linear or higher order interpolation is needed to find the solutions on the departure points. In the example, only one spatial variable was used. However, we believe that Semi-Lagrangian forward methods should work for more spatial variables, but more complicated for numerical analysis. The computing cost will increase a lot, like interpolation for more variables at each step. We will report the similar results for the systems with more spatial variables. Also more terms involving higher order derivative only in spatial variables will be added to the right side function later. More work of applications to Shallow-water equations and Navier-Stokes equations will be reported later. References [1] A. Bourchtein, L. Bourchtein; Semi-Lagrangian semi-implicit time-splitting scheme for the shallow water equations, Int. J. Numer. Meth. Fluids 54 (2007) 453–471. [2] S. Bak; High-order characteristic-tracking strategy for simulation of a nonlinear advection– diffusion equation, Numer Methods Partial Differential Eq. 35 (2019) 1756–1776. [3] S. Bak, P. Kim, D. Kim; A semi-Lagrangian approach for numerical simulation of coupled Burgers’ equations, Commun Nonlinear Sci Numer Simulat 69 (2019) 31–44. [4] J. R. Bates, A. McDonald; Multiply-upstream, semi-Lagrangian advective schemes: analysis and application to a multi-level primitive equation model, Monthly Weather Review, vol 110, (1982) 1831–1842. [5] R. Bermejo, J. Carpio; An adaptive finite element semi-Lagrangian implicit-explicit Runge- Kutta-Chebyshev method for convection dominated reaction-diffusion problems, Applied Nu- merical Mathematics 58 (2009) 16–39. [6] R. L. Burden, J. D. Faires; Numerical Analysis, Thomason Brooks/Cole, 2005. [7] L. Bonaventura, R. Ferretti, L. Rocchi; A fully semi-Lagrangian discretization for the 2D incompressible Navier-Stokes equations in the vorticity-streamfunction formulation, Applied Mathematics and Computation 323 (2018) 132–144. [8] M. P. Calvo, J. de Frutos, J. Novo; Linearly implicit Runge-Kutta methods for advection- reaction-diffusion equations, Applied Numerical Mathematics 37 (2001) 535–549. [9] M. F. Carfora; An unconditionally stable semi-Lagrangian method for the spherical atmo- spherical shallow water equations, Int. J. Numer. Meth. Fluids 34 (2000) 527–558. [10] N. Crouseilles, T. Respaud, E. Sonnendrücker; A forward semi-Lagrangian method for the numerical solution of the Vlasov equation, Computer Physics Communications 180 (2009) 1730–1745. [11] A. Efremov, E. Karepova, V. Shaydurov, A. Vyatkin; A Computational Realization of a Semi- Lagrangian Method for Solving the Advection Equation, Journal of Applied Mathematics Volume 2014, Article ID 610398, 12 pages, https://doi.org/10.1155/2014/610398. [12] C. W. Gear; Numerical Initial Value Problems in Ordinary Differential Equations, Prentice- Hall, Englewood Cliffs, NJ, 1971. [13] D. F. Griffiths, D. J. Higham; Numerical Methods for Ordinary Differential Equations: Ini- tial Value Problems, Springer Undergrauate Mathematics Series, Springer-Verlag London Limited, 2010. EJDE-2018/CONF/26 SEMI-LAGRANGIAN FORWARD METHODS FOR SOME PDES 113 [14] D. X. Guo, J. B. Drake; A global semi-Lagrangian spectral model for the reformulated shallow water equations, Dynamical systems and differential equations, Discrete Contin. Dyn. Syst. suppl., (2003), 375–385. [15] D. X. Guo, J. B. Drake; A global semi-Lagrangian spectral model of the shallow water equa- tions with variable resolution, J. Comput. Phys. 206 (2005) 559–577. [16] D. X. Guo; A semi-Lagrangian Runge-Kutta Method for time-dependent partial differential equations, J. of Applied Analysis & Computation, 3 (3) (2013) 251–263. [17] D. X. Guo; On stability and convergence of semi-Lagrangian methods for the first-order time- dependent nonlinear partial differential equations in 1D, J. of Computational and Applied Mathematics, 324 (2017) 72–84. [18] P. H. Lauritzen, E. Kaas, B. Machenhauer; A massive-conservative semi-implicit semi- Lagrangian limite-area shallow-water model on the sphere, Monthly Weather Review, vol 134, (2006) 1205–1221. [19] M. Mortezazadeh, L. Wang; A high-order backward forward sweep interpolating algorithm for semi-Lagrangian method, Int. J. Numer. Meth. Fluids, 84 (2017) 584–597 [20] R. D. Nair, J. S. Scroggs, F. H. Semazzi; A forward-trajectory global semi-Lagrangian trans- port scheme, J. Comput. Phys., vol 190 (2003), 275–294. [21] A. Robert; A stable numerical integration scheme for the primitive meteorological equations, Atmos. Ocean., vol 19 (1981), 35–46. [22] M. Seaid; An Eulerian-Lagrangian method for parabolic-hyperbolic equations, Applied Nu- merical Mathematics 59 (2009), 745–768. [23] M. Shahrill, M. Chong, & H. Nor; Applying Explicit Schemes to the Korteweg-de Vries Equation, Modern Applied Science, Vol. 9, No. 4 (2015), 200–224. [24] A. Staniforth, J. Cote; Semi-Lagrangian integration schemes for atmospheric models - A review, Monthly Weather Review, vol 119 (1991), 2206–2223. [25] T. R. Taha, M. J. Ablowitz; Analytical and Numerical Aspects of Certain Nonlinear Ecolution Equations. III. Numerical, Korteweg-de Vries Equation, J. Comput. Phys., Vol. 55 (1984), 231–253. [26] R. Temam; Navier-Stokes equations – Theory and numerical analsys, North-Holland Pub- lishing Company, 1977. [27] J. G. Verwer; Explicit Runge-Kutta methods for parabolic partial differential equations, Ap- plied Numerical Mathematics 22 (1996), 359–379. Daniel Guo Department of Mathematics and Statistics, University of North Carolina at Wilming- ton, NC 28403, USA Email address: guod@uncw.edu 1. Introduction 2. Time-dependent partial differential equations 3. Computing the arrival points 4. One-step forward difference method 5. Semi-Lagrangian forward methods 6. Numerical experiments 7. Conclusions References