# Extension of δ_-ziti method in the unit ball: Numerical integration, resolution of Poisson's problem and Heat transfer

Inspired by the Galerkin and particular method, a new approximation approach is recalled in the Cartesian case. In this paper, we are interested specially by constructing this method, when the domain of consideration is a two dimensional ball, to extend this work to the several dimension. We reduce the number of iterations to calculate integrals and numerical solution of Poisson and the Heat problems (elliptic nd parabolic PDEs), in a very fast way.

There are no comments yet.

## Authors

• 1 publication
• 1 publication
03/19/2018

### Numerical Integration on Graphs: where to sample and how to weigh

Let G=(V,E,w) be a finite, connected graph with weighted edges. We are i...
04/19/2022

### Efficient Monte Carlo Method for Integral Fractional Laplacian in Multiple Dimensions

In this paper, we develop a Monte Carlo method for solving PDEs involvin...
05/11/2022

### Poisson Integrators based on splitting method for Poisson systems

We propose Poisson integrators for the numerical integration of separabl...
02/06/2020

### A C^1 Petrov-Galerkin method and Gauss collocation method for 1D general elliptic problems and superconvergence

In this paper, we present and study C^1 Petrov-Galerkin and Gauss colloc...
11/18/2019

### Energetic Stable Discretization for Non-Isothermal Electrokinetics Model

We propose an edge averaged finite element(EAFE) discretization to solve...
10/03/2021

### Numerical computation of Neumann controllers for the heat equation on a finite interval

This paper presents a new numerical method which approximates Neumann ty...
04/24/2012

### Geodesics in Heat

We introduce the heat method for computing the shortest geodesic distanc...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1. Introduction

