Barycentric rational collocation methods for a class of nonlinear parabolic partial differential equations

Applied Mathematics Letters 68 (2017) 13–19

Contents lists available at ScienceDirect

Applied Mathematics Letters www.elsevier.com/locate/aml

Barycentric rational collocation methods for a class of nonlinear parabolic partial differential equations✩ Wei-Hua Luoa,b , Ting-Zhu Huanga, *, Xian-Ming Gua , Yi Liu b a

School of Math. Sci., University of Electronic Science and Technology of China, Chengdu, 611731, China Data Recovery Key Laboratory of Sichuan Province, Neijiang Normal University, Neijiang, Sichuan, 641100, China b

article

abstract

info

Article history: Received 9 October 2016 Received in revised form 16 December 2016 Accepted 17 December 2016 Available online 23 December 2016

In this study, based on barycentric rational interpolation functions, three collocation methods are presented for a class of nonlinear parabolic partial differential equations (PDEs). Corresponding algebraic equations are obtained. Numerical examples involving the Fitzhugh–Nagumo equation, Heat equation and Runge’s example are reported to illustrate the effectiveness of the proposed techniques. © 2016 Elsevier Ltd. All rights reserved.

Keywords: Barycentric rational interpolation Collocation method Runge–Kutta method Parabolic PDEs

1. Introduction For a function u : [a, b] → R, given a sequence of points a = x1 < · · · < xn+1 = b, denote by pi the unique polynomial of degree at most d that interpolates u at the d + 1 points xi , . . . , xi+d . Floater et al. [1] have presented a family of barycentric rational interpolants of the form ∑n−d+1 γi (x)pi (x) , (1) r(x) = i=1 ∑n−d+1 γi (x) i=1 i

(−1) where γi (x) = (x−x )···(x−x , i = 1, . . . , n − d + 1, d ≤ n. They proved this function has no real pole, and i i+d ) the interpolation error decreases as O(hd+1 ) when u ∈ C d+2 [a, b]. Berrut et al. [2] analyzed convergence of derivatives of this interpolant and showed the rate O(hd+1−k ) of the kth derivative for k = 1, 2. Later,

✩ This research is supported by 973 Program (2013CB329404), NSFC (61370147), the Fundamental Research Funds for the Central Universities (No. ZYGX2013Z005), and the Scientific Research Fund of Sichuan Provincial Education Department (15ZA0288). * Corresponding author. Fax: +86 28 61831280. E-mail addresses: [email protected] (W.-H. Luo), [email protected] (T.-Z. Huang), [email protected] (X.-M. Gu), [email protected] (Y. Liu).

http://dx.doi.org/10.1016/j.aml.2016.12.011 0893-9659/© 2016 Elsevier Ltd. All rights reserved.

W.-H. Luo et al. / Applied Mathematics Letters 68 (2017) 13–19

14

Cirillo et al. [3] extended these results and proved that (1) and their derivatives converge at the rate of O(hd+1−k ), where hj is the local mesh size, for any k ≥ 0 and any set of well-spaced nodes. For (1), several j particular cases may be observed: when d = n, (1) becomes the usual Lagrange interpolant; taking d = 0, one gets the first linear rational interpolant [4]; a second interpolant of [4], which reduces to d = 1 for equispaced points, displays exponential convergence for analytic functions sampled at transformed Chebyshev points [5]. Due to its high accuracy, the interpolant (1), as a tool for approximating smooth functions, is becoming increasingly attractive. In fact, before the general result (1), Baltensperger et al. [6] have presented a collocation method which can be seen as the case d = 0 in (1). Berrut et al. [7] constructed a linear rational pseudospectral method for a two-point boundary value problem. For a fourth-order problem arising in thin film flows, Cueto-Felgueroso et al. [8] constructed an adaptive rational spectral method. With general d > 0, Klein [9] and Sharp [10] studied various spectral methods for some boundary value problems (BVPs), coupled system of singularly perturbed linear/nonlinear BVPs. Recently, G¨ uttel et al. [11] presented suggestions for choosing d for stable convergence. The aim of this paper is to use (1) to establish collocation methods for a class of parabolic PDEs ∂2u ∂u ∂u = k2 2 + k3 + f (u, x, t), t ∈ (0, T ], x ∈ (a, b); ∂t ∂x ∂x u(x, 0) = ϕ(x), x ∈ [a, b];

