Latest revision as of 09:34, 23 June 2022

Abstract

One of the greatest challenges in the area of applied mathematics continues to be the design of numerical methods capable of approximating the solution of partial differential equations quickly and accurately. One of the most important equations, due to the hydraulic and transport applications it has, and the large number of difficulties that it usually presents when solving it numerically is the Diffusion Equation.

In the present work, a Method of Lines applied to the numerical solution of the said equation in irregular regions is presented using a scheme of Generalized Finite Differences. The second-order finite difference method uses a central node and 8 neighbor points in order to address the spatial approximation. A series of tests and numerical results are presented, which show the accuracy of the proposed method.

Keywords: Generalized finite differences, method of lines, diffusion equation, irregular regions

1. Introduction

When a substance is been transported with a speed equal to zero, i.e. no flow moves it, then, according to the transport equation [1], the initial profile of the substance should remain the same as time goes by, maintaining the substance concentration at each point in the fluid.

However, in practice, this does not happen due to molecular diffusion; this is due to the molecules being in constant motion, causing collisions and rebounds in different directions. The molecules will tend to separate, or diffuse, within the flow. In general, this movement moves from areas with a higher density of molecules to areas with a lower density; this is called diffusion and can be described by a well-known Partial Differential Equation:

where ${\textstyle \nu }$ represents the diffusion coefficient, which tells how easy a substance diffuses in a medium.

In the past years, many people have worked in different methods to achieve good numerical solutions to this equation, involving a great variety of modern and classical techniques. Nevertheless, even though a large number of numerical methods have been proposed to solve it, a large number of these have a rather limited application to real-world scenarios since they are designed for regular regions.

This is due to the fact that the spatial discretization of the diffusion equation imposes bounds on the time step size in order to achieve numerical stability. However, the use of semidiscretization schemes allows for overcoming many stability issues by using well-known and widely used one-step methods for ordinary initial value problems. For example, in Manshoor et al. [2], a Method of Lines, involving solutions with ${\textstyle a}$ Runge-Kutta method, is presented along with its stability analysis; the results presented show that it is possible to use this kind of method to compute numerical solutions of the equation, yet, the regions where the method is tested are unidimensional regions. On the side of the Generalized Finite Difference methods, Ureña et al. [3] present a scheme to achieve numerical solutions using this method; the results presented in this work show that it is possible to solve different Partial Differential Equations using finite differences over non-regular clouds of points; this work presents that the solutions obtained with this method satisfactorily match with the exact solution of the proposed tests. Some other authors, like Li et al. [4,5,6,7,8,9], and Wang et al. [10] have addressed a large number of applications. Nevertheless, the computational cost of the method as proposed in the aforementioned papers is high, so for each node of the cloud it is required to take up to twenty-six support nodes, even in 2D problems, which increases the computational cost of the implementation.

On the other hand, some variations of the presented generalized finite difference method for several transport equations, which produce satisfactory numerical solutions using a low-cost implementation for the spatial discretization were presented in [11,12,13,14,15]. Even though, the use of several straightforward time integration schemes remained an important issue to take into account.

For the case of interest of this work, it is important to obtain an approximation in generalized finite differences, for the spatial part, to the solution of the problem

where ${\textstyle \Omega }$ is a simply connected planar domain, and its boundary, and ${\textstyle \partial \Omega }$ is a positively oriented Jordan polygon, as the domain shown in Figure 1.

Figure 1. Example of an ${\displaystyle \Omega }$ domain

On the other hand, for the temporal discretization, a Method of Lines (MOL) is proposed [16]. The basic idea of the MOL is to solve a time-dependent Partial Differential Equation (PDE) by discretizing the spatial derivatives and then, integrating the semi-discretized problem as an Ordinary Differential Equations (ODE) system.

2. Proposed scheme

In order to apply a MOL for the case of the diffusion equation

it is possible to discretize the spatial derivatives applying a generalized finite differences method, for that it is convenient to considerate the approximation to the second or-er linear operator

where ${\textstyle A}$, ${\textstyle B}$, ${\textstyle C}$, ${\textstyle D}$, ${\textstyle E}$, and ${\textstyle F}$ are given functions. Within an arbitrary distribution of nodes, like the one presented in Figure 2, it is possible to approximate its value at a node ${\textstyle p_{0}=(x_{0},y_{0})}$ using values of ${\textstyle u}$ at some neighbor nodes ${\textstyle p_{i}=(x_{i},y_{i})}$, ${\textstyle i~=~1,2,\dots ,q}$ [12]. For this work, a finite difference scheme is applied at the node ${\textstyle p_{0}}$, which can be written as a linear combination as