The chemotactic dynamics of a population requires several steps, particularly, aggregation and blow-up. The Keller-Segel model describes this phenomena. It was suggested by Patlak in 1953 [3], Keller-Segel in 1970 [4], which allows for both diffusion and aggregation: depending on the initial data, the solution might exist globally in time or blow up in finite time, depending on the balance of forces between the different parameters involved in the system, the blow-up phenomenon may or may not occur. In fact, the blow-up is a singular behaviour of a Dirac solution. Most of numerical methods ( e.g.  Galerkin, Particular method, spectral method..) does not ensure the transition from a regular behaviour to another singular one (i.e.  detection of blow-up ). Under certain formulations of the Keller-Segel model, the phenomenon of aggregation has been shown to lead to finite-time blow-up. A large body of works has been devoted to determine when blow-up occurs or whether globally solutions exists: Authors of [5] developed a family of new interior penalty discontinuous Galerkin methods for solving the Keller–Segel chemotaxis model. In [6] they investigated non-negativity of exact and numerical solutions to a generalized Keller–Segel model where this model includes the so-called minimal Keller–Segel model. The main aim of [7] is to study the Keller–Segel model of chemotaxis and to develop a composite particle-grid numerical method with adaptive time stepping which allows to resolve and propagate singular solutions. The purpose of [8] is to formulate a phenomenological model from which the existence and properties of migrating bands can be deduced. Authors of [9] detected the blow-up as -function (amoebae aggregation) at the proximity of the origin in dimension one and on a ball in a multidimensional space. Therefore, it was necessary to find a new numerical scheme, which detect this type of singularities easily, without loosing the advantages of classical methods. The method is on the challenge. It was tested on several type of problems (see [1] and [2]), including the Keller-Segel model, but only in the Cartesian case ( segments, rectangle, cube ).The main goal of 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 famous function with compact support, defined by

 Φ(x)={exp(1∣x∣2−R2)si ∣x∣

where R 0, and .
This function is used especially in numerical analysis, distributions and functional analysis. It is characterized by giving the best approximation of the Dirac measure. In [1] and [2], the multi-dimensional Cartesian case was detailed. The main aim of this paper is the construction of method when the domain is a disk in the two-dimensional case, (in general, a multi-dimensional 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 to reconstruct our basis functions in every segment, which means that we inject all the work already done in the mono-dimensional case. To test this strategy, we apply the resulting tools to calculate numerically 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”). The second strategy is a direct use of the polar parametrisation of a disk, we will show that this strategy is also efficient and gives us a good approximation ( Integrals ans two tests of resolving PDEs).

The outline of this paper is as follows. In section 3, we will present the mathematical tools of construction, which permits to apply this method in the Cartesian case, as shown in figure 1, to calculate, numerically, some integrals defined in a disk domain.The 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 method.In the last section 5, we apply our approach to find the numerical solution of Poisson problem and the Heat equation. Our goal is to compare the solution obtained by method, with a given analytical one defined in a disk domain and to calculate the error in . By the next, we present an approximated solution using the finite element method and we compare it with our one.

## 2. Overview of the mono-dimensional construction.

All the results presents in this section, are proved in [1] and [2]. The fundamental results of construction are given as follow:First, we take a uniform mesh of [a,b] with the step , where N is an integer such that
From the function defined in (1), we define by:

 φϵ(x)=cϵΦ(xϵ),    for all  ϵ>0,

where c:= is the constant of normalization.

This sequence converges to Dirac in the sense of distributions, which is often used to detect singularities.

We construct the family as follows:

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩φi(x)=φh(x−xi)=ChΦ(x−xih),    for all  x∈[xi−1,xi+1],   i∈[2,N],φ1(x)=φh(x−x1)=ChΦ(x−x1h),    for all  x∈[x1,x2],φN+1(x)=φh(x−xN+1)=ChΦ(x−xN+1h),    for all  x∈[xN,xN+1].

Let 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

satisfying the following relation

 ⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩~Ψi(x)=φi(x)+λi−1~Ψi−1(x),λ1=−αβ,λi+1=g(λi), (2)

with

 g(X)=λ12−λ1X,  α=(φ1,φ2),  β=(φ1,φ1).

The spectral method applied to find the direct formula of the basis functions gives,

 λi=−(φi,ψi−1)(ψi,ψi−1)\raisebox2.15pt,

and the recurrence application of the definition given in (2) gives the following formula:

 ~Ψi(x)=φi(x)+λi−1φi−1(x)+λi−1λi−2φi−2(x)+⋯+1∏k=i−1λkφ1.

Let the normalization of ( for more details see [1] and [2]).
The method permits to approach a given function and an integral, using the following relations:

 f(x)≃N∑i=1ciΨi(x),ci≃∫baf(x)Ψi(x)dx,∫baf(x)dx≃N∑i=1ciIi, (3)

where . If we take in (3), we obtain:

 ci≃f(ri)Ψi(ri)\raisebox2.15pt,i=1⋯N−1,cN≃f(b)ΨN(b)\raisebox2.15pt,∫baf(x)Ψi(x)dx≃f(ri)Ψi(ri)\raisebox2.15pt,∫baf(x)dx≃N−1∑i=1f(ri)Ψ2i(ri)+f(b)ΨN(b)2\raisebox2.15pt. (4)

To reduce the iterations number, in [11] they proved that,

 |λi+1−λi|<ϵ  as soon as   i≥N0=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣ln(ϵ(2−λ21)λ31−λ1)ln(λ12+λ1)2⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦+1,

where denotes the floor function. In particular for   , we concluded that the parameter is nearly stationary from a certain rank, which reduces considerably the number of iterations. Using as a root of , we can define the parameter by,

## 3. The first strategy: Cartesian coordinates.

### 3.1. Construction of intern nodes

In this section, we are interested in the extension of method, when is a disk centred in the point (or the ball in the multi-dimensional case). As a first step, we start by a general presentation of this new strategy.We present the important steps of the construction, inspired by the mono-dimensional case. For this, suppose that we can 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. (see Figure 1)

An algorithm which explains the steps of the construction will be presented in the next part . The main idea is to fix the number of nodes in every interval and to vary the step of subdivision associated to every interval (horizontally and vertically). Every node is noted by , (see Figure 1).

###### Remark 3.1.

Note that, for every fixed vertical level (respectively the horizontal level ), every internal segment is limited by and , (respectively, and ).
For simplicity, the step of horizontal subdivision will be noted (respectively, the vertical step will be noted ).

Here we present an algorithm to calculate the internal nodes.

### 3.2. Construction of the orthonormal set

For every node (respectively, ) we associate the function , ( noted if there is no ambiguity) (respectively ,the family will be noted ) defined by:

 φi(x):=chjΦ(x−xjihj),  ∀i,j=1⋯N,φj(y):=chiΦ(y−yijhi),  ∀j,i=1⋯N, (5)

where,
is the step of construction in the horizontal interval of indication , which describe the distance between the node and .
is the step of subdivision in the vertical interval of indication , which describe 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,

 Horizontally:⎧⎪ ⎪⎨⎪ ⎪⎩~Ψ1(x)=φ1(x)~Ψi(x)=φi(x)+i−1∑k=1λ(i)k~Ψk(x),    for all  i=2,…,N,
 Vertically:⎧⎪ ⎪⎨⎪ ⎪⎩~Ψ1(y)=φ1(y)~Ψi(y)=φi(y)+i−1∑k=1λ(i)k~Ψk(y),    for all  i=2,…,N,

which will be reduced in the following theorem, already proved in the mono-dimensional case, (see [1] and [2]).

###### Theorem 3.1.

The orthogonal family (vertically and horizontally), verify the following recurrence relation:

 ⎧⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪⎩~Ψ1=φ1~Ψi+1=φi+1+λi~Ψi,    for all  i=1⋯N−1,λi−1=−(φi,φi−1)(~Ψi−1,~Ψi−1)\raisebox2.15pt, (6)

where is the usual scalar product in the Hilbert space

###### Corollary 3.2.

The family and the set defined in the theorem  (3.1), verify the following relations:

 1)   ~Ψi=φi+λi−1φi−1+λi−1λi−2φi−2+⋯+λi−1⋯λ1φ1.2)   ~Ψi(xji)=φi(xji)=ch2je.3)   ~Ψi(yji)=φi(yji)=ch2ie.4)   In every fixed level, (φi,~Ψi−1)=(φi,φi−1).5)   −1<λi=−(φi,φi+1)(~Ψi,~Ψi)<0. (7)