k1

(2)

u(a, t) = ξ(t), u(b, t) = ζ(t), t ∈ (0, T ]. Various methods have been introduced for some special cases of Eq. (2). Taking k3 = 0, f (u, x, t) = −u(1 − u)(ρ − u), Dehghana et al. [12] studied a homotopy perturbation method (HPM), a variational iteration method and an Adomian decomposition method for the Fitzhugh–Nagumo equation. Also, for the Fitzhugh–Nagumo equation with variable coefficients, Bhrawy [13] presented a Jacobi–Gauss–Lobatto collocation method. In the case of k3 = 0, f (u, x, t) = au − buq , (2) is called the Newell–Whitehead–Segel equation, for which Maci´ as-Di´ az et al. [14] proposed a finite-difference scheme. By means of the homotopy analysis, Sayevand et al. [15] studied a nonlinear equation which can be seen as an extension of Eq. (2). Jafari et al. [16] presented a fourth-order formula by using HPM. For k3 = 0, f (u, x, t) = ρu(1 − u) which corresponds to the known Fisher’s equation, Zhao [17] studied a discrete singular convolution algorithm. In this paper, using the interpolation (1), we construct three collocation methods for (2) with general f (u, x, t), give corresponding algebraic equations, and show the efficiency of these schemes by some examples. Compared with some existing results, these methods are especially well designed for the solution of PDEs sampled at equispaced points. The remainder of the paper is arranged as follows: In Section 2, some preliminary results are first given. Then, details about three collocation methods are given by using the R–K method and the barycentric rational quadrature, corresponding algebraic equations are formed. In Section 3, practical models are taken to numerically test the presented methods. Finally, some conclusions about these schemes are drawn. 2. Description of the numerical methods 2.1. Preliminary results Given a two-dimensional function u(x, t) : ([a, b] × [0, T ]) → R and two uniform partitions a = x1 < · · · < xM +1 = b, 0 = t1 < · · · < tN +1 = T. For some positive integers dx ≤ M, dt ≤ N , let Jm = {i ∈ Z, m − d ≤ ρm

λn

t−tn T m ∑Mx−x i ≤ m}, Jn = {i ∈ Z, n − d ≤ i ≤ n}, h = b−a +1 ρk , ψn (t) = ∑N +1 λk , M , τ = N , define φm (x) = k=1 x−xk k=1 t−tk ∑ ∏i+d ∑ ∏i+d 1 with ρm = i∈Jm (−1)i j=i,j̸=m xm1−x , λn = i∈Jn (−1)i j=i,j̸=n tn −t (see [1]). j

j

W.-H. Luo et al. / Applied Mathematics Letters 68 (2017) 13–19

15

For conveniently describing the coming collocation methods, we denote by ⎞ ⎛ ⎞ ⎛ ′ (i) (i) ′ φ1 (x1 ) ··· φM +1 (x1 ) ψ1 (t1 ) ··· ψN +1 (t1 ) ⎟ ⎜ ⎜ ⎟ .. .. .. .. ⎟ , i = 1, 2. W1 = ⎝ ⎠ , Vi = ⎜ . . . . ⎠ ⎝ ′ ′ (i) (i) ψ1 (tN +1 ) · · · ψN +1 (tN +1 ) φ1 (xM +1 ) · · · φM +1 (xM +1 ) ⎛ 2 2 ⎞ ω1 · · · ωq+1 ∫ tm q+1 ∑ λj λl ⎜ .. .. ⎟ m W2 = ⎝ . / ds, j = 1, . . . , q + 1, m = 2, . . . , q + 1. . ⎠ , ωj = s − t s − tl j t1 q+1 q+1 l=1 · · · ωq+1 ω1 Next, we will describe three collocation methods using (1). Without loss of generality, we take k1 = 1. 2.2. The Runge–Kutta-barycentric rational collocation method(RK-C method) ∑M +1 If we approximate solution of (2) by u(x, t) = j=1 φj (x)uj (t), then, considering the boundary conditions, the resulting ODEs corresponding to the points xi , i = 1, . . . , M + 1 are du = Hu + f , dt