where ${\textstyle \Gamma _{0},\Gamma _{1},\cdots ,\Gamma _{q}}$ are adequate weighs.

Figure 2. Arbitrary distribution of ${\displaystyle p_{0}}$ and its neighbors

According to Strikwerda [17] and Thomas [18], a finite difference scheme ${\textstyle L_{0}}$ is consistent with the linear operator ${\textstyle L}$ if the local truncation error ${\textstyle \tau }$ satisfies that

as ${\textstyle p_{1}}$, ${\textstyle p_{2}}$, ${\textstyle \dots }$, ${\textstyle p_{q}\to p_{0}}$.

Using the six first terms of Taylor's series expansion, up to second order, of the consistency condition (4), it is possible to obtain the system

where ${\textstyle \Delta x_{i}=x_{i}-x_{0}}$ and ${\textstyle \Delta y_{i}=y_{i}-y_{0}}$. In order o solve this linear system, it is possible to separate the first equation of the system (5)

and then, the problem defined by

can be solved using the reduced Cholesky factorization of its normal equations, as in Tinoco-Guerrero et al. [19], namely

where

The value of ${\textstyle \Gamma _{0}}$ is then obtained from Eq.(6) assuming, for the case of the diffusion equation, ${\textstyle F(p_{0})=0}$.

Now, the scheme defined by Eq.(7) can be used to approximate the linear operator

taking ${\textstyle A=C=\nu }$, and ${\textstyle B=D=E=0}$. The resulting ${\textstyle \Gamma _{i}}$ coefficients, define the Finite-Difference Scheme

for diffusion equation, where ${\textstyle u_{i}(t)}$ is the approximation to the solution in the point ${\textstyle p_{i}=(x_{i},y_{i})}$, and ${\textstyle p_{j}}$ are the corresponding neighbor nodes of ${\textstyle p_{i}}$.

An important issue to be taken into account is the number of neighbors, ${\textstyle q}$, to use in the scheme. In this paper, 8 neighbor points were taken into account following the stencil shown in Figure 3.

Figure 3. Centered stencil used in the scheme

Once a discretization of the spatial operator is obtained, following the idea presented in Tinoco-Guerrero et al. [1], the semi-discretized PDE can be rewritten as a linear ODE system in time,

where ${\textstyle m}$ is the total amount of grid points.

There exist several ways to approximate the solution of this ODE system, for example, it is possible to use a forward Euler method as

nevertheless, this kind of implementation has proven to be conditionally stable and, in some cases, the stability conditions are difficult to accomplish.

Another way to solve this system is the Method of Lines, where each ${\textstyle u_{i}(t)}$ is solved for a fixed grid node ${\textstyle p_{i}=(x_{i},y_{i})}$, i.e. it is solved by “lines”. As we are now dealing with a set of ODEs, it is possible to use a Runge-Kutta method, to solve the ODE on each line in order to obtain stable and accurate results.

In the present work it is proposed to use Runge-Kutta [20] of 2nd, 3rd, and 4th order, as follows.

Second-order Runge-Kutta method

where

Third-order Runge-Kutta method

where

Fourth-order Runge-Kutta method

where

in all the cases, ${\textstyle k}$ represents the time level.

3. Numerical tests

To show the performance of the proposed scheme, the problem of obtaining a numerical solution to the diffusion equation in four different regions was proposed. The first region, denoted as A, corresponds to the unitary square, for comparison reasons, and the other three regions, denoted as B, C, and D, are non-rectangular planar domains. All the regions were scaled to fit on ${\textstyle [0,1]\times [0,1]}$, and meshed with ${\textstyle 21\times {21}}$ nodes, following a variational procedure implemented in UNAMalla [21], then they were subdivided to obtain meshes with ${\textstyle 41\times {41}}$ nodes. The meshes for with ${\textstyle 21\times {21}}$ nodes for each region can be seen in Figure 4.

Figure 4. Meshes of the test regions with ${\displaystyle 21\times {21}}$ nodes

For all the regions, two different tests were performed:

Test 1. Following the same idea as in Tinoco-Ruiz [22], where a well-known diffusion problem is presented, the initial and boundary conditions were taken from the closed-form solution:

Test 2. For comparison purposes, one of the tests presented in Sánchez et al. [23], was also taken into account; in this case, the initial and boundary conditions were taken from the closed-form solution:

The time interval ${\textstyle [0,1]}$ was subdivided with different discretizations, the number of time steps was chosen to satisfy the classical Courant condition [24]:

where ${\textstyle \Delta x}$ and ${\textstyle \Delta y}$ are the spatial steps on each grid, taken from the meshes with ${\textstyle 41\times {41}}$ nodes, and ${\textstyle \nu =0.2}$.

The norm of the quadratic error, at a ${\textstyle k-}$the time level can be computed as

