A New Efficient Method for the Numerical Solution of Linear Time-Dependent Partial Differential Equations

This paper presents a new efficient method for the numerical solution of a linear timedependent partial differential equation. The proposed technique includes the collocation method with Legendre wavelets for spatial discretization and the three-step Taylor method for time discretization. This procedure is third-order accurate in time. A comparative study between the proposed method and the one-step wavelet collocation method is provided. In order to verify the stability of these methods, asymptotic stability analysis is employed. Numerical illustrations are investigated to show the reliability and efficiency of the proposed method. An important property of the presented method is that unlike the one-step wavelet collocation method, it is not necessary to choose a small time step to achieve stability.

In solving time-dependent problems, Legendre wavelets are often used for spatial discretization.Different techniques are implemented for time discretization.In some articles, Legendre wavelets are also applied for time discretization.Therefore, the collocation points should be defined for both time and spatial variables.Also in this technique, multi-dimensional wavelets should be used to approximate required functions, which deal with large matrices and require large storage space.For example, readers can refer to [9].
There are many contexts that use collocation methods in solving functional equations.For example, Luo et al. [12] presented three collocation methods based on a family of barycentric rational interpolation functions for solving a class of nonlinear parabolic partial differential equations.Furthermore, for solving a class of fractional subdiffusion equation, Luo et al. in 2016 [13] used the quadratic spline collocation method.
Another path for time discretization uses a finite difference method.Islam et al. [10] used a fully implicit scheme, which is based on the first-order Taylor expansion.Yin et al. [11] employed the θ-weighted scheme for nonlinear Klein-Sine-Gordon equations.Stability is the important point in using finite difference methods.Thus, methods that are first-order accurate in time might be inappropriate.
Here, we exploit the three-step finite element method for time discretization [14][15][16].For the suitable differentiable function F(t), these three steps are defined as follows: It can be shown that the above equations are equivalent to the third-order Taylor expansion.Therefore, this method is third-order accurate in t.The first idea of using these three steps has been demonstrated by Jiang and Kawahara [14].Equations ( 1)-( 3) are usually accompanied by the Galerkin finite element method, which is known as the three-step Taylor-Galerkin method [17].Kumar and Mehra [2] proposed a three-step wavelet Galerkin method based on the Daubechies wavelets for solving partial differential equations subject to periodic boundary conditions.In this paper, motivated and inspired by the ongoing research, we develop a new effective method, which combines the Legendre wavelets collocation method for spatial discretization and the mentioned three steps for time discretization in the numerical solution of a linear time-dependent partial differential equation subject to the Dirichlet boundary conditions.We call this method the three-step wavelet collocation method.Furthermore, we explain the asymptotic stability of the proposed method.
The organization of this paper is as follows.In Section 2, fundamental properties of the Legendre wavelets are described.The three-step wavelet collocation method is presented in Section 3. The analysis of asymptotic stability is performed in Section 4. Some numerical examples are presented in Section 5. Finally, Section 6 provides the conclusions of the study.

Basic Properties of Legendre Wavelets
Legendre wavelets are defined on the interval [0, 1] as follows [5]: where k can assume any positive integer, where: and c l,m =< f (x), ψ l,m >, in which < ., .> denotes the inner product.The derivative of the vector Ψ(x) can be expressed by: where D is the 2 k (M + 1) operational matrix.Mohammadi and Hosseini obtained D and the operational matrix for the n-th derivative: in [5].

Three-Step Wavelet Collocation Method
In this section, we explain the main structure of the three-step wavelet collocation method.

Time Discretization
Consider the following linear time-dependent partial differential equation: with the initial condition: and boundary conditions: Assume that n ≥ 0 and ∆t denote the time step such that t n = n∆t, n = 0, 1, • • • N t .By using the Taylor expansion, the value of the function u(x, t) at the time t n+1 can be expressed as follows: where the symbols u n and ( ∂u ∂t ) n represent u(x, t n ) and ∂u ∂t (x, t n ), respectively.
We can use the first-order Taylor expansion for time discretization and Legendre wavelets for spatial discretization [10].We call this method the one-step wavelet collocation method.In addition, the time derivative in the given differential equations is approximated by Euler's formula: and therefore, we have semi-discrete equation: The three-step Taylor method for time discretization is derived by applying a factorization process to the right side of Equation (10) as follows: where the symbol I is the identity operator.Now, using Equation (11) and employing a new notation, the three-step Taylor method is obtained as follows: ) n (12) It should be noted that u n+ 1 3 , u n+ 1 2 and u n+1 represent the computed solution at time level ) and (t n + ∆t), respectively.