(3)

where u = [u1 (t), . . . , uM +1 (t)]T , u1 (t) = ξ(t), uM +1 (t) = ζ(t), H(1, 1 : M + 1) = 0, H(M + 1, 1 : M + 1) = 0, H(2 : M, 1 : M + 1) = k2 V2 (2 : M, 1 : M + 1) + k3 V1 (2 : M, 1 : M + 1), and f = [ξ ′ (t), f (u2 (t), x2 , t), . . . , f (uM (t), xM , t), ζ ′ (t)]T . Denote by ˜f (t, u) the term Hu + f . Since (3) is a system of ODEs, we can use the classical four-order Runge–Kutta formula to iteratively solve it. 2.3. The barycentric rational quadrature-collocation method (B-Q-C method) 2

Let g(u, x, t) = k2 ∂∂xu2 + k3 ∂u ∂x + f (u, x, t), then, the first equation in (2) can be rewritten as ∂u = g(u, x, t), t ∈ (0, T ], x ∈ (a, b). ∂t

(4)

For fast computing (4), we take a fixed q(q = kN ), and divide [t1 , tN +1 ] into kf ix subintervals f ix [t(i−1)q+1 , tiq+1 ], i = 1, . . . , kf ix . In every subintervals [t(i−1)q+1 , tiq+1 ], integrating (4) from t(i−1)q+1 to t, there is ∫ t u(x, t) = u(x, t(i−1)q+1 ) + g(u, x, s)ds (5) t(i−1)q+1

∑iq+1

let u(x, t) = m=(i−1)q+1 ψm (t)um , employing the barycentric rational quadrature (see [18,19]), then, the ∑M +1 use of interpolation um = s=1 φs (x)us,m , m = (i − 1)q + 1, . . . , iq + 1 and boundary values will result in q(M + 1) algebraic equations [I1 − k2 W2 (:, 2 : q + 1) ⊗ V2 − k3 W2 (:, 2 : q + 1) ⊗ V1 ]u − W2 (:, 2 : q + 1) ⊗ I2 f1 = S1 ⊗ u(i−1)q+1 + W2 (:, 1) ⊗ [(k2 V2 + k3 V1 )u(i−1)q+1 + f2 ], where I1 , I2 are respectively q(M + 1) and M + 1-dimensional identity matrices, and u = [u1,(i−1)q+2 , . . . , uM +1,(i−1)q+2 , . . . , u1,iq+1 , . . . , uM +1,iq+1 ]T , u(i−1)q+1 = [u1,(i−1)q+1 , . . . , uM +1,(i−1)q+1 ]T , f1 = [fˆ(i−1)q+2 , . . . , fˆiq+1 ]T , fˆk = [f (u1,k , x1 , tk ), . . . , f (uM +1,k , xM +1 , tk )], k = (i − 1)q + 2, . . . , iq + 1, f2 = [f (u1,(i−1)q+1 , x1 , t(i−1)q+1 ), . . . , f (uM +1,(i−1)q+1 , xM +1 , t(i−1)q+1 )]T , S1 = [1, . . . , 1]Tq×q .

(6)

W.-H. Luo et al. / Applied Mathematics Letters 68 (2017) 13–19

16

Some characteristics about this scheme can be summarized as follows: • by the properties of barycentric rational weights λq , the quadrature weights matrix W2 can be repeatedly used in all the subintervals [t(i−1)q+1 , tiq+1 ], i = 1, . . . , kf ix ; • the initial value u(x, t(i−1)q+1 ), i = 2, . . . , kf ix can always be obtained from the previous step.

2.4. The barycentric rational collocation method (B-C method) Let u(x, t) =

∑M +1 ∑N +1 m=1

n=1

φm (x)ψn (t)um,n . Using the initial condition and boundary conditions, we get

ui,1 = ϕ(xi ), i = 1, . . . , M + 1, u1,j = ξ(tj ), uM +1,j = ζ(tj ), j = 2, . . . , N + 1. For solving the remaining unknowns ui,j , i = 2, . . . , M, j = 2, . . . , N + 1, we rewrite u(x, t) as u(x, t) =

M N +1 ∑ ∑

φm (x)ψn (t)um,n +

m=2 n=2

N +1 ∑

