# A New Numerical Method to Solve SomePDEs in the Unit Ball and Comparison with the Finite Element and the Exact Solution

Authors: Rajae Malek; Chérif Ziti
Journal: International Journal of Differential Equations
 (2021)

---

## 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.


---

## Body

## 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 ℝd 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 PDEs, 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 p-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 δ−ziti 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 (PDEs and ODEs). 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 by(1)Φx=exp1x2−1,ifx<1,0,otherwise,where x∈ℝd and d 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.Figure 1
Illustration of the elements of the familyφi.The outline of this work is as follows: in Section2, 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 δ−ziti 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 δ−ziti method with a given analytical one defined in a disk domain and to calculate the error in L∞Ω. 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 PDEs 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φii

From the functionΦ defined by (1), we define the set φε:(2)φεx=CεΦxε,for allε>0,where C\coloneq 1/∫ℝΦxdx is a constant of normalization. This sequence φε converges to Dirac measure δ0 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 h=b−a/m where m is a given integer such that xi=a+i−1h,∀i∈1,m+1.We construct the familyφi1≤i≤m+1 as follows:(3)φix=φhx−xi=ChΦx−xih,for allx∈xi−1,xi+1,i∈2,…,m,φ1x=φhx−x1=ChΦx−x1h,for allx∈x1,x2,φm+1x=φhx−xm+1=ChΦx−xm+1h,for allx∈xm,xm+1.The curves ofφi are shown in Figure 1.


### 2.2. Construction of the Orthonormal FamilyΨii

Let us consider the Hilbert spaceL2ℝ with the usual scalar product ,. Observe that the family φi1≤i≤m+1 is linearly independent, then using the Gram–Schmidt process, we construct a unique orthogonal family, noted as Ψ˜i satisfying the following relation:(4)Ψ˜ix=φix+λi−1Ψ˜i−1x,λ1=−αβ,λi+1=gλi=λ12−λ1λi,i=1,…,m−1,with α=φ1,φ2,β=φ1,φ1. Using the orthogonality of the family Ψ˜ii, we can find the second definition of λi given by(5)λi=−φi,Ψ˜i−1Ψ˜i−1,Ψ˜i−1,and the recurrence application of definition (4) gives the following formula:(6)Ψ˜ix=φix+λi−1φi−1x+λi−1λi−2φi−2x+⋯+∏k=1i−1λkφ1.LetΨi=Ψ˜i/Ψ˜i the normalization of Ψ˜i. The curves of Ψi are shown in Figure 2.Figure 2
Illustration of orthonormal familyΨi.We can verify the following assumption.Fori=2,…,m+1, Ψi admits i−1 reel root(s) noted as rkk=1,…,i−1 in the interval x1,xi, more precisely, rk∈xk,xk+1.The method permits approaching a given functionf and its integral by the following relations:(7a)fx≃∑i=1m+1ciΨix,(7b)ci≃∫abfxΨixdx,(7c)∫abfxdx≃∑i=1m+1ciIi,whereIi≔∫abΨixdx.The main idea of theδ−ziti method is to use the roots rk of the functions Ψi which satisfy(8)Ψiri>0,Ψirk=0,for alli≠k.If we takex=rk in (7a), we obtain(9a)ci≃friΨiri,i∈1,…,m,(9b)cm+1≃fbΨm+1b,(9c)∫abfxΨixdx≃friΨiri,(9d)∫abfxdx≃∑i=1mfriΨi2ri+fbΨm+1b2.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:(10)λi+1−λi<εas soon asi≥N0=lnε2−λ12/λ13−λ1lnλ1/2+λ12+1,where [.] denotes the integer part. In particular, for ε=10−M, we conclude that the parameters λi are nearly stationary from a certain rank, which reduces considerably the number of iterations. Using ri as a root of Ψi+1, we can define the parameter λi by λi=−φi+1ri/φiri (for more details of the construction of the δ−ziti method, see [1, 2].)


## 2.1. Construction of the Linearly Independent Familyφii

From the functionΦ defined by (1), we define the set φε:(2)φεx=CεΦxε,for allε>0,where C\coloneq 1/∫ℝΦxdx is a constant of normalization. This sequence φε converges to Dirac measure δ0 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 h=b−a/m where m is a given integer such that xi=a+i−1h,∀i∈1,m+1.We construct the familyφi1≤i≤m+1 as follows:(3)φix=φhx−xi=ChΦx−xih,for allx∈xi−1,xi+1,i∈2,…,m,φ1x=φhx−x1=ChΦx−x1h,for allx∈x1,x2,φm+1x=φhx−xm+1=ChΦx−xm+1h,for allx∈xm,xm+1.The curves ofφi are shown in Figure 1.


## 2.2. Construction of the Orthonormal FamilyΨii