Spatial Discretization
After time discretization, the spatial derivatives of u(x, t) are approximated by Legendre wavelets.The collocation method is utilized in this part.Let the unknown solution u(x, t n ) be expanded by: According to Equation (15), we use only one-dimensional Legendre wavelets to approximate the solution.The solution dependence on the time variable is specified by the coefficient c n l,m .In other words, the vector coefficient C n is calculated at time t n .Therefore, the approximation solution at time t n+ 1 3 can be written as follows: We can also approximate f (x, t) at time t n as: where the vector (F n ) T is given by Equation (4) at time t n .Substituting Equation ( 6) into Equation ( 12) results in: Now, by using the operation matrix of the derivative and Equation ( 5), we have: Then, substituting Equations ( 15)-( 17) and ( 19) into Equation ( 18) yields: By using boundary conditions and Equation ( 16), the following equalities are satisfied: Considering the initial Condition (7), we have: By substituting Equation ( 20) in (2 and using Equations ( 21) and ( 22), we can obtain a linear system of equations with 2 k (M + 1) unknown variables, c n+ 1 3 l,m , which can be written in matrix form: where A and B are 2 k (M + 1) × 2 k (M + 1) and 2 k (M + 1) × 1 matrices, respectively.Since the vector C 0 is obtained from Equation ( 23), all entries of B are known.
The above system of linear equations can be solved by numerical methods.Here, for square matrix A, we use LU decomposition to solve the linear System (24) with partial pivoting.In this method, the square matrix A can be decomposed into two square matrices L and U such that A = LU, where U is an upper triangular matrix formed as a result of applying the Gauss elimination method on A, and L is a lower triangular matrix with diagonal elements being equal to one.Solving the system AC n+ 1 3 = B is then equivalent to solving the two simpler systems Ly = B and UC n+ 1 3 = y.The first system can be solved by forward substitution, and the second system can be solved by backward substitution.Solving the linear system with triangular matrices makes it easy to do calculations in the process of finding the solution.Since the Gaussian elimination can produce bad results for small pivot elements, we adopt the partial pivoting strategy.In this strategy, when we are choosing the pivot element on the diagonal at position a ii , locate the element in column i at or below the diagonal that has the largest absolute value, and make it as the pivot at that step by interchanging two rows.Applying this strategy to our matrix avoids any distortion due to the pivots being small.For more details, readers can refer to [18,19].
After solving this system and determining C n+ 1 3 , we exploit Equation ( 13) to find C n+ 1 2 .In addition, there is a similar process that results in: and: where are the same collocation points used in the previous step.A matrix form of Equations ( 25)-( 27) can be displayed as: where the dimension of the square matrix P and column vector Q is 2 k (M + 1).Since the vector C n+ 1 3 is obtained from the previous step, all entries of Q are known.Similarly, we use LU decomposition to solve the above system.Finally, by implementing similar analysis in the two previous steps, using Equation ( 14) with boundary conditions and exploiting C n+ 1 2 , the vector C n+1 can be specified in each step for n = 0, 1, 2, • • • .Therefore, we can obtain the numerical solution, u(x, t n ), in any time t = t n .

Stability Analysis
For stability analysis, we use the asymptotic (or absolute) stability of a numerical method, which is defined in [20].In a numerical scheme, when we fix the final time t = n∆t and let n → ∞, we want the corresponding numerical solution to remain bounded; a scheme satisfying this property is called stable.Therefore, a stability analysis needs a restriction on the mesh size ∆t.In practice, we can only choose a finite and proper mesh size.It is then important to study the region of absolute stability in order to to choose the proper mesh size in practical computation.
Let us start from the typical evolution equation: where the non-linear operator f contains the spatial part of the partial differential equation.Let us abbreviate u(x j , t n ) by u n j .We shall approximate u n j by U n j .Following the general formulation of the proposed method, the semi-discrete version is: where u j is the spectral approximation to u, f j denotes the spectral approximation to the operator f and Q j is the projection operator, which characterizes the scheme.Let us set U(t) = Q j u j (t).Then, the previous discrete problem can be written in the form: As is often done, we confine our discussion of time-discretizations to the linearized version of (28): where L is the diagonalizable matrix resulting from the implementation of spectral method on the spatial variable of the partial differential equation.According to different contexts, the time discretization is said to be stable if U n , the computed solution at the time t n = n∆t, has been upper bounded, i.e., there exists a constant M such that: In many problems, the solution is bounded in some norm for all t > 0. In these cases, a method that produces the exponential growth allowed by Estimate (30) is not practical for long-time integrations.For such problems, the notion of asymptotic (or absolute) stability is useful.Definition 1.The region of absolute stability of a numerical method is defined for the scalar model problem: to be the set of all λ∆t such that U n is bounded as t → ∞ [20].
Finally, we say that a numerical method is asymptotically stable for a particular problem if, for sufficiently small ∆t, the product of ∆t times every eigenvalue of L lies within the region of absolute stability.In the following items, we summarize some remarkable characteristics of absolute stability [21]: 1.An absolutely stable method is one that generates a solution u n that tends to zero as t n tends to infinity, 2. A method is said to be A-stable, if it is absolutely stable for any possible choice of the time-step, ∆t, otherwise a method is called conditionally stable.3. Absolutely stable methods keep the perturbation controlled, 4. The analysis of absolute stability for the linear model problem can be exploited to find stability conditions on the time step when considering some nonlinear problems.
Since the three-step Equations ( 12)-( 14) are equivalent to the third-order Taylor expansion, to demonstrate the stability region and achieve the stability condition, we use Equation (10).For simplicity, consider Equation (6), where µ = 0 and f (x, t) = 0.Then, successive differentiations of the obtained equation indicate that: In Equation (32), we use Euler's formula to avoid the third-order space derivatives, as it is used in the finite element context [22].By rearranging Equation ( 10) and substitution of Equations ( 31) and (32), we have the semi-discrete equation: After applying the wavelet collection method, Equation (33) transforms into the following equation: where: and {x i } are the collocation and boundary points.Here, the matrix L, which is introduced in Equation (29), is defined as There is a similar process to the one-step method.Lambert provided an explanation for how to draw the stability region.Readers can refer to [23], Chapter 3. Briefly, we can plot the region of absolute stability, R L , by the meaning of the first and second characteristic polynomials.If we set, ĥ = λ∆t, the region of absolute stability is a function of the method and the complex parameter ĥ only, so that we are able to plot the region R L in the complex ĥ-plane.
First of all, we can write Equation (34) as a usual linear multi-step method given by: where k is the number of steps required for the method, and α j and β j are constants subject to the conditions: According to Equations ( 34) and (35), we have: Afterward, the first and second characteristic polynomials are defined as follows, respectively: where ξ ∈ C is a dummy variable.Using the values of k and {α j , β j } 1 j=0 in (36), for the proposed method, we have: Then, we plot the boundary of R L , which consists of the contour ∂R L .The contour ∂R L in the complex ĥ-plane is defined by the requirement that for all ĥ ∈ ∂R L , one of the roots of π(r, ĥ) := ρ(r) − ĥσ(r) has modulus one, that is, it is of the form r = exp(iθ).Thus, for all ĥ in ∂R L , the identity: π(exp(iθ), ĥ) = ρ(exp(iθ)) − ĥσ(exp(iθ)) = 0, must hold.This equation is readily solved for ĥ, and we have that the locus of ∂R L is given by: ĥ = ĥ Finally, we use (37) to plot ĥ(θ) for a range of θ ∈ [0, 2π] and link consecutive plotted points by straight lines to get a representation of ∂R L .
Therefore, according to Lambert's book and the above explanations, the stability region of the three-step and one-step wavelet collocation methods is the circle with center (−1, 0) and radius one.Therefore, these methods will be stable if the eigenvalues of the corresponding system and ∆t satisfy Re(λ j ∆t) ∈ [−2, 0].

Numerical Examples
In this section, some numerical examples in the form of Equation ( 6) with initial and boundary Conditions ( 7)-( 9) are discussed.The error function is defined as the maximum error L ∞ : where u is the approximate solution, which is obtained by the proposed method and {x i } are the Gauss-Legendre and boundary points.All programs have been performed in MATLAB 2016.
In general, the numerical results are sensitive to the selection of parameters such as time step ∆t, final time t and parameters of wavelet order M and k.In the following examples, we choose t = 0.5.Although the one-step wavelet collocation method needs less calculation, the three-step wavelet collocation method is more successful in finding the numerical solution.Furthermore, we compare our method with the three-step method proposed in [17].They used the finite element method with standard linear interpolation functions for spatial discretization and the same three-step formula for time discretization.Numerical results show that utilizing Legendre wavelets with these three steps gives higher accuracy.
Numerical results for M = 6 and M = 8 are reported in Tables 1 and 2, respectively.The first two rows of Table 1 show the obtained results for k = 2.As can be seen from these rows, a change in the length of the time step makes the one-step method fail to find an approximate solution.In other words, the one-step method is unstable for ∆t = 0.01, while the three-step method gives high accuracy.There is a similar analysis for other parameters.
The exact and three-step approximate solutions of u(x, t) are shown in Figure 1, where M = 8 and k = 2. Figure 2 shows the absolute stability region based on the three-step and one-step methods with the position of λ j ∆t, where {λ j } 2 k (M+1) j=1 are the eigenvalues of corresponding matrix L. This figure is drawn for M = 6 and ∆t = 0.01.As can be seen in this figure, there are some eigenvalues for the one-step method with Re(λ j ∆t) / ∈ [−2, 0]; however, the stability region of the three-step method includes all {λ j ∆t} . Therefore, the one-step method is not stable for ∆t = 0.01, while the three-step method is stable.Example 2. In this example, we consider Equation (6) with ν = 1/π 2 , µ = −4, f (x, t) = 0, g(x) = sin(πx), h 0 (t) = 0 and h 1 (t) = 0.The exact solution of this problem is u exact (x, t) = e −5t sin(πx) [3].Table 3 gives the comparison between the three-step wavelet collocation method and one-step wavelet collocation method for M = 8.We can see from this table that the one-step wavelet collocation method tends to be unstable with a small change in time length.However, the three-step method keeps its stability for bigger ∆t.
The exact and approximate solutions for the three-step wavelet collocation method are shown in Figure 3.The stability region for both three-step and one-step methods by choosing M = 8, k = 1 and ∆t = 0.006 is shown in Figure 4.As can be seen from this figure, there are some eigenvalues in the system of the one-step collocation method with Re(λ j ∆t) / ∈ [−2, 0].Therefore, this method is not stable for ∆t = 0.006.In general, for k = 1, the one-step collocation method shows a stable and accurate result if ∆t 0.001, while the three-step collocation method is stable for ∆t 0.006.Example 3.For the last example, consider Equation (6) with ν = 1/π 2 , µ = 0, f (x, t) = sin(πx), g(x) = sin(πx) + cos(πx), h 0 (t) = e −t and h 1 (t) = −e −t .The exact solution of this problem is u exact (x, t) = sin(πx) + e −t cos(πx) [3].
Table 4 shows the maximum error for some different values using the one-step and three-step wavelet collocation methods.It is clear from this table that the one-step wavelet collocation method is unstable, while the three-step wavelet collocation method has a more precise response.
Figure 5 shows the three-step approximate solution and the exact solution.The stability region for both three-step and one-step methods by choosing M = 8, k = 1 and ∆t = 0.006 is shown in Figure 6.

Conclusions
In this paper, we proposed a new numerical method for a linear time-dependent partial differential equation.We called this method the three-step wavelet collocation method.In this method, time discretization was performed prior to the spatial discretization.These steps are equivalent to the third-order Taylor expansion; therefore, this method is third-order accurate in time.For spatial discretization, Legendre wavelets were used, which resulted in good spatial accuracy and spectral resolution.A comparison between the proposed method and other methods was presented.The theoretical aspect of absolute stability was discussed.This stability is based on λ j ∆t, where {λ j } are the eigenvalues of the corresponding system.Numerical performance shows that the three-step method leads to an effective time-accurate scheme with an improved stability property.
The proposed method can be easily implemented for other cases of time-dependent partial differential equations.For example, extending our results with Legendre wavelets in solving nonlinear partial differential equations or fractional equations is worthwhile for future contribution.

Figure 1 .Figure 2 .
Figure 1.The exact and approximate solutions of Example 1 in the case M = 8, k = 2 and ∆t = 0.005.

Figure 3 .Figure 4 .
Figure 3.The exact and approximate solutions of Example 2 in the case M = 8, k = 2 and ∆t = 0.0016.

Table 1 .
The L ∞ error of Example 1 in M = 6.

Table 2 .
The L ∞ error of Example 1 in M = 8.

Table 3 .
The L ∞ error of Example 2 in M = 8.