### 3.3. Fundamental results: Numerical integrations

In this paragraph, we are interested by 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:

 Ψij(x,y)=Ψji(x).Ψij(y),   ∀i,j=1⋯N,rij=(rji,sij), (8)

where,
are the basis functions in the horizontal dimension, (respectively, are the basis functions in the vertical dimension).
are the roots of (respectively, are the roots of ).
In this section, we are interested by the approximation of a double integral, defined in a disk domain, using (3), we obtain the following results:

###### Theorem 3.3.

Let , a given function in , and N denotes the roots number, therefore, we have the following approximations:

 ∫Ωg(x,y) dxdy≃N∑i,j=1g(rji,sij)Ψji(rji).Ψii(sij)∫bjajΨji(x)dx∫diciΨij(y)dy,∫Ωg(x,y)Ψij(x,y) dxdy≃g(rji,sij)(Ψji(rji).Ψii(sij))2\raisebox2.15pt. (9)

here we take, and .

In the table (1), we present some numerical tests of integration. We compare the exact value, with the numerical approximation, obtained by method, in the Cartesian case.

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

 ∀(x,y)∈B(0,1),x=rcos(θ),y=rsin(θ),(r,θ)∈[0,1]×[0,2π].

### 4.1. Construction of the method’s tools using polar coordinates

This first part of the algorithm, compute and in
The construction’s algorithm using the two variables and is given as follows::

The second part of the algorithm permits to compute and in

The polar set is defined by,

 Ψij(r,θ):=Ψi(r).Ψj(θ). (10)

### 4.2. Fundamental results: Numerical Integration

To test the previous strategy, we present in the following table, some numerical tests. We compare the exact value with the numerical one, using method in the polar case.

Using the polar coordinates and for some types of integrals, the following table shows us the error between exact and approximated solution founded using method.

