An algorithm based on a new DQM with modified exponential cubic B-splines for solving hyperbolic telegraph equation in (2 + 1) dimension

1 Introduction

The hyperbolic partial differential equations have a great attention due to its useful understanding in various physical phenomena of applied sciences and engineering, for instance, hyperbolic partial differential equation models fundamental equations in atomic physics [1], vibrations of structures (e.g. buildings, machines and beams).

Consider second order two-space dimensional linear hyperbolic telegraph equation

together with the initial conditions

and the boundary conditions either of Dirichlet type (3) or Neumann type (4) as follows

where ψ, ϕ, ψ_i, ϕ_i (i = 1,2,3,4) are known smooth functions.

where ∂Ω is the boundary of the computational domain Ω = [0, 1]², α > 0, β are arbitrary constants. Eq. (1) with β = 0 is a damped wave equation while for β > 0 it becomes telegraph equation which is more convenient than diffusion equation in modeling reaction-diffusion for such branches of sciences [2], and mostly used in wave propagation of electric signals in a cable transmission line [3], etc.

Bellman et al. [15] developed “differential quadrature method (DQM)” for numeric study of partial differential equations (PDEs). After seminal work of Bellman et al., Quan and Chang [16, 17], the DQM has been employed with different base functions such as polynomial based DQM (PDQM) [7], cubic B-spline DQM [18, 19], MCB-DQM [9, 20, 22, 23], DQM based on fourier expansion and Harmonic function [24], sinc DQM [25], generalized DQM [26], quartic B-spline based DQM [27, 28], Quartic and quintic B-spline methods [29, 30] and exponential cubic B-spline DQM [31].

In the past years, various techniques have been developed for numerical study of PDEs [6, 7, 12]. 2D telegraph equation has been studied numerically via. Taylor matrix method [5], meshless local weak-strong (MLWS) method and meshless local Petrov-Galerkin (MLPG) method [2], higher order implicit collocation method [6], PDQM [7, 8], modified cubic B-spline DQM (MCB-DQM) [9], an alternating direction implicit scheme [10], hybrid method by Dehghan and Salehi [11], modified extended cubic B-spline DQM (mECDQ) [12], compact finite difference scheme with accuracy of order four in both space and time [13] and a meshless scheme with radial basis functions [14].

In this paper our aim is to proposed modified exponential cubic-B-spline differential quadrature method (mExp-DQM) for hyperbolic PDEs. The mExp-DQM is used to convert the initial- boundary value system of the telegraph equation into a initial value system of ODEs, in time. Among various time integration algorithm SSP-RK54 algorithm [32] is adopted to solve the resulting system of ODEs.

The rest of the paper is organized into five more sections. Section 2 deals with the description of mExp-DQM. Section 3 presents the procedure for the implementation of mExp-DQM for the problem (1) with the initial conditions (2) and the boundary conditions (3) or (4). The stability analysis of mExp-DQM is studied in Section 4. The main goal, the numeric study of six test problems to establish the accuracy of the proposed method is carried out in terms of the relative error (R_e), L₂ and L_∞ error norms in Section 5. Finally, Section 6 concludes the paper with reference to critical analysis and research perspectives.

2 Description of mExp-DQM

The DQM is an approximation to derivatives of a function is the weighted linear sum of the functional values at certain grid points[15], where the weighting coefficients depend on grids only [26], and so, a uniform partition P[Ω] = {(x_i, y_j) ∈ Ω: h_x = x_i+1 – x_i, h_y = y_j+1 – y_j, i ∈ Δ_x, j ∈ Δ_y}, of the domain Ω = {(x, y) ∈ R²: 0 ≤ x, y ≤ 1} of the problem is distributed with the following grid points:

where Δ_x = {1, 2, …, N_x}, Δ_y = {1, 2, …, N_y}, and h_x = 1Nx−1 and h_y = 1Ny−1 are the discretization steps in both x and y directions, respectively.