ψn (t)(φ1 (x)u1,n + φM +1 (x)uM +1,n ) +

n=2

M +1 ∑

φm (x)ψ1 (t)um,1 ,

m=1

and let U = [u2,2 , . . . , u2,N +1 , u3,2 , . . . , u3,N +1 , . . . , uM,2 , · · · , uM,N +1 ]T , U1 = [u1,2 , . . . , u1,N +1 ]T , U2 = [uM +1,2 , . . . , uM +1,N +1 ]T , U3 = [u1,1 , . . . , uM +1,1 ]T , F = [f1 , . . . , fM −1 ]T , where fi = [f (ui+1,2 , xi+1 , t2 ), . . . , f (ui+1,N +1 , xi+1 , tN +1 )], i = 1, . . . , M − 1. Substituting (xi , tj ), i = 2, . . . , M, j = 2, . . . , N + 1 into the first equation in (2), we obtain the collocation equations RU + F + G = 0

(7)

where P1 = IM −1,M −1 , P2 = V2 (2 : M, 2 : M ), P3 = V1 (2 : M, 2 : M ), Q1 = W1 (2 : N + 1, 2 : N + 1), Q2 = IN,N , Q3 = IN,N , A1 = IM +1,M +1 (2 : M, 1 : M + 1), A2 = V2 (2 : M, 1), A3 = V1 (2 : M, 1), B1 = W1 (2 : N + 1, 1), B2 = IN,N , B3 = IN,N , C1 = V2 (2 : M, M + 1), C2 = V1 (2 : M, M + 1), D1 = B2 , D2 = B3 , R = −k1 (P1 ⊗ Q1 ) + k2 (P2 ⊗ Q2 ) + k3 (P3 ⊗ Q3 ), G = −k1 (A1 ⊗ B1 )U3 + k2 [(A2 ⊗ B2 )U1 + (C1 ⊗ D1 )U2 ] + k3 [(A3 ⊗ B3 )U1 + (C2 ⊗ D2 )U2 ]. 3. Numerical examples In this section, we give some numerical examples to illustrate the efficiency of the above three schemes. For the RK-C method, we use “ODE45” from Matlab r2010b to solve (3). For examining the convergence in space, we set both ‘RelTol’ and ‘AbsTol’ as 10−10 . The total number of iterations Tstep is automatically determined by the Matlab software. For the B-Q-C method (Eq. (6)), by the theoretical results in [1,2,18], the accuracy in space and time are expected to be O(hdx −1 ) and O(τ dt +1 ), hence, we compute the errors at all the points (xi , tj ), i = 1, . . . , M +1, j = 1, . . . , N +1 by taking h = τ = N1 and dx = dt +2. Similarly, by [1,2], the parameters in the B-C method (Eq. (7)) are taken as dx = dt + 1. In the listed tables, we compute the errors using L∞ norm, and denote by “En ” the errors at all the collocation points. The notations “Ordern ”, n (N ) computed by log(En (N/2))/E , show the convergence rate of the errors at the collocation points. log2 In order to check the convergence rate corresponding to various dx , dt , we take different values among these examples.

W.-H. Luo et al. / Applied Mathematics Letters 68 (2017) 13–19

17

Table 1 Results for Example 1. RK-C (dx = 6) M

En

Ordern

Tstep

Time (s)

16 32 64 128

7.9832e−7 2.1474e−8 5.8848e−10 1.5902e−11

5.2163 5.1895 5.2097

2 11 47 192

5.3242 34.4641 269.3531 2.1675e+3

B-Q-C (q = 8, dt = 4, dx = 6, M = N )

533 353 269 209

B-C method (dt = 5, dx = 6, M = N )

M

En

Ordern

Time (s)

En

Ordern

Time (s)

16 32 64 128

8.1913e−7 1.9421e−8 4.8682e−10 1.3280e−11

5.3984 5.3181 5.1961

0.3840 0.4638 1.3319 9.1236

8.1702e−7 1.9434e−8 4.9592e−10 1.2767e−11

5.3937 5.2923 5.2796

0.0339 0.2412 5.6792 170.8573

Table 2 Results for Example 2. RK-C method (dx = 8) M

En

16 32 64 128

