Efficient BTCS + CTCS Finite Difference Scheme for General Linear Second Order PDE ()
1. Introduction
We are interested in the following second order linear Partial Differential Equation (PDE):
(1)
where
designates an scalar function depending on the variables x and t. Usually x states for the space variable and t for the time. The Equation (1) can be brought in the following form:
(2)
The resolution of such an equation presents a great interest, because it governs several phenomenas in physics, chemistry, mathematics, economy, etc. [1]. It is a general form that deals with:
· Elliptic equations: It treats the Poisson equation which describes electrostatic and magnetostatic phenomenas, flow of perfect incompressible fluids, vortex and viscous incompressible flow, filtration of fluids through porous materials, etc. It concerns also the Euler-Tricomi equation that allows the study of transonic flow.
· Parabolic equations: It handles problems governed by the extended diffusion equation or the Diffusion Advection Reaction equation (of concentration, matter, temperature or electromagnetic field) or the BlackScholes equation in mathematical finance [2] [3].
· Hyperbolic equation: Equations such the extended complete wave equation and the extended complete telegraph equation are contained in (1). These two equations have important applications in electromagnetic and telecommunications [4] [5] [6] [7].
We propose to solve this equation with a simple, accurate and efficient method, by reducing the costs in time, memory space, and method’s heaviness. We aim at a good compromise between simplicity and result accuracy.
So, we will first, formulate the problem. Then, we choose a hybrid approximation scheme, using two finite differences approaches: Backward Time Centered Space (BTCS) and Centered Time Centered Space (CTCS) [8] [9] [10]. This heterogeneous scheme has the advantages of both. The first method BTCS is implicit and the second one is explicit. With their superposition, the resulting method becomes implicit. Therefore, we will use different methods for processing the matrices that result from the discretization. These methods are the algorithm of Usmani [11], which inverts any regular tridiagonal matrix directly; and the algorithm of Thomas which inverts a diagonal dominant tridiagonal matrix using the Right Hand Side (RHS) of the equation.
Subsequently, we will treat numerical experimentation to validate our proposed method. Finally, we discuss the results and give outlook for further work.
2. Problem Formulation and Meshing
We consider a scalar function
which depends on the real variables x and t. This function satisfies the partial differential Equation (3) with the given conditions.
(3)
Here, The coefficients A, B, C, D, E and F are constant or known functions depending on the variables x and t. L is a real constant. The function
is a known excitation. The functions
,
,
and
have been also given.
The following mesh has been considered: the spatial interval
is discretized in
points
, with
. The considered instants are:
. The time step is
and must be sufficiently small to allow a good, accurate and efficient resolution of the problem. An approximated value of
is to be found at point
:
. We define:
,
,
,
and
. It holds:
,
,
,
,
and
.
3. Schema BTCS + CTCS
We consider the Finite Difference method in Time Domain (FDTD); more precisely the two schemes: BTCS and CTCS [8] [9] [10]. Then, we superpose the two approaches in order to obtain a better approximation of the solution of the treated differential equation. In this way we combine the advantages of the two schemes. As one can remark, these two schemes have the same spatial discretization (Centered Space CS). Thus, we have for the derivative in x direction the following approximations:
(4)
For the BTCS scheme, we get the first and second order backward time derivatives, for
:
(5)
Combining the Equations (3)-(5), one gets:
(6)
With the CTCS scheme, we have the following first and second order derivatives, for
:
(7)
Associating the Equations (3), (4), and (7); one gets:
(8)
We define
,
,
,
,
,
,
and
; for
. Thus, we obtain the following matrix equation, considering the mesh points
and
:
(9)
This previous equation is equivalent to:
(10)
with
(11)
So, we will discuss the inversion methods of the matrix
, which permits to get the solution
.
4. General Solution
An important discussion is to be done with respect to Equation (9), particularly concerning the matrix
. Principally, this matrix must be inverted for each time
in order to get the solution at this time:
. But, if the coefficients of Equation (1) do not dependent on the time
then one and exactly one inversion is sufficient.
The algorithm of Thomas can be used for the inversion (see Appendix). With this method, we do not have an idea about the regularity of the matrix. But, it is clear that if the tridiagonal matrix
is diagonal dominant then it is regular. Of course, this property is sufficient but not necessary for the regularity.
We prefer an inversion with the Usmani’s algorithm [11], which is a general and stable method. It allows a direct inversion that does not use the right hand side (
) in the inversion process; contrarily to the algorithm of Thomas. Usmani has presented a direct and exact method to invert a tridiagonal regular matrix [11].
We can apply it for
, defining:
Then, the coefficients of the inverted matrix
are obtained with:
(12)
We get the solution:
(13)
The solution
, at point
and time
, is given by a simple matrix- vector multiplication:
(14)
5. Solution for Time-Depending Coefficients
In the case that the coefficients of Equation (1) do not dependent on the space:
,
,
,
,
and
; then the matrix (A) dependent only on the time. Its coefficients are constant at a fixed time
. The formula of the invert matrix
can be simplified [12] [13] [14]:
Defining a real number
and a complex number
in following manner:
(15)
we get
, which is the determinant of the submatrix of order i of (A); and which is of dimension (
). In the case, where
, the determinant of the matrix (A) is:
One can verify that this determinant, for
, is given by the following relations [12] [13] [14]:
(16)
The inverted matrix (B) exists when d is different of one of the following values [14]:
(17)
Outside these values of
, the matrix (A) is regular and its inverts (B) is given by the following formula:
(18)
6. Numerical Verification
For the numerical verification, we choose the equation of telegraph equation, which has been treated by several authors [4] [5] [6] [7]. We compare our results with those obtained by [4].
The following problem treated by [4] was considered and our method was applied:
(19)
By choosing dx = 0.02 and dt = 0.001, the following results were obtained, showing the
and
errors:
The results are very satisfactory because our method is not heavy and leads to a precise solution. Compared to the one used in [4], it could be said to be less precise. But it should be emphasized that the method, used in [4], is a Cubic B-splines collocation method, which is expensive in calculation. On the other hand, we used the finite difference method.
7. Conclusions
In this work, a method of solving a general linear partial differential equation has been presented. The finite difference hybrid approach (BTCS + FTCS) that has been used is simple, accurate and efficient; and is economical in terms of calculation and occupancy of the memory space.
This study can allow numerous applications of this method to several phenomena of the sciences and techniques governed by this equation. In addition, other methods could be explored to improve performance.
Appendix
For a resolution of Equation (10) with the algorithm of Thomas, the vector
,
, and
can be introduced, in order to proceed to a forward elimination. Then the solutions are obtained by backward restitution. This algorithm is presented below.