Abstract
In this paper, we give a new strategy to extend a numerical approximation method for two-dimensional reaction-diffusion problems. We present numerical results for this type of equations with a known analytical solution to qualify errors for the new method. We compare the results obtained using this approach to the standard finite element approach. The proposed method is adequate even with the singular right-hand side of type Dirac.
1. Introduction
This work is motivated by the extension of a numerical approximation method constructed in [1, 2] and called ziti. Note that when the mathematical model is used to approximate a real concrete problem, most of the numerical schemes have undesirable oscillations, especially near the domain’s boundary or near the physical phenomena, such as shock waves, relaxation waves, composite waves, blowup [3–6], and boundary layer which can exist in mathematical models, depending on the parameters involved, the nonlinearity [7, 8], the coupling, the boundary, and the nature of the domain. Note that the link between different parameters can affect the degree of regularity of the solution (local, global, explodes in time, singular, etc.) [9–11]. The ziti method has been constructed to exhibit such solutions regardless of their degree of regularity. This method has been tested on several models, which present or do not present singularities, and has been compared with exact solutions and classical methods (finite elements, finite volume, finite differences, particular method, and spectral method), but only in the case of a cubic domain (interval, square, cube, or even Cartesian domain in in the general case), especially the differential problems for which the application of other classical methods does not guarantee the transition from a regular behavior to another singular (i.e., giving a solution of the problem, by detecting the singularity when it exists. In this case, the method ziti has shown its efficiency, robustness, and its convergence. It could detect the blowup, the boundary layer, and also it approaches the singular solution avoiding unwanted oscillations [1, 2]. In this paper, we construct another strategy cited in [2], but for a domain more complicated than the Cartesian one, for example, a disk in the two-dimensional case, or any ball in the general case. The major goal of our investigation is to apply this numerical method in two types of , even if the right-hand side could present a singularity of type Dirac. In fact, despite the relatively simple model represented by the transient diffusion equation (especially with bad sign), its numerical solution can still be a serious challenge when the solution diffuses with steep boundary layers, or when we start with very sensitive initial data. The numerical study of the diffusion problems has been the aim of several papers. The third boundary-value problem for a loaded heat equation in a -dimensional parallelepiped is considered in [12]. Considering that the problem is time-dependent, evaluating the solution process over refined meshes and for many hundreds or thousands of time steps can become a serious numerical burden. This can be a significant numerical challenge in two-dimensional problems and even more so in three dimensions [13–16]. Our paper mainly completes the investigations of [12–18], in the case of time-dependent diffusion problems. The difference is the use of a new numerical method, which gives us an admissible approximated solution. Considering that the diffusion problem is time-dependent, evaluating the solution process over refined meshes and for many hundreds or thousands of time steps can become a serious numerical burden. This can be a significant numerical challenge in two-dimensional problems and even more. In this work, we will present a numerical solution. Our objective is to solve such kind of problems with a very simple algorithm, which is a big profit, especially for the CPU time. The main goal of the method is to approach a function with several variables, to integrate it in a given domain, and to resolve numerically Partial and Ordinate Differential Equations ( and ). This method is based on the classical variation formulation of Galerkin and the most important step is the construction of our orthonormal family from the test function with compact support, defined bywhere and is the space dimension. This function is especially used in numerical analysis, distributions, and functional analysis. In [1, 2, 19], all the mathematical tools of ziti construction in the multidimensional Cartesian case were detailed. The ziti method is composed of two steps: it starts like the Galerkin (or spectral) method with a variational formulation and hence there is the need for an orthonormal basis (total or not) of an adequate space. Then, we use the roots of the base instead of the nodes of the mesh, which will give ziti a resemblance with the particle method. At first glance, the ziti method looks as if it is the finite element method, but after using the roots, the method will follow the technique of the particular one.
The main aim of this paper is the construction of the ziti method when the domain is a disk in the two-dimensional case (in general, a multidimensional ball). To generalize this method, we opt for two strategies: the first one consists in sweeping all the disk with segments, in the two directions, as shown in Figure 1 and reconstructing our basis functions in every segment. To test this strategy, we apply the resulting tools to calculate numerically the integrals and to solve partial differential equations (two tests will be detailed: an elliptic equation “the Poisson problem” and a parabolic one “the heat equation with bad sign”). The second strategy is a direct use of the polar parametrization of a disk; we will show that this strategy is also efficient and gives a good approximation even when a singularity is forced from the right-hand side term.
The outline of this work is as follows: in Section 2, we present our approximation method in the monodimensional case and for several dimensions. In Section 3, we present the mathematical tools of construction, which permits the application of this method in the Cartesian case and the calculation, numerically, of some integrals defined in a disk domain. Section 4 is devoted to the construction of the method’s fundamental elements, using polar coordinates. Like the previous section, one of the most important parts is the numerical integration using our method and in the two cases, we will compare the exact value of an integral by the numerical one, obtained by the method. In Section 5, we apply our approach to find the numerical solution of the Poisson problem and the heat equation. Our goal is to compare the solution obtained by the method with a given analytical one defined in a disk domain and to calculate the error in . In the next, we present an approximated solution using the finite element method and we compare it with ours. To see that the method used in the current work is efficient, we treated two types of with an exact solution which converges to the Dirac mass (exact solution for the Poisson model and initial data for the heat transfer with bad sign). The function used is an approximation of the Dirac mass, and we will see that the method can detect this kind of singularities (always in a disk domain). Finally, in Section 6, we finish with some concluding remarks.
2. Overview of the Monodimensional Construction
In this section, we will present the necessary mathematical tools for the construction of our method, in the monodimensional case.
2.1. Construction of the Linearly Independent Family
From the function defined by (1), we define the set :where C\coloneq is a constant of normalization. This sequence converges to Dirac measure in the sense of distributions. Fundamental results of construction in the monodimensional case are given as follows:
First, we take a uniform mesh of [a, b] with the step where is a given integer such that .
We construct the family as follows:
The curves of are shown in Figure 1.
2.2. Construction of the Orthonormal Family
Let us consider the Hilbert space with the usual scalar product . Observe that the family is linearly independent, then using the Gram–Schmidt process, we construct a unique orthogonal family, noted as satisfying the following relation:with . Using the orthogonality of the family , we can find the second definition of given byand the recurrence application of definition (4) gives the following formula:
Let the normalization of . The curves of are shown in Figure 2.
We can verify the following assumption.
For , admits reel root(s) noted as in the interval , more precisely, .
The method permits approaching a given function and its integral by the following relations:
The main idea of the ziti method is to use the roots of the functions which satisfy
If we take in (7a), we obtain
Relations (9a)–(9d) will be used frequently to construct the numerical scheme associated with the nonlinear system studied. To reduce the iterations number, the authors of [1] proved the following optimality result:where [.] denotes the integer part. In particular, for , we conclude that the parameters are nearly stationary from a certain rank, which reduces considerably the number of iterations. Using as a root of , we can define the parameter by (for more details of the construction of the ziti method, see [1, 2].)
3. The First Strategy: Cartesian Coordinates
In this section, we are interested in the extension of the method, when is a disk centered in the origin (or the ball in the multidimensional case). Our strategy is inspired by the monodimensional case. We sweep the inside of the domain by a set of intervals horizontally and vertically (see Figures 3 and 4 ). Note that the distance between every two successive chords is constant. We choose a fixed number of nodes to subdivide every interval. As the length of each chord varies according to the position of the segment inside the disk, we will adapt our subdivision such that each interval has the same number of nodes already fixed. Therefore, for each node, we define the associated basis function (denoted ); after that, we compute the roots (denoted by ) of every element of this basis (; see Figure 5).
Remark 1. Note that, for every fixed vertical level (resp., the horizontal level ), every internal segment is limited by and (resp., and ; see Figure 3). For simplicity, the step of the horizontal subdivision will be noted (resp., the vertical step will be noted ).
We present an algorithm to calculate the internal nodes (Algorithm 1).
| ||||||||||||||||||||||||||||||||
3.1. Construction of the Orthonormal Set
For every node (resp., ), we associate the function (noted as if there is no ambiguity) (resp., the family will be noted as ) defined bywhere is the step of construction in the horizontal interval of indication , which describes the distance between the nodes and is the step of subdivision in the vertical interval of indication , which describes the distance between and
It is simple to see that the family is linearly independent, so we can construct an orthogonal family by using the Gram–Schmidt process, in the space (construction in every internal interval of the domain , horizontally and vertically), verifying the following relation:
which will be reduced in the following theorem, already proved in the monodimensional case (see [1, 2]).
Theorem 1. The orthogonal family (vertically and horizontally) verifies the following recurrence relation:where is the usual scalar product in the Hilbert space .
Corollary 1. The family and the set defined in Theorem 1 verify the following relations:
3.2. Fundamental Results. Numerical Integrations
In this paragraph, we are interested in the approximation of integrals, where the domain is the unit disk , using the horizontal and vertical test functions, as well as the roots, verifying the following relations:where are the basis functions in the horizontal dimension (resp., are the basis functions in the vertical dimension) are the roots of (resp., are the roots of )
In this section, we are interested in the approximation of a double integral, defined in a disk domain; using (7a), we obtain the following results:
Theorem 2. Let , let be a given function in , and let N denote the roots number; therefore, we have the following approximations:where we take and .
Proof. In [1, 2], the authors approximated integral formulas in the one-dimensional case using the ziti method as follows:We recall the principal points of the proof which is inspired by the spectral method:with . Using the roots of , we obtainbecause of the following orthogonality result:Multiplying (18) by and integrating it over , we obtainUsing the orthonormal property of , (21) becomesand, therefore,In the bidimensional case, using the Cartesian strategy,so if we take , where is the nodes number in every segment, horizontally and vertically, we will have the following approximation formulas:Here, we take, and .
In Table 1, we present some numerical tests of integration. We compare the exact value noted , with the numerical approximation obtained by the method in the Cartesian case. Each interval is subdivided into 50 nodes.
4. Second Strategy: Polar Coordinates
In this section, we built all the necessary elements for the approximation , using polar coordinates. The domain is represented using the polar coordinates, with the following parametrization:
The polar set is defined by
4.1. Fundamental Results: Numerical Integration
To test the previous strategy, we present in Tables 2 and 3 some numerical tests. We compare the exact value with the numerical one, using the method in the polar case.
In Table 3, we can treat even the generalized integrals, for example. , which is in fact an operation of Riemann integral; we found a good approximation using the roots. For the two last examples, and , other approximation methods (e.g., Simpson, Trapeze\“enleadertwodots”) did not give any result, which is an important point for our construction.
5. Numerical Applications
5.1. Elliptic PDE Case: Nonevolutive Diffusion Problem in a Unit Disk
5.1.1. The Cartesian Case
In this section, let us consider a partial differential equation, which admits an exact solution and we will compare it with the numerical one, using our method in the Cartesian case. Let . The problem studied is given byand, for simplicity of numerical simulations, we take , with a given analytical solution , where the source term takes the value .
5.1.2. The Strong Discretization
The first step to approach the previous problem is to multiply equation (7a) by a test function and to integrate the result over the domain , which gives
Using Theorem 2, we obtain the following equality:
The next step consists of approaching the second derivative, which gives us the following scheme:where denotes the approximation at the root and is the nodes number in every internal segment (horizontally and vertically). In the end, we will have a global matrix, with lines and columns, defined as follows:where is a tridiagonal matrix, defined byand is a diagonal matrix defined bywhereand, therefore, we should resolve a simple system in the form , when is the global matrix defined previously, is the unknown vector of size , and is the source term vector of size .
Remark 2. To complete the resolution of the previous system, we must add boundary conditions (homogeneous Dirichlet in this case).
5.1.3. Numerical Results
Let . In this case, we fix the points number in every single segment and we vary the subdivision step. It is clear that the minimum of all the steps is obtained at the first segment (horizontally or vertically) and the maximum is on the segment confused with the diameter of the disk (i.e., for two different intervals, horizontally or vertically, the associated step is not the same). We are interested in the shape of the approximated solution with the scheme, using the Cartesian coordinates and the segments approach. For a fixed node’s number in every segment (horizontal or vertical), , the numerical implementation of the scheme gives us an approximated solution, which is near the exact one, given by , when . The two solutions are illustrated in Figures 6 and 7.
In Table 4, we present the error between the exact and approximated solution, with different values of nodes number .In Table 4,
We remark that the committed error between the exact solution and the computed one using the ziti method decreases when the number of points on each interval increases, which shows the convergence of the solution.
We present, in the following subsection, a comparison between approximated solutions using the finite element and ziti methods.
5.1.4. Comparison with the Finite Elements Method
The finite element method (FEM) is a widely used analogy to resolve some types of partial differential equations. A large class of works was already done to resolve the Poisson problem using FEM (see [20, 21]). The starting point for the FEM is a PDE expressed in a variational form. The basic recipe for turning a PDE into a variational problem is to multiply the equation by a test function and to integrate the resulting expression over all the domain : it is the common step between the Galerkin analogy and ziti. In this part, we present the approximated solution of Poisson’s problem defined in (7a), using the finite element method, which is illustrated in Figure 8.
5.1.5. The Polar Case
Now, we consider the same partial differential equation defined before, which admits a polar analytical solution; we will compare it with the approximated one founded using the method in the polar case.
Let ; the strategy presented consists in using the results of approximation in the monodimensional case and taking into consideration the following function basis:
The problem presented in the previous Subsection 5.1.1 is equivalent to the polar one, expressed as follows:withwhere is the source term defined in Cartesian problem, with an exact polar solution . Note that, the roots of the basic functions will be noted; and are those associated with .
To obtain a numerical scheme using the method, we should multiply equation (38a) by a test function and after that, we use the strong result of approximation (16), which giveswith
Therefore, the goal is to find solution to the polar problem given in (38a) and (38b). Like the Cartesian analogy, we resolve in this case a simple system in the form ; when is the global matrix defined previously, we should just add the polar terms and in the corresponding terms of the matrices and , in the polar unknown vector of size , and is the source term vector of size .
5.1.6. Numerical Tests
Let be the root’s number for the radius and the other one for the angles. Two solutions, exact and approximated, are illustrated in Figures 9 and 10 .
5.1.7. Numerical Error
Table 5 shows us the error between exact and approximated solution, using different and . In Figure 11, we present the evolution of the error between exact and approximated solution, with several values of the nodes number for the diffusion problem.
(a)
(b)
(c)
(d)
Table 5 shows how the error is reduced when the number of nodes increases from to (resp. ) which is illustrated in Figure 11. The results obtained in this study suggest that having just 60 radial nodes leads to an efficient solution.
5.2. Detection of Blowup in the Stationary Diffusion Problem
To evaluate the efficiency of the ziti approach, we consider the following diffusion problem:where the source term is given bywhere is the normalization constant defined by (1) and is the regular function given by
It is easy to see that the analytical solution of (42) and (43) is given by
This function converges to Dirac mass in the weak sense when goes to 0. We present in Figure 12 the exact solution (resp., Figure 13 presents the approximated solution obtained by the scheme presented in (31)), to see that the ziti method detects this kind of singularity.
5.3. Parabolic PDE Case: Heat Equation
This section is devoted to the application of the method on a diffusion equation, in a domain , where . The heat equation describes the distribution of heat (or variation in temperature) in a given region over time. For a function, (resp., ) of two spatial variables in the Cartesian case ( in the polar case) and the time variable t, the heat equation is given by
where is a given time interval and is the unit disk with boundary . Here, denotes the space variables, is the time variable, is the diffusion coefficient, and represents the effect of internal source terms. Using the same analogy applied in the previous sections, we multiply equation (46a) by a test function , after we integrate over the domain . It remains just the direct application of our approximations formulas given in (3.3). First, the time domain is divided into uniform subintervals with duration for . To denote the value of in the root at time , we use the notation for . At , we find directly the following explicit scheme:where
Our goal is to find an approximated solution, near the exact one, which verifies the boundary conditions. The functionis an exact solution of the heat equation, with
The above example provides a starting point for numerically solving diffusion problems on the disk. These ideas can also be used on the sphere. Figures 14 and 15 show the allure of the exact and approximated solution, at a given finite time. Note that, the stability condition CFL is numerically well verified.
In Table 6, we present the error between exact and approximated global solution at a given finite time. Note that the CPU time necessary to obtain an approximated solution at (when the number of nodes in every interval, horizontally and vertically, is 200) was around 150 seconds, which means that the method is fast and robust and gives us the numerical results in a very short period of time. Table 6 shows how the error is reduced when we refine the mesh from a number of nodes to .
For checking the stability of our numerical scheme, we use the following strategy:in which is the matrix of numerical solution at iteration. Errors obtained corresponding to using the ziti method with for the heat transfer problem are reported in Table 7.
5.3.1. Comparison with Finite Element Method
Since the low order FEM is often used in industrial applications to solve heat transfer problems, comparisons between the linear FEM and the ziti are also included in this section. We present the approximated solution given by the finite element method. In this direction, a large body of works was already done; see [20]. The solution of the heat transfer using the finite element method at is illustrated in Figure 16.
5.4. Detection of Blowup in the Heat Transfer with Bad Sign
In the current study, we are interested in the two-dimensional transient diffusion problem with bad sign (negative diffusion coefficient), governed by the following boundary-value equation:where is a given time interval and is the unit disk with boundary . Here, denotes the space variables, is the time variable, is a negative diffusion coefficient, and represents the effect of internal source terms, given by the following equality:
We can verify that the exact solution is given bywhich converges to Dirac in the weak sense when goes to 0 as shown in Figure 17. The objective is to follow the behavior of the approximated solution, obtained by the same scheme presented in (47)–(48), which is illustrated in Figure 18.
6. Conclusion
As a conclusion, the current study is an extension of the ziti method when the domain of consideration is a disk, to be able to generalize it in the case of a multidimensional ball in . We started by presenting the fundamental tools of construction using two strategies: the Cartesian case, which is based on the hypothesis of being able to sweep the inside of the domain by a set of intervals, horizontally and vertically; therefore, all the work resides in the construction of the nodes in every interval. The second strategy is the polar one, using the polar coordinates in all the steps of construction. The first application in this work is the generalized integration; we obtained a good approximation even if some of the classical methods did not give any result. We applied our method for two types of . Where the exact solution is regular, we also compared our results with the finite element method, and the result was impressive. On the other hand, the right-hand side is obtained such that the exact solution presents a singularity of type Dirac, where the parameter of regularization is near 0 and less than the step of discretization. For the stationary diffusion problem, we compare the exact solution which is a very sensitive function approximating the Dirac mass, with the approximated one which goes also to Dirac. The numerical problem faced in this case is how to find a numerical scheme which detects this type of singularity and gives a solution having the same regular behavior of the exact one: our method can achieve this objective. For the heat transfer, we treated two types of singularities: the negative diffusion parameter (bad sign) and the case when the initial data is in function of to test the ability to detect such types of singularity. The case of a negative diffusion presents a great difficulty numerically, since most of the numerical schemes tend to explode, even if the domain is Cartesian and therefore the solution is lost. The Improved CFL condition saves time for advanced previsions (CFL around 0.9, which requires a fewer number of iterations to reach the final time). We can conclude that, with the extension and development of ziti, we have reached our goal: to approach differential problem defined in a disk, to find the results of the classical methods in the best conditions (speed, flexibility), even in the presence of highly singular terms (the numerical scheme detect the singularity, and solution continues regularly). Note that the maximum errors are calculated, finding that the errors are small. Also, a comparison of numerical and analytical solutions is made, finding that the proposed scheme has good accuracy.
Data Availability
No data were used to support this study.
Conflicts of Interest
The authors declare that they have no conflicts of interest.