2.1309e−4 5.4498e−6 4.8183e−8 3.8593e−10

Ordern

Tstep

Time (s)

5.2891 6.8215 6.9640

6 15 47 192

14.2155 51.2954 312.3647 2.3944e+3

373 409 965 073

B-Q-C method (q = 8, dt = 6, dx = 8, M = N )

B-C method (dt = 7, dx = 8, M = N )

M

En

16 32 64 128

5.2212e−3 1.4050e−4 1.1272e−6 8.8261e−9

Ordern

Time (s)

En

Ordern

Time (s)

5.2157 6.9617 6.9968

0.1222 0.1768 0.8546 7.5962

5.5418e−3 1.0410e−4 8.1921e−7 6.4257e−9

5.7343 6.9895 6.9942

0.0305 0.1754 3.6840 285.3472

Example 1. The Fitzhugh–Nagumo equation [12]: √ exact solution u = Example 2.

∂u ∂t

1 2

=

+ ∂2u ∂x2

1 2

tanh(

2x+(1−2α)t ), 4

∂u ∂t

=

∂2u ∂x2

+ u(u − α)(1 − u), t ∈ (0, 1], x ∈ (0, 1), with

here we choose α = 31 .

+ u3 + f (x, t), t ∈ (0, 1], x ∈ (0, 1), with exact solution u =

Example 3. The heat equation

∂u ∂t

=

∂2u ,t ∂x2

1 . (1+25x2 )(1+25t2 )

∈ (0, 1], x ∈ (0, 1), with exact solution u = exp(x + t).

From the results listed in these tables we can make the following observations: 1. For any fixed dh , dt , the RK-C method achieves an accuracy of O(hdx −1 ) in space, meanwhile, the experimental convergence rates of the B-Q-C method and the B-C method respectively attain O(hdx −1 + τ dt +1 ) and O(hdx −1 + τ dt ) in time and space. Compared with some other classical schemes, such as finite difference methods [14] (accuracy of O(τ + h2 )), HPM [16] (accuracy of O(τ 4 + h4 )), these presented methods possess higher accuracy in space and time. 2. From Table 2, we can see that, for appropriately large d, the RK-C method, B-Q-C method and B-C method can all well approximate the solution which represents a Runge function. 3. The rate of convergence of the RK-C method, B-Q-C method and B-C method is stable as the meshes are refined. 4. Concerning CPU time, Tables 1–3 show that for small enough steps in space, the B-Q-C scheme is the absolute winner comparing with both the B-C scheme and the RK-C method. 5. For the RK-C method, some remarkable features must be observed. With the refinement of the grids in space, the meshes in time and the involved time increase sharply. This can be attributed to the requirements for the stability of the RK-C method. In fact, for the ODEs (3), it is easy to obtain the stability domain

W.-H. Luo et al. / Applied Mathematics Letters 68 (2017) 13–19

18

Table 3 Results for Example 3. RK-C method (dx = 5) M

En

Ordern

Tstep

Time (s)

16 32 64 128

5.4376e−5 1.5014e−6 5.1210e−8 1.5475e−9

5.1786 4.8737 5.0484

5 12 46 192

9.3527 31.6875 270.3672 2.0157e+3

004 768 735 210

B-Q-C method (q = 8, dt = 3, dx = 5, M = N )

B-C method (dt = 4, dx = 5, M = N )

M

En

Ordern

Time (s)

En

Ordern

Time (s)

16 32 64 128

5.0014e−5 2.5504e−6 1.3435e−7 7.1526e−9

4.2935 4.2467 4.2314

0.3300 0.3394 0.4161 1.3998

6.0132e−5 3.2454e−6 1.6972e−7 9.0457e−9

4.2117 4.2572 4.2298

0.0343 0.0998 1.6237 81.0432

1 64 ,

1 128 ).