In the previous table, we remark that even we choose a generalised integral, like the example , which is in fact an operation of Riemann integral, we found a good approximation using roots. For the two last examples, and , other approximation methods (e.g.  Simpson, Trapeze..) didn’t give any result, which is an important point for our construction.

## 5. Numerical applications.

### 5.1. Elliptic PDE case : Poisson problem

#### 5.1.1. The Cartesian case

In this section, let 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 by:

 −Δu  =f    in  Ω, (11a) u(x,y)=0    in  ∂Ω, (11b)

with a given analytical solution , and the source term function is defined by .

#### The strong discretization

The first step to approach the previous problem, is to multiply the equation (11a) by a test function and to integrate the result over the domain , which gives,

 −∫ΩΔu(x,y).Ψij(x,y)dxdy=∫Ωf(x,y).Ψij(x,y)dxdy. (12)

Using the theorem (3.3), we obtain the following scheme:

 −Δu(rji,sij)Ψij(rji,sij)=f(rji,sij)Ψij(rji,sij),∀(rji,sij)∈Ω. (13)

The next step, consists to approach the second derivative, which gives us the following scheme:

 ⎧⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪⎩ui−1,j+2uij+ui+1,j(rji+1−rji)(rji−rji−1)+ui,j−1+2uij+ui,j+1(sij+1−sij)(sij−sij−1)=f(rji,sij),  i,j=2⋯N−1,u1,j=uM,j=0,  j=1⋯N,ui,1=ui,M=0,  i=1⋯N, (14)

where, is the nods number in every internal segment (horizontally and vertically). At the end, we will have a global matrix, with lines and columns, defined as follows::

 M=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣D2A30⋯0A2D3A4⋯0⋮0⋯An−3Dn−2An−100⋯An−2Dn−1⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦

where is a tri-diagonal matrix, defined by:

 Di+1k,k=−2Ψi,k+1dxk.dxk−1−2Ψi,k+1dyk.dyk−1, i,k=1⋯N−2, Di+1k,k+1=−2Ψi,k+2dyk.dyk−1, i,k=1⋯N−2, Di+1k−1,k=−2Ψi,kdyk.dyk−1, i,k=1⋯N−2,

and is a diagonal matrix defined by:

 Aik,k=−Ψi,kdxk.dxk−1,  i=3⋯N−1,  and  k=2⋯N−1,

where,

 dxk=rjk−rjk−1\raisebox2.15pt, dyk=rik−rik−1\raisebox2.15pt,

therefore, we should resolve a simple system in the form , when is the global matrix defined previously,

in the unknown vector of size

and is the source term vector of size .

###### Remark 5.1.

To complete the resolution of the previous system, we must add boundary conditions ( homogeneous Dirichlet in this case).

#### 5.1.2. 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 by the shape of the approximated solution with 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 to the exact one, given by and .

In the table 4, we present the error between exact and approximated solution, with different values of nodes number .

where,

 Ermax :=maxi,j(|uex(i,j)−u(i,j))|, Ermean :=1N∑i,j|uex(i,j)−u(i,j))|.

We present in the following subsection, a comparison between approximated solutions using Finite Element and methods.

#### 5.1.3. Comparison with Finite elements method

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 [10] and [12]). The starting point for the FEM is a PDE expressed in 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 Galerkin analogy and . In this part, we present the approximated solution of the Poisson’s problem defined in (11a), using Finite Element Method.

#### 5.1.4. 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 method in the polar case.

Let , the strategy presented consists in using the results of approximation in the mono-dimensional case and taking into consideration the following function basis:

 Ψij(r,θ)=Ψi(r)Ψj(θ),  ∀i,j=1⋯N.

The problem presented in the previous subsection 5.1.1 is equivalent of the polar one, expressed as follows:

 −Δu  =f    in  Ω, (15a) u(r=1,θ)=0   ∀θ∈[0,2π], (15b)

with

 Δu:=∂2u∂r2+1r∂u∂r+1r2∂2u∂θ2, f(r,θ)=f(rcosθ,rsinθ)=4

where id 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 method, we should multiply the equation 15a by a test function and after, we use the strong result of approximation 3.3, which gives:

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩Δiju=f(ri,θj),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,