Let (x_i, y_j) be the generic grid point and

The approximation for r-th order derivative of u(x, y, t), for r ∈ {1, 2}, with respect to x, y at (x_i, y_j) for i ∈ Δ_x, j ∈ Δ_y is given by

where time dependent unknown quantities aiℓ(r)andbjℓ(r) are called weighting functions of the rth-order derivative, which are computed using a set of base functions.

The exponential cubic B-splines function ζ_i = ζ_i(x) at node i in x direction, reads [31, 33]:

where

The set {ζ₀, ζ₁, ζ₂, …, ζ_Nx, ζ_Nx+1} forms a basis over the interval [a, b]. The values of ζ_i and its first and second derivatives in the grid point x_j, denoted by ζ_ij := ζ_i(x_j), ζij′:=ζi′(xj)andζij″:=ζi″(xj), respectively, read:

Analogous to [20, 23], modified exponential cubic B-splines base functions are obtained by modifying exponential cubic B-splines (6) as follows:

The set {ψ₁, ψ₂, …, ψ_Nx} forms a basis over [a, b] in direction of x. Analogously procedure is followed for y direction.

3 The mExp-DQM for the telegraph equation

First, we set u_t = v and thus u_tt = v_t. Keeping all above in mind, telegraph equation (1) with the initial condition reduces to

Set f(x_i, y_j, t) = f_ij. mExp-DQM converts Equation (12) to

No further simplification is required whenever the boundary conditions are of Dirichlet type. The solution on the boundaries can be read directly from the conditions (3) as:

On the other hand if the boundary conditions are of Neumann or mixed type, further simplification is required. The procedure to find solutions at the boundary is discussed in [9, 34], in detail. The mExp-DQM on the boundaries yields a system of linear equations. The solutions at the boundaries are obtained by solving these resulting system of linear equations as follows:

The solution for u_1j, u_Nxj is obtained from Eq. (5) with r = 1 and conditions (4) at x = 0 and x = 1 as

where S1(j)=ψ1(j)−∑ℓ=2Nx−1a1ℓ(1)uℓjandS2(j)=ψ2(j)−∑ℓ=2Nx−1aNxℓ(1)uℓj and the solutions for the boundary values u_i1 and u_iNy are obtained from Eq. (5) with r = 1 and the conditions (4) at y = 0 and y = 1 as

where S3(i)=ψ3(i)−∑ℓ=2Ny−1b1ℓ(1)uiℓandS4(i)=ψ4(i)−∑ℓ=2Ny−1bNyℓ(1)uiℓ.

Once the solutions u_1j, u_Nxj, u_i1 and u_iNy on boundary are obtained from either boundary condition (Dirichlet type (3) or Neumann type (4)), Eq. (13) can be rewritten as follows:

where 2 ≤ i ≤ N_x – 1, 2 ≤ j ≤ N_y – 1 and

The initial valued system (17) can be solved by using various existing time integration schemes. The SSP-RK algorithm allows low storage and large region of absolute property [32]. In what follows, SSP-RK54 algorithm, strongly stable for nonlinear hyperbolic differential equations, has been adopted.

4 Stability analysis

In compact form, the system (17) can be rewritten as follows:

where

1. A=OIB−2αI,G=O1K,U=uv,andU0=ϕψ

2. O and O₁ are null matrices;

3. I is the identity matrix of order (N_x – 2)(N_y – 2);

4. U = (u, v)^T: u = (u₂₂, u₂₃, …, u_{2(Ny – 1)}, u₃₂, u₃₃, …, u_{3(Ny – 1)}, …, u_{(Nx – 1)2}, …, u_{(Nx – 1)(Ny – 1)}), v = (v₂₂, v₂₃, …, v_{2(Ny – 1)}, v₃₂, v₃₃, …, v_{3(Ny – 1)}, …, v_{(Nx – 1)2}, …, v_{(Nx – 1)(Ny – 1)}).