Fig. 1. Eigenvalue distributions of H in the Example 3 (dx = 5, left: h =

right: h =

1 ∥I + τ H + 12 (τ H)2 + 16 (τ H)3 + 24 (τ H)4 ∥ ≤ 1. This means bigger ∥H∥ will lead to smaller time steps. Fig. 1 shows that all the eigenvalues of the matrix H have negative real parts, and the spectral radius σ increases sharply. Hence, the rapid increase of the time steps Tstep can be well explained by using the usual inequity ∥H∥m ≥ σ, where ∥H∥m is any consistent matrix norm of H.

4. Conclusions Based on a family of barycentric rational interpolation functions, we have introduced three methods for a class of nonlinear parabolic PDEs. The efficiency of these schemes has been investigated by practical numerical examples, including the Fitzhugh–Nagumo equation which models the transmission of nerve impulses. References [1] M.S. Floater, K. Hormann, Barycentric rational interpolation with no poles and high rates of approximation, Numer. Math. 107 (2007) 315–331. [2] J.-P. Berrut, M.S. Floater, G. Klein, Convergence rates of derivatives of a family of barycentric rational interpolants, Appl. Numer. Math. 61 (2011) 989–1000. [3] E. Cirillo, K. Hormann, J. Sidon, Convergence rates of derivatives of Floater-Hormann interpolants for well-spaced nodes, Appl. Numer. Math. (2016) (http://dx.doi.org/10.1016/j.apnum.2016.07.008). [4] J.-P. Berrut, Rational functions for guaranteed and experimentally well-conditioned global interpolation, Comput. Math. Appl. 15 (1988) 1–16.

19

W.-H. Luo et al. / Applied Mathematics Letters 68 (2017) 13–19

[5] R. Baltensperger, J.-P. Berrut, B. No¨ el, Exponential convergence of a linear rational interpolant between transformed Chebyshev points, Math. Comp. 68 (1999) 1109–1120. [6] R. Baltensperger, J.-P. Berrut, The linear rational collocation method, J. Comput. Appl. Math. 134 (2001) 243–258. [7] J.-P. Berrut, R. Baltensperger, The linear rational pseudospectral method for boundary value problems, BIT 41 (2001) 868–879. [8] L. Cueto-Felgueroso, R. Juanes, Adaptive rational spectral methods for the linear stability analysis of nonlinear fourth-order problems, J. Comput. Phys. 228 (2009) 6536–6552. [9] G. Klein, Une extension de la m´ ethode de collocation rationnelle lin´ eaire pour la r´ esolution de probl` emes de valeurs aux limites (Masters thesis), Department of Mathematics, University of Fribourg, Switzerland, 2008. [10] N. Sharp, Barycentric rational interpolation and collocation methods (Masters thesis), Simon Fraser University, 2014. [11] S. G¨ uttel, G. Klein, Convergence of linear barycentric rational interpolation for analytic functions, SIAM J. Numer. Anal. 50 (2012) 2560–2580. [12] M. Dehghana, J.M. Heris, A. Saadatmandi, Application of semi-analytic methods for the Fitzhugh–Nagumo equation, which models the transmission of nerve impulses, Math. Methods Appl. 33 (2010) 1384–1398. [13] A.H. Bhrawy, A Jacobi–Gauss–Lobatto collocation method for solving generalized Fitzhugh–Nagumo equation with time-dependent coefficients, Appl. Math. Comput. 222 (2013) 255–264. [14] J.E. Maci´ as-Di´ az, J. Ruiz-Ram´ırez, A non-standard symmetry-preserving method to compute bounded solutions of a generalized Newell–Whitehead–Segel equation, Appl. Numer. Math. 61 (2011) 630–640. [15] K. Sayevand, D. Baleanu, M. Fardi, Travelling wave solutions: A new approach to the analysis of nonlinear physical phenomena. Centr, Eur. J. Phys. 12 (7) (2014) 480–489. [16] K. Sayevand, H. Jafari, On systems of nonlinear equations: some modified iteration formulas by homotopy perturbation method with accelerated fourth- and fifth-order convergence, Appl. Math. Model. 40 (2) (2015) 1467–1476. [17] S. Zhao, G.W. Wei, Comparison of the discrete singular convolution and three other numerical schemes for solving Fisher’s equation, SIAM J. Sci.Comput. 25 (2003) 127–147. [18] J.P. Berrut, S.A. Hosseini, G. Klein, The linear barycentric rational quadrature method for Volterra integral equations, SIAM J. Sci. Comput. 36 (36) (2014) A105–A123. [19] G. Klein, J.P. Berrut, Linear barycentric rational quadrature, BIT 52 (2) (2012) 407–424.