Let us consider the Hilbert spaceL2ℝ with the usual scalar product ,. Observe that the family φi1≤i≤m+1 is linearly independent, then using the Gram–Schmidt process, we construct a unique orthogonal family, noted as Ψ˜i satisfying the following relation:(4)Ψ˜ix=φix+λi−1Ψ˜i−1x,λ1=−αβ,λi+1=gλi=λ12−λ1λi,i=1,…,m−1,with α=φ1,φ2,β=φ1,φ1. Using the orthogonality of the family Ψ˜ii, we can find the second definition of λi given by(5)λi=−φi,Ψ˜i−1Ψ˜i−1,Ψ˜i−1,and the recurrence application of definition (4) gives the following formula:(6)Ψ˜ix=φix+λi−1φi−1x+λi−1λi−2φi−2x+⋯+∏k=1i−1λkφ1.LetΨi=Ψ˜i/Ψ˜i the normalization of Ψ˜i. The curves of Ψi are shown in Figure 2.Figure 2
Illustration of orthonormal familyΨi.We can verify the following assumption.Fori=2,…,m+1, Ψi admits i−1 reel root(s) noted as rkk=1,…,i−1 in the interval x1,xi, more precisely, rk∈xk,xk+1.The method permits approaching a given functionf and its integral by the following relations:(7a)fx≃∑i=1m+1ciΨix,(7b)ci≃∫abfxΨixdx,(7c)∫abfxdx≃∑i=1m+1ciIi,whereIi≔∫abΨixdx.The main idea of theδ−ziti method is to use the roots rk of the functions Ψi which satisfy(8)Ψiri>0,Ψirk=0,for alli≠k.If we takex=rk in (7a), we obtain(9a)ci≃friΨiri,i∈1,…,m,(9b)cm+1≃fbΨm+1b,(9c)∫abfxΨixdx≃friΨiri,(9d)∫abfxdx≃∑i=1mfriΨi2ri+fbΨm+1b2.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:(10)λi+1−λi<εas soon asi≥N0=lnε2−λ12/λ13−λ1lnλ1/2+λ12+1,where [.] denotes the integer part. In particular, for ε=10−M, we conclude that the parameters λi are nearly stationary from a certain rank, which reduces considerably the number of iterations. Using ri as a root of Ψi+1, we can define the parameter λi by λi=−φi+1ri/φiri (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δ−ziti method, when Ω is a disk centered in the origin O=0,0 (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 Ψi); after that, we compute the roots (denoted by ri) of every element of this basis (Ψirk=0,for alli≠k; see Figure 5).Figure 3
Example of horizontal and vertical segments and the repartition of internal nodes.Figure 4
Sweeping the interior of the domain by a set of horizontal and vertical segments and position of the nodes.Figure 5
Illustration of the basis functionΨ4 with 3 roots, localized in the segment [−1, 1].Remark 1.
Note that, for every fixed vertical levelj (resp., the horizontal level i), every internal segment is limited by aj=−1−yj2 and bj=1−yj2 (resp., ci=−1−xi2 and di=1−xi2; see Figure 3). For simplicity, the step of the horizontal subdivision will be noted hj (resp., the vertical step will be noted hi).
We present an algorithm to calculate the internal nodes (Algorithm1).Algorithm 1:Construction of the nodes in the Cartesian case.
Data:N = the number of nodes, in every intervalFix an interval [a, b]h1=b−a/Nfori=1,…,Ndoxi1=a+i−1h1 horizontal nodes for the first interval.yi1=a+i−1h1 vertical nodes for the first intervalendforj=2,…,Ndofori=1,…,Ndoa = −1−−1+j⋅h2b = −ahj=b−a/Nxij=a+i−1hjendend


### 3.1. Construction of the Orthonormal Set

For every nodexij (resp., yji), we associate the function φij (noted as φi if there is no ambiguity) (resp., the family φji will be noted as φj) defined by(11)φix≔chjΦx−xijhj,∀i,j=1,…,N,φjy≔chiΦy−yjihi,∀j,i=1,…,N,wherehj is the step of construction in the horizontal interval of indication j, which describes the distance between the nodes xij and xi+1jhi is the step of subdivision in the vertical interval of indication i, which describes the distance between yji and yj+1iIt is simple to see that the familyφi is linearly independent, so we can construct an orthogonal family Ψ˜ii=1,…,N by using the Gram–Schmidt process, in the space L2a,b (construction in every internal interval of the domain Ω=B0,1, horizontally and vertically), verifying the following relation:(12)horizontally:Ψ˜1x=φ1x,Ψ˜ix=φix+∑k=1i−1λkiΨ˜kx,for alli=2,…,N,vertically:Ψ˜1y=φ1y,Ψ˜iy=φiy+∑k=1i−1λkiΨ˜ky,for alli=2,…,N,which will be reduced in the following theorem, already proved in the monodimensional case (see [1, 2]).Theorem 1.
The orthogonal familyΨ˜ii=1,…,N (vertically and horizontally) verifies the following recurrence relation:(13)Ψ˜1=φ1,Ψ˜i+1=φi+1+λiΨ˜i,for alli=1,…,N−1,λi=−φi,φi+1Ψ˜i,Ψ˜i,where , is the usual scalar product in the Hilbert space L2a,b.Corollary 1.
The familyΨ˜i and the set λi defined in Theorem 1 verify the following relations:(14)1Ψ˜i=φi+λi−1φi−1+λi−1λi−2φi−2+⋯+λi−1⋯λ1φ1,2Ψ˜ixij=φixij=chj2e,3Ψ˜iyij=φiyij=chi2e,4In every fixed level,φi,Ψ˜i−1=φi,φi−1.5−1<λi=−φi,φi+1Ψ˜i,Ψ˜i<0.


### 3.2. Fundamental Results. Numerical Integrations

In this paragraph, we are interested in the approximation of integrals, where the domain is the unit diskΩ=B0,1, using the horizontal and vertical test functions, as well as the roots, verifying the following relations:(15)Ψijx,y=Ψijx⋅Ψjiy,∀i,j=1,…,N,rij=rij,sji,whereΨijx are the basis functions in the horizontal dimension (resp., Ψjiy are the basis functions in the vertical dimension)rij are the roots of Ψijx (resp., sji are the roots of Ψjiy)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Ω=B0,1, let g be a given function in L2Ω, and let N denote the roots number; therefore, we have the following approximations:(16)∫Ωgx,ydxdy≃∑i,j=1Ngrij,sjiΨijrij⋅Ψjisji∫ajbjΨijxdx∫cidiΨjiydy,∫Ωgx,yΨijx,ydxdy≃grij,sjiΨijrij⋅Ψjisji2,where we take rNj=bj and siN=bi.Proof.
In [1, 2], the authors approximated integral formulas in the one-dimensional case using the δ−ziti method as follows:(17)∫abfxΨixdx≃friΨiri,∫abfxdx≃∑i=1mfriΨi2ri+fbΨm+1b2.
We recall the principal points of the proof which is inspired by the spectral method:(18)fx≃∑i=1m+1αiΨix,with αi=∫abfxΨixdx. Using the roots rk of Ψk, we obtain(19)frk≃∑i=1m+1αiΨirk=αkΨkrk,because of the following orthogonality result:(20)Ψirk=0,ifi≠k,Ψkrk,ifi=k.
Multiplying (18) by Ψk and integrating it over a,b, we obtain(21)∫abfxΨkxdx≃∑i=1m+1friψiri∫abψixψkxdx.
Using the orthonormal property ofΨi, (21) becomes(22)∫abfxψkxdx≃∑i=1m+1friψiriδik,and, therefore,(23)∫abfxψkxdx≃frkψkrk.
In the bidimensional case, using the Cartesian strategy,(24)Ψijx,y≔ΨijxΨjiy,so if we take fx,y≃∑i,j=1NαijΨijx,y, where N is the nodes number in every segment, horizontally and vertically, we will have the following approximation formulas:(25)∫Ωfx,ydxdy≃∑i,j=1Nfrij,sjiΨijrij⋅Ψjisji∫ajbjΨijxdx∫cidiΨjiydy,∫Ωfx,yΨijx,ydxdy≃frij,sjiΨijrij⋅Ψjisji2.
Here, we take,rNj=bj and siN=bi.
In Table1, we present some numerical tests of integration. We compare the exact value noted Ex, with the numerical approximation obtained by the δ−ziti method in the Cartesian case. Each interval is subdivided into 50 nodes.Table 1
Comparison between numerical integration usingδ−ziti and the exact value in the Cartesian case.
∫Ωfx,ydxdyExδ−ZitiError∫Ω1/x2+y2+11/4dxdy2.8558871302.8558348910.000052239∫Ω1/x2+y2+1dxdy2.6025805692.6020557810.000524788∫Ωexp1/x2+y2+2dxdy4.7177508864.7171385200.000612366∫Ωln1/x2+y2+2dxdy−2.857412051−2.8573464430.000065608


## 3.1. Construction of the Orthonormal Set

For every nodexij (resp., yji), we associate the function φij (noted as φi if there is no ambiguity) (resp., the family φji will be noted as φj) defined by(11)φix≔chjΦx−xijhj,∀i,j=1,…,N,φjy≔chiΦy−yjihi,∀j,i=1,…,N,wherehj is the step of construction in the horizontal interval of indication j, which describes the distance between the nodes xij and xi+1jhi is the step of subdivision in the vertical interval of indication i, which describes the distance between yji and yj+1iIt is simple to see that the familyφi is linearly independent, so we can construct an orthogonal family Ψ˜ii=1,…,N by using the Gram–Schmidt process, in the space L2a,b (construction in every internal interval of the domain Ω=B0,1, horizontally and vertically), verifying the following relation:(12)horizontally:Ψ˜1x=φ1x,Ψ˜ix=φix+∑k=1i−1λkiΨ˜kx,for alli=2,…,N,vertically:Ψ˜1y=φ1y,Ψ˜iy=φiy+∑k=1i−1λkiΨ˜ky,for alli=2,…,N,which will be reduced in the following theorem, already proved in the monodimensional case (see [1, 2]).Theorem 1.
The orthogonal familyΨ˜ii=1,…,N (vertically and horizontally) verifies the following recurrence relation:(13)Ψ˜1=φ1,Ψ˜i+1=φi+1+λiΨ˜i,for alli=1,…,N−1,λi=−φi,φi+1Ψ˜i,Ψ˜i,where , is the usual scalar product in the Hilbert space L2a,b.Corollary 1.
The familyΨ˜i and the set λi defined in Theorem 1 verify the following relations:(14)1Ψ˜i=φi+λi−1φi−1+λi−1λi−2φi−2+⋯+λi−1⋯λ1φ1,2Ψ˜ixij=φixij=chj2e,3Ψ˜iyij=φiyij=chi2e,4In every fixed level,φi,Ψ˜i−1=φi,φi−1.5−1<λi=−φi,φi+1Ψ˜i,Ψ˜i<0.


## 3.2. Fundamental Results. Numerical Integrations

In this paragraph, we are interested in the approximation of integrals, where the domain is the unit diskΩ=B0,1, using the horizontal and vertical test functions, as well as the roots, verifying the following relations:(15)Ψijx,y=Ψijx⋅Ψjiy,∀i,j=1,…,N,rij=rij,sji,whereΨijx are the basis functions in the horizontal dimension (resp., Ψjiy are the basis functions in the vertical dimension)rij are the roots of Ψijx (resp., sji are the roots of Ψjiy)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Ω=B0,1, let g be a given function in L2Ω, and let N denote the roots number; therefore, we have the following approximations:(16)∫Ωgx,ydxdy≃∑i,j=1Ngrij,sjiΨijrij⋅Ψjisji∫ajbjΨijxdx∫cidiΨjiydy,∫Ωgx,yΨijx,ydxdy≃grij,sjiΨijrij⋅Ψjisji2,where we take rNj=bj and siN=bi.Proof.
In [1, 2], the authors approximated integral formulas in the one-dimensional case using the δ−ziti method as follows:(17)∫abfxΨixdx≃friΨiri,∫abfxdx≃∑i=1mfriΨi2ri+fbΨm+1b2.
We recall the principal points of the proof which is inspired by the spectral method:(18)fx≃∑i=1m+1αiΨix,with αi=∫abfxΨixdx. Using the roots rk of Ψk, we obtain(19)frk≃∑i=1m+1αiΨirk=αkΨkrk,because of the following orthogonality result:(20)Ψirk=0,ifi≠k,Ψkrk,ifi=k.
Multiplying (18) by Ψk and integrating it over a,b, we obtain(21)∫abfxΨkxdx≃∑i=1m+1friψiri∫abψixψkxdx.
Using the orthonormal property ofΨi, (21) becomes(22)∫abfxψkxdx≃∑i=1m+1friψiriδik,and, therefore,(23)∫abfxψkxdx≃frkψkrk.
In the bidimensional case, using the Cartesian strategy,(24)Ψijx,y≔ΨijxΨjiy,so if we take fx,y≃∑i,j=1NαijΨijx,y, where N is the nodes number in every segment, horizontally and vertically, we will have the following approximation formulas:(25)∫Ωfx,ydxdy≃∑i,j=1Nfrij,sjiΨijrij⋅Ψjisji∫ajbjΨijxdx∫cidiΨjiydy,∫Ωfx,yΨijx,ydxdy≃frij,sjiΨijrij⋅Ψjisji2.
Here, we take,rNj=bj and siN=bi.
In Table1, we present some numerical tests of integration. We compare the exact value noted Ex, with the numerical approximation obtained by the δ−ziti method in the Cartesian case. Each interval is subdivided into 50 nodes.Table 1
Comparison between numerical integration usingδ−ziti and the exact value in the Cartesian case.
∫Ωfx,ydxdyExδ−ZitiError∫Ω1/x2+y2+11/4dxdy2.8558871302.8558348910.000052239∫Ω1/x2+y2+1dxdy2.6025805692.6020557810.000524788∫Ωexp1/x2+y2+2dxdy4.7177508864.7171385200.000612366∫Ωln1/x2+y2+2dxdy−2.857412051−2.8573464430.000065608


## 4. Second Strategy: Polar Coordinates

In this section, we built all the necessary elements for the approximationδ−ziti, using polar coordinates. The domain Ω=B0,1 is represented using the polar coordinates, with the following parametrization:(26)∀x,y∈B0,1,x=ρcosθ,y=ρsinθ,ρ,θ∈0,1×0,2π.The polar setΨijρ,θ is defined by(27)Ψijρ,θ≔Ψiρ.Ψjθ.


### 4.1. Fundamental Results: Numerical Integration

To test the previous strategy, we present in Tables2 and 3 some numerical tests. We compare the exact value with the numerical one, using the δ−ziti method in the polar case.Table 2
Comparison between numerical integration using theδ−ziti and the exact value in the polar case.
1/2π∫Ωfx,ydxdyExδ−zitiError∫Ω1/x2+y2+11/4dxdy0.45452855370.4544964599186500.0000320937813496069∫Ω1/x2+y2+1dxdy0.41421356240.4144406924675260.0002271300675258380∫Ωexp1/x2+y2+2dxdy0.75085337380.7507720373200430.0000813364799567839∫Ωln1/x2+y2+2dxdy0.4547712524−0.4536232978390540.0011479545609460200∫Ωxydxdy00.0001048932849246290.0001∫Ωlnx2+y2/x2+y2dxdy−1−1.000112995319610.000112995319609510Table 3
Generalized integrals using the polar parametrization.
∫Ωfx,ydxdyExδ−zitiError∫Ωy/x2+y2x2+y2+1arctany/x1/3dxdy0.22043663480.2233581847629060.00292154∫Ωy/x2+y2arctany/x1/2dxdy00.000796460.00079646∫Ω1/2πx2+y23/4dxdy1/π=0.31818181820.31657675840.0016057∫Ωlnx2+y2/2πx2+y2dxdy−0.500/2π=0.15909090910.15609090910.003In Table3, we can treat even the generalized integrals, for example. ∫Ωy/x2+y2arctany/x1/2dxdy, which is in fact an operation of Riemann integral; we found a good approximation using the δ−ziti roots. For the two last examples, ∫Ω1/2πx2+y23/4dxdy and ∫Ωlnx2+y2/2πx2+y2dxdy, other approximation methods (e.g., Simpson, Trapeze\“enleadertwodots”) did not give any result, which is an important point for our construction.


## 4.1. Fundamental Results: Numerical Integration

To test the previous strategy, we present in Tables2 and 3 some numerical tests. We compare the exact value with the numerical one, using the δ−ziti method in the polar case.Table 2
Comparison between numerical integration using theδ−ziti and the exact value in the polar case.
1/2π∫Ωfx,ydxdyExδ−zitiError∫Ω1/x2+y2+11/4dxdy0.45452855370.4544964599186500.0000320937813496069∫Ω1/x2+y2+1dxdy0.41421356240.4144406924675260.0002271300675258380∫Ωexp1/x2+y2+2dxdy0.75085337380.7507720373200430.0000813364799567839∫Ωln1/x2+y2+2dxdy0.4547712524−0.4536232978390540.0011479545609460200∫Ωxydxdy00.0001048932849246290.0001∫Ωlnx2+y2/x2+y2dxdy−1−1.000112995319610.000112995319609510Table 3
Generalized integrals using the polar parametrization.
∫Ωfx,ydxdyExδ−zitiError∫Ωy/x2+y2x2+y2+1arctany/x1/3dxdy0.22043663480.2233581847629060.00292154∫Ωy/x2+y2arctany/x1/2dxdy00.000796460.00079646∫Ω1/2πx2+y23/4dxdy1/π=0.31818181820.31657675840.0016057∫Ωlnx2+y2/2πx2+y2dxdy−0.500/2π=0.15909090910.15609090910.003In Table3, we can treat even the generalized integrals, for example. ∫Ωy/x2+y2arctany/x1/2dxdy, which is in fact an operation of Riemann integral; we found a good approximation using the δ−ziti roots. For the two last examples, ∫Ω1/2πx2+y23/4dxdy and ∫Ωlnx2+y2/2πx2+y2dxdy, 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Ω=B0,1. The problem studied is given by(28a)−divK.∇u=finΩ,(28b)ux,y=0in∂Ω,and, for simplicity of numerical simulations, we take K=I2, with a given analytical solution uex=1−x2−y2, where the source term takes the value f=4.


#### 5.1.2. The Strong Discretization

The first step to approach the previous problem is to multiply equation (7a) by a test function Ψij and to integrate the result over the domain Ω=B0,1, which gives(29)−∫ΩΔux,y⋅Ψijx,ydxdy=∫Ωfx,y⋅Ψijx,ydxdy.Using Theorem2, we obtain the following equality:(30)−Δurij,sjiΨijrij,sji=frij,sjiΨijrij,sji,∀rij,sji∈Ω.The next step consists of approaching the second derivative, which gives us the following scheme:(31)ui−1,j+2uij+ui+1,jri+1j−rijrij−ri−1j+ui,j−1+2uij+ui,j+1sj+1i−sjisji−sj−1i=fij,i,j=2,…,N−1,u1,j=uN,j=0,j=1,…,N,ui,1=ui,N=0,i=1,…,N,where uij denotes the approximation at the root rij,sji and N is the nodes number in every internal segment (horizontally and vertically). In the end, we will have a global matrix, with N−2×N−2 lines and N−2×N−2 columns, defined as follows:(32)M=D2A30⋯0A2D3A4⋯0⋮0⋯An−3Dn−2An−100⋯An−2Dn−1,where Di is a N−2×N−2 tridiagonal matrix, defined by(33)Dk,ki+1=−2Ψi,k+1dxk⋅dxk−1−2Ψi,k+1dyk⋅dyk−1,i,k=1,…,N−2,Dk,k+1i+1=−2Ψi,k+2dyk⋅dyk−1,i,k=1,…,N−2,Dk−1,ki+1=−2Ψi,kdyk⋅dyk−1,i,k=1,…,N−2,and Ai is a N−2×N−2 diagonal matrix defined by(34)Ak,ki=−Ψi,kdxk⋅dxk−1,i=3,…,N−1,k=2,…,N−1,where(35)dxk=rkj−rk−1j,dyk=rki−rk−1i,and, therefore, we should resolve a simple system in the form MX=F, when M is the global matrix defined previously, X is the unknown vector of size N−2×N−2, and F is the source term vector of size N−2×N−2.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Ω=B0,1. 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 δ−ziti scheme, using the Cartesian coordinates and the segments approach. For a fixed node’s number in every segment (horizontal or vertical), N=100, the numerical implementation of the scheme gives us an approximated solution, which is near the exact one, given by uexx,y=1−x2−y2, when fx,y=4. The two solutions are illustrated in Figures 6 and 7.Figure 6
The exact solution.Figure 7
Theδ−ziti solution.In Table4, we present the error between the exact and approximated solution, with different values of nodes number N.In Table 4,(36)Ermax≔maxi,juexi,j−ui,j,Ermean≔1N∑i,juexi,j−ui,j.Table 4
The committed error using several values of nodes numberN.
NhminhmaxErmaxErmean600.0085340.033330.011740.007581000.0039790.020000.00460.002881500.002170.013330.002040.00131692000.0014100.010000.001670.0007514823We 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 v 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.Figure 8
Approximated solution using finite element method.


#### 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δ−ziti method in the polar case.LetΩ=B0,1; the strategy presented consists in using the results of approximation in the monodimensional case and taking into consideration the following function basis:(37)Ψijρ,θ=ΨiρΨjθ,∀i,j=1,…,N.The problem presented in the previous Subsection5.1.1 is equivalent to the polar one, expressed as follows:(38a)∂2u∂ρ2+1ρ∂u∂r+1ρ2∂2u∂θ2=finΩ,(38b)uρ=1,θ=0,∀θ∈0,2π,with(39)uρ,θ≔uρcosθ,ρsinθ,fρ,θ≔fρcosθ,ρsinθ=4,where f is the source term defined in Cartesian problem, with an exact polar solution uex=1−ρ2. Note that, the roots of the basic functions Ψiρ will be noted; ri and θj are those associated with Ψjθ.To obtain a numerical scheme using theδ−ziti method, we should multiply equation (38a) by a test function Ψijρ,θ and after that, we use the strong result of approximation (16), which gives(40)Δiju=fij,i,j=2,…,N−1,u1,j=u2,j,j=1,…,N,uM,j=0,j=1,…,N,ui,1=ui,2,i=1,…,N,ui,M=ui,M−1,i=1,…,N,with(41)Δiju=ui−1,j−2uij+ui+1,jri+1−riri−ri−1+1riui+1,j−uijri+1−ri+1ri2ui,j−1+2uij+ui,j+1θj+1−θjθj−θj−1.Therefore, the goal is to finduρ,θ 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 MX=F; when M is the global matrix defined previously, we should just add the polar terms 1/ρ and 1/ρ2 in the corresponding terms of the matrices Ai and Di, X in the polar unknown vector of size N−2×N−2, and F is the source term vector of size N−2×N−2.


#### 5.1.6. Numerical Tests

LetNr=100 be the root’s number for the radius and Nθ=100 the other one for the angles. Two solutions, exact and approximated, are illustrated in Figures 9 and 10 .Figure 9
The exact solution.Figure 10
Theδ−ziti solution.


#### 5.1.7. Numerical Error

Table5 shows us the error between exact and approximated solution, using different Nr and Nθ. In Figure 11, we present the evolution of the error between exact and approximated solution, with several values of the nodes number Nr for the diffusion problem.Table 5
The infinite error, using several values ofNr and Nθ.
NrNθhrhθError60600.016660.0333π8.86.10−41001000.010.02π3.82.10−41501500.006660.0133π1.91.10−42002000.0050.01π1.15.10−4Figure 11
Evolution of the error with several values ofNr. (a) The error for Nr=60. (b) The error for Nr=100. (c) The error for Nr=150. (d) The error for Nr=200.
(a)(b)(c)(d)Table5 shows how the error is reduced when the number of nodes increases from Nr=60 to Nr=200 (resp. Nθ) 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:(42)−div∇u=finΩ,ux,y=0in∂Ω,where the source term is given by(43)fx,y=−2C−2x4−4x2y2+2ε4−2ε2y2−2y4−2ε2x2ε2−x2−y24Φxε,yε,where C is the normalization constant defined by (1) and Φ is the regular function given by(44)Φx,y=exp1x2+y2−1,ifx2+y2<1,0,otherwise.It is easy to see that the analytical solution of (42) and (43) is given by(45)uexx,y=Cε2Φxε,yε.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.Figure 12
The exact solution withε=0.01, and maxuex=2230.Figure 13
Theδ−ziti solution with ε=0.01, and maxuap=2450.


### 5.3. Parabolic PDE Case: Heat Equation

This section is devoted to the application of theδ−ziti method on a diffusion equation, in a domain Σ=0,T×Ω, where Ω=B0,1. The heat equation describes the distribution of heat (or variation in temperature) in a given region over time. For a function, ut,x,y (resp., ut,r,θ) of two spatial variables x,y in the Cartesian case (r,θ in the polar case) and the time variable t, the heat equation is given by(46a)∂u∂t−DΔu=ft,x,t,x∈Σ,(46b)ut,x=0,x∈∂Ω,(46c)u0,x=u0x≥0,x∈Ω,where0,T is a given time interval and Ω=B0,1 is the unit disk with boundary ∂Ω. Here, x=x,yt denotes the space variables, t is the time variable, D is the diffusion coefficient, and f represents the effect of internal source terms. Using the same analogy applied in the previous sections, we multiply equation (46a) by a test function Ψij, after we integrate over the domain Ω. It remains just the direct application of our approximations formulas given in (3.3). First, the time domain 0,tmax is divided into Nt uniform subintervals tn,tn+1 with duration dt=tn+1−tn for n=0,…,Nt. To denote the value of u in the root rij at time tn, we use the notation ui,jn for i,j=1,…,N. At t=tn, we find directly the following explicit scheme:(47)uijn+1=uijn+dt⋅D⋅Δijnu+dt⋅fi,jn,i,j=2,…,N−1,u1,jn+1=0,j=1,…,N,uN,jn+1=0,j=1,…,N,ui,1n+1=0,i=1,…,N,ui,Nn+1=0,i=1,…,N,where(48)Δijn=ui−1,jn−2uijn+ui+1,jnri+1−riri−ri−1+ui,j−1n+2uijn+ui,j+1nθj+1−θjθj−θj−1.Our goal is to find an approximated solution, near the exact one, which verifies the boundary conditions. The function(49)uext,x,y=1−x2−y2exptis an exact solution of the heat equation, with(50)D=1,ft,x,y=−3−x2−y2expt.The above example provides a starting point for numerically solving diffusion problems on the disk. These ideas can also be used on the sphere. Figures14 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.Figure 14
The exact solution att=0.01, using the Cartesian discretization.Figure 15
Theδ−ziti solution at t=0.01, after 345-time iteration, using the Cartesian discretization.In Table6, 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 t=0.01 (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 N=60 to N=200.Table 6
The committed error att=0.01.
NhminhmaxErmaxErmeanCPU600.0085340.033332.22431.10−51.413039.10−5571000.0039790.020001.0174501.10−65.453807.10−7831500.002170.013331.00359253.10−65.36991485.10−71222000.0014100.010009.9973690.10−75.3444320.10−7150For checking the stability of our numerical scheme, we use the following strategy:(51)Edtn=un:,:−un−1:,:∞un:,:∞,in which un:,: is the matrix of numerical solution at nth iteration. Errors obtained corresponding to Edtn using the δ−ziti method with N=200 for the heat transfer problem are reported in Table 7.Table 7
Errors obtained for the present method withN=200 at different final times for the heat transfer problem.
IterationsEdtndt=1×10−7dt=0.5×10−7dt=1×10−8dt=0.5×10−854.4282×10−62.2152×10−64.4335×10−72.2171×10−7204.4284×10−62.2155×10−64.4338×10−72.2173×10−7304.4286×10−62.2159×10−64.4340×10−72.2177×10−71124.4298×10−62.2163×10−64.4351×10−72.2189×10−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 t=1 is illustrated in Figure 16.Figure 16
Approximated solution using the finite element method.


### 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:(52a)∂u∂t−DΔu=ft,x,t,x∈0,T×Ω,(52b)ut,x=0,x∈∂Ω,(52c)u0,x=u0x≥0,x∈Ω,where 0,T is a given time interval and Ω=B0,1 is the unit disk with boundary ∂Ω. Here, x=x,yt denotes the space variables, t is the time variable, D is a negative diffusion coefficient, and f represents the effect of internal source terms, given by the following equality:(53)ft,x,y=C1ε2−2D−2x4−4x2y2+2ε4−2ε2y2−2y4−2ε2x2ε2−x2−y24Φxε,yεexpt.We can verify that the exact solution is given by(54)uext,x,y=Cε2Φxε,yεexpt,which 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.Figure 17
The exact solution att = 0.01, with D = −0.25, ε=0.01, maxuext,.=710.Figure 18
Theδ−ziti solution at t = 0.01, with D = −0.25, ε=0.01maxuapt,.=780, after 345-time iterations, using the Cartesian discretization.


## 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Ω=B0,1. The problem studied is given by(28a)−divK.∇u=finΩ,(28b)ux,y=0in∂Ω,and, for simplicity of numerical simulations, we take K=I2, with a given analytical solution uex=1−x2−y2, where the source term takes the value f=4.


### 5.1.2. The Strong Discretization

The first step to approach the previous problem is to multiply equation (7a) by a test function Ψij and to integrate the result over the domain Ω=B0,1, which gives(29)−∫ΩΔux,y⋅Ψijx,ydxdy=∫Ωfx,y⋅Ψijx,ydxdy.Using Theorem2, we obtain the following equality:(30)−Δurij,sjiΨijrij,sji=frij,sjiΨijrij,sji,∀rij,sji∈Ω.The next step consists of approaching the second derivative, which gives us the following scheme:(31)ui−1,j+2uij+ui+1,jri+1j−rijrij−ri−1j+ui,j−1+2uij+ui,j+1sj+1i−sjisji−sj−1i=fij,i,j=2,…,N−1,u1,j=uN,j=0,j=1,…,N,ui,1=ui,N=0,i=1,…,N,where uij denotes the approximation at the root rij,sji and N is the nodes number in every internal segment (horizontally and vertically). In the end, we will have a global matrix, with N−2×N−2 lines and N−2×N−2 columns, defined as follows:(32)M=D2A30⋯0A2D3A4⋯0⋮0⋯An−3Dn−2An−100⋯An−2Dn−1,where Di is a N−2×N−2 tridiagonal matrix, defined by(33)Dk,ki+1=−2Ψi,k+1dxk⋅dxk−1−2Ψi,k+1dyk⋅dyk−1,i,k=1,…,N−2,Dk,k+1i+1=−2Ψi,k+2dyk⋅dyk−1,i,k=1,…,N−2,Dk−1,ki+1=−2Ψi,kdyk⋅dyk−1,i,k=1,…,N−2,and Ai is a N−2×N−2 diagonal matrix defined by(34)Ak,ki=−Ψi,kdxk⋅dxk−1,i=3,…,N−1,k=2,…,N−1,where(35)dxk=rkj−rk−1j,dyk=rki−rk−1i,and, therefore, we should resolve a simple system in the form MX=F, when M is the global matrix defined previously, X is the unknown vector of size N−2×N−2, and F is the source term vector of size N−2×N−2.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Ω=B0,1. 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 δ−ziti scheme, using the Cartesian coordinates and the segments approach. For a fixed node’s number in every segment (horizontal or vertical), N=100, the numerical implementation of the scheme gives us an approximated solution, which is near the exact one, given by uexx,y=1−x2−y2, when fx,y=4. The two solutions are illustrated in Figures 6 and 7.Figure 6
The exact solution.Figure 7
Theδ−ziti solution.In Table4, we present the error between the exact and approximated solution, with different values of nodes number N.In Table 4,(36)Ermax≔maxi,juexi,j−ui,j,Ermean≔1N∑i,juexi,j−ui,j.Table 4
The committed error using several values of nodes numberN.
NhminhmaxErmaxErmean600.0085340.033330.011740.007581000.0039790.020000.00460.002881500.002170.013330.002040.00131692000.0014100.010000.001670.0007514823We 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 v 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.Figure 8
Approximated solution using finite element method.


### 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δ−ziti method in the polar case.LetΩ=B0,1; the strategy presented consists in using the results of approximation in the monodimensional case and taking into consideration the following function basis:(37)Ψijρ,θ=ΨiρΨjθ,∀i,j=1,…,N.The problem presented in the previous Subsection5.1.1 is equivalent to the polar one, expressed as follows:(38a)∂2u∂ρ2+1ρ∂u∂r+1ρ2∂2u∂θ2=finΩ,(38b)uρ=1,θ=0,∀θ∈0,2π,with(39)uρ,θ≔uρcosθ,ρsinθ,fρ,θ≔fρcosθ,ρsinθ=4,where f is the source term defined in Cartesian problem, with an exact polar solution uex=1−ρ2. Note that, the roots of the basic functions Ψiρ will be noted; ri and θj are those associated with Ψjθ.To obtain a numerical scheme using theδ−ziti method, we should multiply equation (38a) by a test function Ψijρ,θ and after that, we use the strong result of approximation (16), which gives(40)Δiju=fij,i,j=2,…,N−1,u1,j=u2,j,j=1,…,N,uM,j=0,j=1,…,N,ui,1=ui,2,i=1,…,N,ui,M=ui,M−1,i=1,…,N,with(41)Δiju=ui−1,j−2uij+ui+1,jri+1−riri−ri−1+1riui+1,j−uijri+1−ri+1ri2ui,j−1+2uij+ui,j+1θj+1−θjθj−θj−1.Therefore, the goal is to finduρ,θ 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 MX=F; when M is the global matrix defined previously, we should just add the polar terms 1/ρ and 1/ρ2 in the corresponding terms of the matrices Ai and Di, X in the polar unknown vector of size N−2×N−2, and F is the source term vector of size N−2×N−2.


### 5.1.6. Numerical Tests

LetNr=100 be the root’s number for the radius and Nθ=100 the other one for the angles. Two solutions, exact and approximated, are illustrated in Figures 9 and 10 .Figure 9
The exact solution.Figure 10
Theδ−ziti solution.


### 5.1.7. Numerical Error

Table5 shows us the error between exact and approximated solution, using different Nr and Nθ. In Figure 11, we present the evolution of the error between exact and approximated solution, with several values of the nodes number Nr for the diffusion problem.Table 5
The infinite error, using several values ofNr and Nθ.
NrNθhrhθError60600.016660.0333π8.86.10−41001000.010.02π3.82.10−41501500.006660.0133π1.91.10−42002000.0050.01π1.15.10−4Figure 11
Evolution of the error with several values ofNr. (a) The error for Nr=60. (b) The error for Nr=100. (c) The error for Nr=150. (d) The error for Nr=200.
(a)(b)(c)(d)Table5 shows how the error is reduced when the number of nodes increases from Nr=60 to Nr=200 (resp. Nθ) 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.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Ω=B0,1. The problem studied is given by(28a)−divK.∇u=finΩ,(28b)ux,y=0in∂Ω,and, for simplicity of numerical simulations, we take K=I2, with a given analytical solution uex=1−x2−y2, where the source term takes the value f=4.


## 5.1.2. The Strong Discretization

The first step to approach the previous problem is to multiply equation (7a) by a test function Ψij and to integrate the result over the domain Ω=B0,1, which gives(29)−∫ΩΔux,y⋅Ψijx,ydxdy=∫Ωfx,y⋅Ψijx,ydxdy.Using Theorem2, we obtain the following equality:(30)−Δurij,sjiΨijrij,sji=frij,sjiΨijrij,sji,∀rij,sji∈Ω.The next step consists of approaching the second derivative, which gives us the following scheme:(31)ui−1,j+2uij+ui+1,jri+1j−rijrij−ri−1j+ui,j−1+2uij+ui,j+1sj+1i−sjisji−sj−1i=fij,i,j=2,…,N−1,u1,j=uN,j=0,j=1,…,N,ui,1=ui,N=0,i=1,…,N,where uij denotes the approximation at the root rij,sji and N is the nodes number in every internal segment (horizontally and vertically). In the end, we will have a global matrix, with N−2×N−2 lines and N−2×N−2 columns, defined as follows:(32)M=D2A30⋯0A2D3A4⋯0⋮0⋯An−3Dn−2An−100⋯An−2Dn−1,where Di is a N−2×N−2 tridiagonal matrix, defined by(33)Dk,ki+1=−2Ψi,k+1dxk⋅dxk−1−2Ψi,k+1dyk⋅dyk−1,i,k=1,…,N−2,Dk,k+1i+1=−2Ψi,k+2dyk⋅dyk−1,i,k=1,…,N−2,Dk−1,ki+1=−2Ψi,kdyk⋅dyk−1,i,k=1,…,N−2,and Ai is a N−2×N−2 diagonal matrix defined by(34)Ak,ki=−Ψi,kdxk⋅dxk−1,i=3,…,N−1,k=2,…,N−1,where(35)dxk=rkj−rk−1j,dyk=rki−rk−1i,and, therefore, we should resolve a simple system in the form MX=F, when M is the global matrix defined previously, X is the unknown vector of size N−2×N−2, and F is the source term vector of size N−2×N−2.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Ω=B0,1. 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 δ−ziti scheme, using the Cartesian coordinates and the segments approach. For a fixed node’s number in every segment (horizontal or vertical), N=100, the numerical implementation of the scheme gives us an approximated solution, which is near the exact one, given by uexx,y=1−x2−y2, when fx,y=4. The two solutions are illustrated in Figures 6 and 7.Figure 6
The exact solution.Figure 7
Theδ−ziti solution.In Table4, we present the error between the exact and approximated solution, with different values of nodes number N.In Table 4,(36)Ermax≔maxi,juexi,j−ui,j,Ermean≔1N∑i,juexi,j−ui,j.Table 4
The committed error using several values of nodes numberN.
NhminhmaxErmaxErmean600.0085340.033330.011740.007581000.0039790.020000.00460.002881500.002170.013330.002040.00131692000.0014100.010000.001670.0007514823We 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 v 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.Figure 8
Approximated solution using finite element method.


## 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δ−ziti method in the polar case.LetΩ=B0,1; the strategy presented consists in using the results of approximation in the monodimensional case and taking into consideration the following function basis:(37)Ψijρ,θ=ΨiρΨjθ,∀i,j=1,…,N.The problem presented in the previous Subsection5.1.1 is equivalent to the polar one, expressed as follows:(38a)∂2u∂ρ2+1ρ∂u∂r+1ρ2∂2u∂θ2=finΩ,(38b)uρ=1,θ=0,∀θ∈0,2π,with(39)uρ,θ≔uρcosθ,ρsinθ,fρ,θ≔fρcosθ,ρsinθ=4,where f is the source term defined in Cartesian problem, with an exact polar solution uex=1−ρ2. Note that, the roots of the basic functions Ψiρ will be noted; ri and θj are those associated with Ψjθ.To obtain a numerical scheme using theδ−ziti method, we should multiply equation (38a) by a test function Ψijρ,θ and after that, we use the strong result of approximation (16), which gives(40)Δiju=fij,i,j=2,…,N−1,u1,j=u2,j,j=1,…,N,uM,j=0,j=1,…,N,ui,1=ui,2,i=1,…,N,ui,M=ui,M−1,i=1,…,N,with(41)Δiju=ui−1,j−2uij+ui+1,jri+1−riri−ri−1+1riui+1,j−uijri+1−ri+1ri2ui,j−1+2uij+ui,j+1θj+1−θjθj−θj−1.Therefore, the goal is to finduρ,θ 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 MX=F; when M is the global matrix defined previously, we should just add the polar terms 1/ρ and 1/ρ2 in the corresponding terms of the matrices Ai and Di, X in the polar unknown vector of size N−2×N−2, and F is the source term vector of size N−2×N−2.


## 5.1.6. Numerical Tests

LetNr=100 be the root’s number for the radius and Nθ=100 the other one for the angles. Two solutions, exact and approximated, are illustrated in Figures 9 and 10 .Figure 9
The exact solution.Figure 10
Theδ−ziti solution.


## 5.1.7. Numerical Error

Table5 shows us the error between exact and approximated solution, using different Nr and Nθ. In Figure 11, we present the evolution of the error between exact and approximated solution, with several values of the nodes number Nr for the diffusion problem.Table 5
The infinite error, using several values ofNr and Nθ.
NrNθhrhθError60600.016660.0333π8.86.10−41001000.010.02π3.82.10−41501500.006660.0133π1.91.10−42002000.0050.01π1.15.10−4Figure 11
Evolution of the error with several values ofNr. (a) The error for Nr=60. (b) The error for Nr=100. (c) The error for Nr=150. (d) The error for Nr=200.
(a)(b)(c)(d)Table5 shows how the error is reduced when the number of nodes increases from Nr=60 to Nr=200 (resp. Nθ) 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:(42)−div∇u=finΩ,ux,y=0in∂Ω,where the source term is given by(43)fx,y=−2C−2x4−4x2y2+2ε4−2ε2y2−2y4−2ε2x2ε2−x2−y24Φxε,yε,where C is the normalization constant defined by (1) and Φ is the regular function given by(44)Φx,y=exp1x2+y2−1,ifx2+y2<1,0,otherwise.It is easy to see that the analytical solution of (42) and (43) is given by(45)uexx,y=Cε2Φxε,yε.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.Figure 12
The exact solution withε=0.01, and maxuex=2230.Figure 13
Theδ−ziti solution with ε=0.01, and maxuap=2450.


## 5.3. Parabolic PDE Case: Heat Equation

This section is devoted to the application of theδ−ziti method on a diffusion equation, in a domain Σ=0,T×Ω, where Ω=B0,1. The heat equation describes the distribution of heat (or variation in temperature) in a given region over time. For a function, ut,x,y (resp., ut,r,θ) of two spatial variables x,y in the Cartesian case (r,θ in the polar case) and the time variable t, the heat equation is given by(46a)∂u∂t−DΔu=ft,x,t,x∈Σ,(46b)ut,x=0,x∈∂Ω,(46c)u0,x=u0x≥0,x∈Ω,where0,T is a given time interval and Ω=B0,1 is the unit disk with boundary ∂Ω. Here, x=x,yt denotes the space variables, t is the time variable, D is the diffusion coefficient, and f represents the effect of internal source terms. Using the same analogy applied in the previous sections, we multiply equation (46a) by a test function Ψij, after we integrate over the domain Ω. It remains just the direct application of our approximations formulas given in (3.3). First, the time domain 0,tmax is divided into Nt uniform subintervals tn,tn+1 with duration dt=tn+1−tn for n=0,…,Nt. To denote the value of u in the root rij at time tn, we use the notation ui,jn for i,j=1,…,N. At t=tn, we find directly the following explicit scheme:(47)uijn+1=uijn+dt⋅D⋅Δijnu+dt⋅fi,jn,i,j=2,…,N−1,u1,jn+1=0,j=1,…,N,uN,jn+1=0,j=1,…,N,ui,1n+1=0,i=1,…,N,ui,Nn+1=0,i=1,…,N,where(48)Δijn=ui−1,jn−2uijn+ui+1,jnri+1−riri−ri−1+ui,j−1n+2uijn+ui,j+1nθj+1−θjθj−θj−1.Our goal is to find an approximated solution, near the exact one, which verifies the boundary conditions. The function(49)uext,x,y=1−x2−y2exptis an exact solution of the heat equation, with(50)D=1,ft,x,y=−3−x2−y2expt.The above example provides a starting point for numerically solving diffusion problems on the disk. These ideas can also be used on the sphere. Figures14 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.Figure 14
The exact solution att=0.01, using the Cartesian discretization.Figure 15
Theδ−ziti solution at t=0.01, after 345-time iteration, using the Cartesian discretization.In Table6, 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 t=0.01 (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 N=60 to N=200.Table 6
The committed error att=0.01.
NhminhmaxErmaxErmeanCPU600.0085340.033332.22431.10−51.413039.10−5571000.0039790.020001.0174501.10−65.453807.10−7831500.002170.013331.00359253.10−65.36991485.10−71222000.0014100.010009.9973690.10−75.3444320.10−7150For checking the stability of our numerical scheme, we use the following strategy:(51)Edtn=un:,:−un−1:,:∞un:,:∞,in which un:,: is the matrix of numerical solution at nth iteration. Errors obtained corresponding to Edtn using the δ−ziti method with N=200 for the heat transfer problem are reported in Table 7.Table 7
Errors obtained for the present method withN=200 at different final times for the heat transfer problem.
IterationsEdtndt=1×10−7dt=0.5×10−7dt=1×10−8dt=0.5×10−854.4282×10−62.2152×10−64.4335×10−72.2171×10−7204.4284×10−62.2155×10−64.4338×10−72.2173×10−7304.4286×10−62.2159×10−64.4340×10−72.2177×10−71124.4298×10−62.2163×10−64.4351×10−72.2189×10−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 t=1 is illustrated in Figure 16.Figure 16
Approximated solution using the finite element method.


## 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 t=1 is illustrated in Figure 16.Figure 16
Approximated solution using the finite element method.


## 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:(52a)∂u∂t−DΔu=ft,x,t,x∈0,T×Ω,(52b)ut,x=0,x∈∂Ω,(52c)u0,x=u0x≥0,x∈Ω,where 0,T is a given time interval and Ω=B0,1 is the unit disk with boundary ∂Ω. Here, x=x,yt denotes the space variables, t is the time variable, D is a negative diffusion coefficient, and f represents the effect of internal source terms, given by the following equality:(53)ft,x,y=C1ε2−2D−2x4−4x2y2+2ε4−2ε2y2−2y4−2ε2x2ε2−x2−y24Φxε,yεexpt.We can verify that the exact solution is given by(54)uext,x,y=Cε2Φxε,yεexpt,which 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.Figure 17
The exact solution att = 0.01, with D = −0.25, ε=0.01, maxuext,.=710.Figure 18
Theδ−ziti solution at t = 0.01, with D = −0.25, ε=0.01maxuapt,.=780, after 345-time iterations, using the Cartesian discretization.


## 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 ℝd. 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 PDEs. 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.



---

*Source: 6696165-2021-04-12.xml*