5. K = (K₂₂, K₂₃, …, K_{2(Ny – 1)}, K₃₂, …, K_{3(Ny – 1)}, … K_{(Nx – 1) 2}, … K_{(Nx – 1) (Ny – 1)}, where K_ij, for i ∈ Δ_x, j ∈ Δ_y is calculated from Eq. (18).

6. B = –β²I + B_x + B_y, where B_x and B_y are the following matrices (of order (N_x – 2)(N_y – 2)) of the weighting coefficients aij(2)andbij(2):

   Bx=a22(2)Ixa23(2)Ix…a2(Nx−1)(2)Ixa32(2)Ixa33(2)Ix…a3(Nx−1)(2)Ix⋮⋮⋱⋮a(Nx−1)2(2)Ixa(Nx−2)3(2)Ix…a(Nx−1)(Nx−1)(2)Ix,By=MyOy…OyOyMy…Oy⋮⋮⋱⋮OyOy…My,(20)

   where identity matrix, I_x, and null matrix, O_y, both are of order (N_y – 2) and

   My=b22(2)b23(2)…b2(Ny−1)(2)b32(2)b33(2)…b3(Ny−1)(2)⋮⋮⋱⋮b(Ny−1)2(2)b(Ny−1)3(2)…b(Ny−1)(Ny−1)(2).

It is worth mentioning that the stability of mExp-DQM for equation (1) depends upon the stability of system (19) and the proposed method for temporal discretization may not converge to the exact solution whenever system (19) is unstable. The stability of system (19) depends upon the eigenvalues of the coefficient matrix A [35]. In fact, the stability region is the set 𝓢 = {z ∈ C :| R(z)| ≤ 1, z = λ_A △t}, where R(.) is the stability function, λ_A →eigenvalue of the coefficient matrix A. The stability region 𝓢 for SSP-RK54 algorithm is depicted in Fig 1, see [36, Fig. 5], which confirm that the sufficient condition for the stability of system (19) is that λ_A △t ∈ 𝓢 for all λ_A, i.e., the real part of each eigenvalue is necessarily either zero or negative.

Fig. 1

Stability region for SSP-RK54 algorithm (left) and eigenvalues of B_x + B_y for different grid sizes h = 0.1, 0.01, 0.025, 0.016

Let λ_A be an eigenvalue associated with eigenvector (X₁, X₂)^T, where each component is a vector of order (N_x – 2)(N_y – 2), then Eq. (19) can be re-written as

which yields

Eq. (22) implies that eigenvalue λ_B of B is given by λ_B = λ_A (λ_A + 2α), where B = –β²I + B_x + B_y.

The sufficient condition for the stability of system (19) is that the real part Re(λ_A) for each eigenvalue λ_A of A is either zero or negative, i.e., Re(λ_A) ≤ 0. Fig. 1 shows that computed eigenvalue λ of B_x + B_y for p = 1 with different grid sizes (h_x = h_y = h = 0.1, 0.05, 0.025, 0.016) are real and negative, and so, each eigenvalue λ_B = –β² + λ of B is real and negative, i.e.,

where Re(z) = real part of z, Im(z) = imaginary part of z.

Setting λ_A = x+ι y, and so, λ_B = x²–y² + 2αx + 2ι(x+α)y, then Eq. (23) yields

The above equations are hold simultaneously only when either x = –α or (x + α)² < α² implying that Re(λ_A) < 0 for each α > 0.

5 Numerical experiments and discussion

This section deals with the main goal of the paper, the computation of numerical solutions of 2D telegraph equation using mExp-DQM with SSP-RK54 algorithm. The accuracy and the efficiency of this method is measured in terms of various error norms (L₂, relative error (R_e),L_∞) for six examples.

where uj∗ represent the numerical solution at node j. Throughout this section, equal grid size is considered in each direction, i.e., h_x = h_y = h.

Example 5.1

Consider telegraph equation(1)with f(x, y, t) = (–2α + β² – 1)exp(–t)sinh x sinh y, ϕ(x, y) = sinh x sinh y, ψ(x, y) = –sinh x sinh y in Ω;ϕ₁(y, t) = 0; ϕ₂(y, t) = exp(–t)sinh(1)sinh y for 0 ≤ y ≤ 1 andψ₃(x, t) = 0; ψ₄(x, t) = exp(–t)sinh x sinh(1) for 0 ≤ x ≤ 1. The exact solution [7] is given by

u(x,y,t)=exp(−t)sinh⁡xsinh⁡y.(25)

The solutions are computed for α = 10, β = 5, α = 10, β = 0 with parameters △t = 0.01, 0.001, h = 0.1, 0.05 and p = 1. The L₂, L_∞errors and CPU time for different time levels t ≤ 10 are compared with the error norms by Mittal and Bhatia [9] taking △t = 0.01, h = 0.1 and are reported in Table 2. In Table 3, the computed results are compared with Mittal and Bhatia [9] and Jiwari et al. [7] for △t = 0.001 and h = 0.05. The findings from the above tables confirms that the proposed results are better than [7, 9]. The CPU time is slightly more than [9] due to selection of SSP-RK54 algorithm instead of SSP-RK43 algorithm, for time integration. The surface plots of numerical at t = 1, 2, 3 with △t = 0.001 and h = .05 are depicted in Fig. 2.

Table 2

Comparison of the mExp-DQM solutions of Example 5.1 with α = 10, β = 5, △t = 0.01, h = 0.1 and p = 1

t	MCB-DQM [9]				mExp-DQM
	L2	L∞	Re	CPU(s)	L2	L∞	Re	CPU(s)
0.5	8.3931E-04	3.3019E-03	2.8902E-03	0.13	6.8998E-06	1.0168E-05	2.8749E-04	0.016
1	6.0254E-04	2.0597E-03	3.4208E-03	0.16	5.3522E-06	7.0133E-06	3.6767E-04	0.046
2	2.4167E-04	7.6531E-04	3.7297E-03	0.19	2.2337E-06	2.8534E-06	4.1711E-04	0.078
3	8.9534E-05	2.7920E-04	3.7937E-03	0.24	8.3375E-07	1.0585E-06	4.2747E-04	0.141
5	1.2168E-05	3.7800E-05	3.8097E-03	0.34	1.1352E-07	1.4389E-07	4.3005E-04	0.218

Table 3

Comparison of mExp-DQM solutions of Example 5.1 with Δt = 0.001, α = 10, β = 0,5, p = 1 and h = 0.05

t	mExp-DQM				MCB-DQM [9]				PDQM [7]
ß = 5	L2	L∞	Re	CPU(s)	L2	L∞	Re	CPU(s)	Re	CPU(s)
0.5	8.1273E-07	1.3152E-06	6.6847E-05	1.279	1.0690E-04	2.4738E-04	1.1088E-04	0.47	1.1185E-04	6
1	5.8429E-07	8.3976E-07	7.9233E-05	2.496	1.5293E-05	3.3082E-04	1.3266E-04	1.10	1.8051E-04	12
2	2.3507E-07	3.2200E-07	8.6737E-05	5.896	4.6468E-05	1.1380E-05	3.1954E-04	1.10	4.7289E-04	25
3	8.8032E-08	1.1937E-07	8.8297E-05	7.488	2.1994E-05	4.3577E-05	1.3024E-04	2.80	1.2656E-04	37
5	1.1979E-08	1.6202E-08	8.8694E-05	12.402	2.7151E-06	5.4141E-06	1.4439E-04	4.30	9.2770E-04	62
ß = 0
0.5	7.2898E-07	9.0815E-07	5.9958E-05	1.279	9.2959E-05	4.2348E-04	3.4675E-04	0.52	1.1198E-04	6
1	8.0739E-07	1.0270E-06	1.0949E-04	2.511	6.3652E-05	2.5838E-04	3.9146E-04	0.98	1.8635E-04	12
2	5.7525E-07	7.2622E-07	2.1226E-04	5.007	2.5540E-05	9.5843E-05	4.2739E-04	1.80	5.1797E-04	25
3	3.1155E-07	3.9340E-07	3.1248E-04	7.394	9.9234E-06	3.5340E-05	4.5140E-04	2.20	1.4412E-04	37
5	6.7799E-08	8.5767E-08	5.0198E-04	12.470	1.5116E-06	4.8043E-06	5.0758E-04	4.50	1.0883E-04	62

Fig. 2

Plots of numerical solution at different time levels for Example 5.1

Example 5.2

Consider telegraph equation(1)in the region Ω with α = β = 1, f(x, y, t) = –2 exp(x + y – t), ϕ(x, y) = exp(x + y),ψ(x, y) = –exp(x + y) in Ω and the mixed boundary conditionsϕ₁(y, t) = exp(y – t), ϕ₂(y, t) = exp(1 + y – t) for 0 ≤ y ≤ 1 andψ₃(x, t) = exp(x – t), ϕ₄(x, t) = exp(1 + x – t) for 0 ≤ x ≤ 1. The exact solution [4] is given by

u(x,y,t)=exp(x+y−t).(26)

The computed results and CPU time are compared with the results by Mittal and Bhatia [9] for different space step size h = 0.1, 0.05 and time step size △t = 0.01, 0.001, p = 1 and reported in Table 4 and Table 5. The surface plots of the mExp-DQM solutions at different time levels t = 1, 2, 4 is depicted in Fig. 3. The findings show that the proposed results are more accurate than the results by Mittal and Bhatia [9].

Table 4

Comparison of the mExp-DQM solutions of Example 5.2 with h = 0.1, Δ = 0.01, α = 1, β = 1, p = 1

t	mExp-DQM			MCB-DQM [9]
	L2	L∞	CPU(s)	L2	L∞	CPU(s)
1	3.9796E-04	6.7076E-04	0.031	1.4441E-02	2.9996E-02	0.03
2	4.5099E-05	1.1091E-04	0.063	1.3898E-03	3.9711E-03	0.05
3	4.0589E-05	7.4545E-05	0.109	1.3018E-03	2.2178E-03	0.08
5	4.1078E-06	8.6460E-06	0.187	1.1112E-04	2.0618E-04	0.11
7	4.6749E-07	1.0452E-06	0.234	1.3695E-05	3.0052E-05	0.14
10	3.8692E-08	7.0454E-08	0.312	1.4408E-06	2.5354E-06	0.19

Table 5

Comparison of the mExp-DQM solutions of Example 5.2 with h = .05, Δ = 0.001, α = β = 1, p = 1

t	mExp-DQM			MCB-DQM [9]
	L2	L∞	CPU(s)	L2	L∞	CPU(s)
0.5	1.28E-04	2.67E-04	1.045	3.4808E-03	9.5129E-03	0.50
1	1.05E-04	1.82E-04	2.074	3.2351E-03	7.4749E-03	0.70
2	1.04E-05	3.07E-05	4.181	2.8518E-04	1.0361E-03	1.30
3	1.09E-05	2.05E-05	6.255	3.1028E-04	5.7859E-04	1.90
5	1.03E-06	2.40E-06	10.358	2.4495E-05	6.7234E-05	3.30
7	1.03E-07	2.59E-07	14.446	2.5376E-06	8.2203E-06	3.90
10	1.18E-08	2.18E-08	20.498	3.6505E-06	8.5897E-06	5.20

Fig. 3

Plots of mExp-DQM solutions solutions at different time levels for Example 5.2

Example 5.3

Consider telegraph equation(1)in Ω with α = β = 1, f(x, y, t) = 2(cos t – sin t)sin x sin y, ϕ(x, y) = sin x. sin y; ψ(x, y) = 0 and the Dirichlet boundary conditions: ϕ₁(y, t) = 0, ϕ₂(y, t) = cos t sin(1)sin y, 0 ≤ y ≤ 1, ϕ₃(x, t) = 0, ϕ₄(x, t) = cos t sin x sin(1), 0 ≤ x ≤ 1.

The exact solution [4] is

u(x,y,t)=cos⁡tsin⁡xsin⁡y.(27)

The computed relative error (R_e), L₂, L_∞error norms are compared with the recent results of Mittal and Bhatia [9] at different time levels t ≤ 10, reported in Table 6 and Table 7 with the parameters △t = 0.01, h = 0.1, p = 1 and △t = 0.001 and h = 0.05, p = 0.15, 1, respectively. The physical solution behavior at t = 1, 2, 3 is depicted in Fig. 4. The findings shows that the proposed solution are much better than the results by Mittal and Bhatia [9], and are in excellent agreement with the exact solutions. The computation time is slightly more than Mittal and Bhatia [9] for large t due to selection of SSP-RK54 algorithm instead of SSP-RK43 algorithm in time integration.

Table 6

Comparison of the mExp-DQM solutions of Example 5.3 with Δt = 0.01 and h = 0.1, p = 1

t	mExp-DQM				MCB-DQM [9]
	L2	L∞	Re	CPU(s)	L2	L∞	Re	CPU(s)
1	3.7330E-06	4.5492E-06	2.7069E-04	0.031	9.9722E-04	2.2746E-03	5.9762E-03	0.08
2	4.4842E-06	5.6294E-06	4.2217E-04	0.062	1.0926E-03	2.8706E-03	8.5019E-03	0.11
3	3.7742E-06	6.3374E-06	1.4916E-04	0.109	2.2877E-04	6.0818E-04	7.4720E-04	0.14
5	4.4186E-06	5.3912E-06	5.9036E-04	0.203	1.1562E-03	2.9942E-03	1.2767E-03	0.20
7	3.2109E-06	3.7239E-06	1.6834E-04	0.312	7.2867E-04	1.8781E-03	3.1572E-03	0.26
10	3.1806E-06	3.7506E-06	1.4949E-04	0.374	5.8889E-04	1.5158E-03	2.2874E-03	0.34

Table 7

Comparison of the mExp-DQM solutions of Example 5.3 with Δt = 0.001; h = 0.05 and p = 1,0.15

t	mExp-DQM (p=1)				mExp-DQM (p=0.15)				MCB-DQM [9]
	L2	L∞	Re	CPU(s)	L2	L∞	Re	CPU(s)	L2	L∞	Re	CPU(s)
1	3.5715E-07	5.8718E-07	5.0162E-05	2.26	3.5729E-07	5.8736E-07	5.0182E-05	2.26	9.8870E-05	2.4964E-04	6.2977E-04	0.78
2	4.4969E-07	6.7211E-07	8.1823E-05	4.52	4.4977E-07	6.7225E-07	8.1838E-05	4.52	1.2148E-04	3.2296E-04	1.0025E-03	1.30
3	7.8128E-07	1.2228E-06	5.9879E-05	6.79	7.8151E-07	1.2231E-06	5.9897E-05	6.79	3.7627E-05	9.9310E-05	1.3078E-04	1.70
5	3.6743E-07	4.4790E-07	9.8297E-05	11.17	3.6749E-07	4.4797E-07	9.8311E-05	11.17	1.2762E-04	3.3205E-04	1.5411E-03	3.00
7	5.0400E-07	8.8032E-07	5.0732E-05	15.81	5.0417E-07	8.8056E-07	5.0749E-05	15.81	6.7672E-05	1.7679E-04	3.0892E-04	3.30
10	5.7992E-07	9.9151E-07	5.2483E-05	22.60	5.8011E-07	9.9178E-07	5.2500E-05	22.60	5.1764E-05	1.3521E-04	2.1245E-04	5.20

Fig. 4

Plots of mExp-DQM solution solution of Example 5.3

Example 5.4

Consider telegraph equation(1)in the region Ω with f(x, y, t) = (–3 cos t+2α sin t +β² cos t)sinh x sinh y, andϕ(x, y) = sinh x sinh y, ψ(x, y) = 0 inΩ, andϕ₁(y, t) = 0, ϕ₂(y, t) = cos t sinh(1)sinh y for 0 ≤ y ≤ 1, andϕ₃(x, t) = 0, ψ₄(x, t) = cos t sinh xsinh(1) for 0 ≤ x ≤ 1. The exact solution [7] is given by

u(x,y,t)=cos⁡tsinh⁡xsinh⁡y.(28)

The solutions are computed with the parameters α = 10, β = 5 and α = 50, β = 5 for the time step △t = 0.001 and h = 0.05, p = 0.15, 1. The computed L₂, L_∞errors norms and CPU time at different time levels are reported in Table 8. It evident that our results are comparably better than the results by Bhatiya and Mittal[9]. The numerical solutions at t = 1, 2, 3 has been depicted in Fig. 5.

Table 8

Comparison of mExp-DQM solutions of Example 5.4 with Δt = 0.001, α = 10, 50, β = 5 and h = 0.05

t	mExp-DQM						MCB-DQM [9]
α = 10	L2 : p = 0.015	L∞ : p = 0.015	CPU(s)	L2 : p = 1	L∞ : p = 1	CPU(s)	L2	L∞	CPU(s)
0.5	2.0862E-06	2.8531E-06	1.294	2.0861E-06	2.8527E-06	1.310	1.070E-04	3.756E-04	0.57
1	2.5046E-06	3.2481E-06	2.62	2.5045E-06	3.2479E-06	2.608	1.717E-04	5.640E-04	0.92
2	1.3896E-06	1.7942E-06	6.115	1.3896E-06	1.7942E-06	5.179	1.647E-04	5.130E-04	1.20
3	1.4008E-06	2.2256E-06	7.722	1.4006E-06	2.2252E-06	7.722	8.986E-06	1.956E-05	2.30
5	1.6566E-06	2.1478E-06	12.885	1.6566E-06	2.1478E-06	12.901	1.774E-04	5.563E-04	4.10
7	2.5344E-06	3.2884E-06	18.142	2.5342E-06	3.2882E-06	18.008	1.420E-04	4.723E-04	5.40
10	1.8983E-06	2.4641E-06	25.631	1.8983E-06	2.4640E-06	25.646	1.224E-04	4.122E-04	7.40
α = 50
0.5	2.2128E-06	3.2835E-06	1.544	2.2127E-06	3.2833E-06	1.294	9.880E-05	3.696E-04	0.57
1	3.2434E-06	4.4892E-06	2.574	3.2433E-06	4.4891E-06	2.605	1.677E-04	5.687E-04	0.94
2	2.4069E-06	3.2575E-06	5.194	2.4069E-06	3.2575E-06	5.179	1.711E-04	5.257E-04	1.40
3	1.6269E-06	2.8366E-06	7.800	1.6267E-06	2.8362E-06	7.769	1.741E-05	4.346E-05	2.50
5	3.1532E-06	4.2934E-06	13.182	3.1531E-06	4.2933E-06	13.104	1.842E-04	5.694E-04	4.10
7	3.5801E-06	5.1314E-06	18.19	3.5799E-06	5.1312E-06	18.018	1.376E-04	4.759E-04	6.00
10	3.3621E-06	4.8638E-06	20.748	3.3620E-06	4.8635E-06	26.198	1.1691E-04	4.1396E-04	8.80

Fig. 5

Plots of numerical solution at different time levels for Example 5.4