where ${\textstyle U_{i,j}^{k}}$ and ${\textstyle u_{i,j}^{k}}$ are the approximated and theoretical solutions, respectively, at the ${\textstyle i,j}$-th grid node, and ${\textstyle A_{i,j}}$ is the area of the polygon defined by the points ${\textstyle P_{i+1,j}}$, ${\textstyle P_{i,j+1}}$, ${\textstyle P_{i-1,j}}$ and ${\textstyle P_{i,j-1}}$.

Figure 5 presents a comparison of the numerical results obtained with the proposed scheme using the second-order Runge-Kutta approximation and the exact solutions for the test region A at time ${\textstyle t=0.66}$s. The approximation is presented on the left, and the exact solution on the right.

Figure 5. Comparison between numerical results (left) and exact solution (right) for test region A. Test 1 (Top) and Test 2 (Bottom)

Similarly, Figures 6, 7, and 8 present the respective comparisons for the regions B, C, and D.

Figure 6. Comparison between numerical results (left) and exact solution (right) for test region B. Test 1 (Top) and Test 2 (Bottom)

Figure 7. Comparison between numerical results (left) and exact solution (right) for test region C. Test 1 (Top) and Test 2 (Bottom)

Figure 8. Comparison between numerical results (left) and exact solution (right) for test region D. Test 1 (Top) and Test 2 (Bottom)

Figures 9 to 12 present the maximum value of the error ${\textstyle \|e^{k}\|_{2}}$ computed for all the regions. For each figure, both tests were performed in the regions meshed with ${\textstyle 21\times {21}}$ nodes and ${\textstyle 41\times {41}}$ nodes.

Figure 9. Maximum ${\displaystyle \|e^{k}\|_{2}}$ computed for test region A with ${\displaystyle 21\times {21}}$ nodes (Left) and ${\displaystyle 41\times {41}}$ nodes (Right)

Figure 10. Maximum ${\displaystyle \|e^{k}\|_{2}}$ computed for test region B with ${\displaystyle 21\times {21}}$ nodes (Left) and ${\displaystyle 41\times {41}}$ nodes (Right)

Figure 11. Maximum ${\displaystyle \|e^{k}\|_{2}}$ computed for test region C with ${\displaystyle 21\times {21}}$ nodes (Left) and ${\displaystyle 41\times {41}}$ nodes (Right)

Figure 12. Maximum ${\displaystyle \|e^{k}\|_{2}}$ computed for test region D with ${\displaystyle 21\times {21}}$ nodes (Left) and ${\displaystyle 41\times {41}}$ nodes (Right)

4. Conclusions

The numerical results show that the proposed method of Lines applied to the diffusion equation produces satisfactory numerical solutions. In the tests carried out, no spurious oscillations or instabilities were perceived. In addition, the results show that it is not necessary to “work more” by using higher Runge-Kutta methods; this scheme accomplishes better results when second-order Runge-Kutta is used to solve the Ordinary Differential Equations system.

It is possible to appreciate that with the presented discretizations, even with the fourth-order Runge-Kutta method, acceptable numerical results can be obtained. Additionally, the numerical results shown in Figures 9 to 12 show that a refinement of the spatial mesh can improve the approximations carried out with the method; as expected in these kinds of methods. It is worth mentioning that, for all the tests, a minimum number of time steps was used, taking into account the classical theory for the diffusion equation, which makes the computational cost of this method low, since it does not require making very small temporary discretizations, as in other cases.

Furthermore, the proposed scheme can produce stable results for non-standard initial conditions in highly irregular regions. For example, in the following videos

• https://youtu.be/soUxe5uUY_U
• https://youtu.be/479YeBjqmSs

solutions of the diffusion equation are presented on domains that are geometrical approximations of real geographical locations, where the boundary conditions are fixed as ${\textstyle 0}$ and the initial condition is stated as:

It is possible to see that even for these conditions in these regions, the scheme produces stable results that show the expected behavior.

An important remark is that, even when in this work structured convex grids were used, the development of the method doesn't take into account a particular data structure, i.e. this method can be used not only on structured meshes but also as a meshless method, that would be developed as future work.

Acknowledgments

We want to thank AULA CIMNE-Morelia and CIC-UMSNH for the financial support for this work. Thanks to Universidad Vasco de Quiroga for lending the facilities to perform the final part of the work done for the paper.

Special thanks to the ``Laboratorio de Matemáticas Aplicadas of the ``Facultad de Ciencias Físico-Matemáticas at ``Universidad Michoacana de San Nicolás de Hidalgo for lending their computational infrastructure.

The authors would also like to thank the reviewers for their valuable feedback and suggestions to improve the present work.

References