Abstract: A meshless method with ridge basis functions for solving the 2D two-flow domain model problem is proposed in this paper, which is using ridge basis functions to construct the approximate functions, and using collocation method to discretize the governing equation. The existence and uniqueness of the solution to the mothod is analyzed. Compared with the traditional finite difference method, the new method reduces the calculation error. Numerical results show that the proposed method is convergent and can improve the error precision.

Keywords: Two-flow domain model, Ridge basis functions, Meshless method, Collocation method.

# 1. Introduction

The migration of solutes in the soil is an important part of the biosphere cycle. By generalizing the flow velocity characteristics of soil pores, various mathematical models are proposed for description. For undisturbed soil, there are a certain number of lardge pores. The convective dispersion model[1] cannot explain the early penetration and tailing of the solute breakthrough curve. Therefore, Coats and Smith established a mobile and immobile model[2] to better describe the characteristics of solute migration. The mobile and immobile model exposed its limitations with the emergence of double peaks in the solute breakthrough curves and preferential flow in the soil. To this end, Skopp proposed a two-flow domain model based on the mobile and immobile model[3]. The model problem is often a partial differential equation system, and it is generally difficult to obtain the corresponding analytical solution. Thus, the numerical method is currently the most effective method to solve it.

Meshless method, which is based on point approximation and does not need initial division and reconstruction of grid, can completely or partially eliminate the grid. It is an efficient numerical method for solving partial differential equations in recent years[4]. With the development of meshless method, the current research based on radial basis functions has a series of results[5-8], but the research based on ridge basis functions is rare in theory and application. Gordon[9] considered the approximation of the norm space ridge basis functions; Ismailov[10] gave the best approximation of multivariate functions in L2 space by linear combination of ridge basis functions on some sets in Rn; After that, Zhang Liwei[11-12] gave the existence of the interpolation problem of the ridge basis functions in 2D space, estimated its error, and proved ridge basis functions approximation of the linear combination of finite plane waves; Shu Hengmu et al.[13] constructed an approximation function satisfying the interpolation theory of ridge basis functions, adopted the meshless method of the weighted least squares discrete equation and applied it to the elastic static problem; In recent years, Wang Zhigang et al.[14] discussed the collocation meshless method of ridge basis functions and gave the existence and uniqueness of the numerical solution. Both ridge basis functions and radial basis functions have the advantages of simple form and isotropy, but the application of the ridge basis functions is more extensive, such as the development equation, the solution of the dynamic system problem. In addition, the selection of the basis function and its parameters, the arrangement of the nodes have direct impact on the interpolation accuracy of radial basis functions and ridge basis functions. However, the direction selection has no effect on radial basis functions, but will affect the interpolation accuracy of ridge basis functions. This is one of the advantages of its application.

This paper mainly considers the 2D two-flow domain model problem, and its equation is expressed

 ${\displaystyle {\begin{cases}{\begin{array}{l}{\frac {\partial C_{A}}{\partial t}}=D_{A}\Delta C_{A}-V_{A}\left({\frac {\partial C_{A}}{\partial x}}+{\frac {\partial C_{A}}{\partial y}}\right)-{\frac {\alpha }{{\theta }_{A}}}\left(C_{A}-C_{B}\right)+f_{a},\left(x,y\right)\in {\Omega }_{2},t\in \left[0,T\right],\\{\frac {\partial C_{B}}{\partial t}}=D_{B}\Delta C_{B}-V_{B}\left({\frac {\partial C_{B}}{\partial x}}+{\frac {\partial C_{B}}{\partial y}}\right)-{\frac {\alpha }{{\theta }_{B}}}\left(C_{B}-C_{A}\right)+f_{b},\left(x,y\right)\in {\Omega }_{2},t\in \left[0,T\right],\\C_{A}\left(x,y,0\right)=C_{A0},C_{B}\left(x,y,0\right)=C_{B0},\left(x,y\right)\in {\Omega }_{2},\\C_{A}\left(x,y,t\right)=f_{1},C_{B}\left(x,y,t\right)=f_{2},\left(x,y\right)\in {\Omega }_{2},t\in \left[0,T\right].\end{array}}\end{cases}}}$
(1)

where i is the flow rate area, subscripts A, B are two flow regions, ${\displaystyle \Omega =\left[0,L\right]\times \left[0,L\right]}$ , ${\displaystyle D_{A}}$ , ${\displaystyle D_{B}}$ , ${\displaystyle V_{A}}$ , ${\displaystyle V_{B}}$ , ${\displaystyle {\theta }_{A}}$ , ${\displaystyle {\theta }_{B}}$ , ${\displaystyle \alpha }$ are constants greater than zero.

The remainder of this paper is organized as follows. Section 2 mainly introduces the meshless method with ridge basis functions. Firstly, definitions of the ridge basis function and its interpolation problem are given. Then, the implementation process of the collocation method is described. Finally, the discrete process of 2D two-flow domain model is deduced. Section 3 analyzes the existence and uniqueness of the solution to the proposed method. Section 4 gives numerical examples of 1D and 2D two-flow domain model to verify the effectiveness of the ridge basis function meshless method. Conclusions and future work are given at the end of the paper.

# 2. Meshless method with ridge basis functions

## 2.1 Ridge basis functions

Definition 2.1 If there exsit unary function ${\displaystyle \phi :R\rightarrow R}$ and vector ${\displaystyle a\in R^{d}\backslash \left\{0\right\}}$ such that ${\displaystyle f\left(x\right)=\phi \left(a\cdot x\right)}$ holds for any ${\displaystyle x\in R^{d}}$ , then ${\displaystyle f\left(x\right):R^{d}\rightarrow R}$ is a ridge basis function. If ridge basis functions ${\displaystyle f_{1},f_{2},\cdots ,f_{n}}$ are basises of the space, then the space is called the ridge basis function space.

The ridge basis function uses the unary function ${\displaystyle \phi }$ to act on the inner product ${\displaystyle a\cdot x}$ . It can be seen that the function space is a linear combination of some plane waves. Define the ridge basis function interpolation as follows:

Definition 2.2 Give ${\displaystyle n}$ different nodes ${\displaystyle {\left\{x_{i}\right\}}_{i=1}^{n}}$ and ${\displaystyle n}$ data points ${\displaystyle {\left\{d_{i}\right\}}_{i=1}^{n}}$ in ${\displaystyle R^{d}}$ , the function ${\displaystyle F(x)}$ statisfying

Definition 2.2 Give ${\displaystyle n}$ different nodes ${\displaystyle {\left\{x_{i}\right\}}_{i=1}^{n}}$ and ${\displaystyle n}$ data points ${\displaystyle {\left\{d_{i}\right\}}_{i=1}^{n}}$ in ${\displaystyle R^{d}}$ , the function ${\displaystyle F(x)}$ statisfying

 ${\displaystyle F\left(x_{j}\right)=d_{j},j=1,2,\cdots ,n,}$
(2)
 ${\displaystyle F\left(x\right)=\sum _{k=1}^{n}{\lambda }_{k}\sum _{v=1}^{m}\phi \left(a{\cdot }_{}^{v}\left(x-\right.\right.}$${\displaystyle \left.\left.x_{k}\right)\right),}$
(3)

where ${\displaystyle {\lambda }_{k}}$ is unknown coefficient, ${\displaystyle a_{v}}$ is fixed direction, ${\displaystyle m}$ is the total number of the disk directions, and ${\displaystyle \phi }$ is a ridge basis function, then the ridge basis function interpolation problem becomes the following system of linear equations:

 ${\displaystyle \sum _{k=1}^{n}{\lambda }_{k}\sum _{v=1}^{m}\phi \left(a_{v}\cdot \left(x_{j}-\right.\right.}$${\displaystyle \left.\left.x_{k}\right)\right)=F\left(x_{j}\right),j=}$${\displaystyle 1,2,\cdots ,n}$ .
(4)

This system has a unique solution if and only if its coefficient matrix is non-singular, which is related with the distribution of nodes and the selection of direction ${\displaystyle a_{v}}$ [14].

In the numerical calculation, the ridge basis functions commonly used are

Gaussians: ${\displaystyle \phi \left(r\right)=exp\left(-c{\left(a\cdot r\right)}^{2}\right)}$

Multi-Quadrics(MQ): ${\displaystyle \phi \left(r\right)={\left(c^{2}+{\left(a\cdot r\right)}^{2}\right)}^{\beta }}$

Markoff distribution functions: ${\displaystyle \phi \left(r\right)=exp\left(-c\vert a\cdot r\vert \right)}$

Thin-Plate Spline(TPS): ${\displaystyle \phi \left(r\right)={\left(a\cdot r\right)}^{2\beta }log\left(a\cdot r\right)}$

Inverse Multi-Quadrics(IMQ): ${\displaystyle \phi \left(r\right){\mbox{=}}{\left(c^{2}+{\left(a\cdot r\right)}^{2}\right)}^{-\beta }}$

where ${\displaystyle c}$ and ${\displaystyle \beta }$ are constants ( ${\displaystyle c>0}$ ), ${\displaystyle r=\Vert x-x_{k}\Vert }$ ..

## 2.2 Collocation method

The collocation method is a true meshless method, which does not need any grid to calculate the integral. Its basic idea is that each node in the domain satisfies the governing equation, and each node on the boundary satisfies the uniform boundary condition. Generally, the governing equation is set as

 ${\displaystyle L\left(u\left(x\right)\right)=f,x\in \Omega }$
(5)

Dirichlet boundary condition is

 ${\displaystyle u-{\overline {u}}=0,x\in {\Gamma }_{1}}$
(6)

Neumann boundary condition is

 ${\displaystyle B\left(u\right)=g,x\in {\Gamma }_{2}}$
(7)

where ${\displaystyle L}$ , ${\displaystyle B}$ are differential operators, ${\displaystyle \partial \Omega ={\Gamma }_{1}\cup {\Gamma }_{2}=\varnothing }$ is boundary, ${\displaystyle {\overline {u}}}$ is the specified value of the unknown function on the boundary ${\displaystyle {\Gamma }_{1}}$ , ${\displaystyle u\left(x\right)}$ is an unknown function that needs to be solved.

Let the approximate expression of the unknown function ${\displaystyle u\left(x\right)}$ be ${\displaystyle u^{h}\left(x\right)}$ , using the weighted residual value method, Eq.(5)~Eq.(7) can be replaced by the following formula

${\displaystyle {\int }_{\Omega }w_{1}\left(L\left(u^{h}\left(x\right)\right)-\right.}$${\displaystyle \left.f\right)d\Omega +{\int }_{{\Gamma }_{1}}w_{2}\left(u^{h}\left(x\right)-\right.}$${\displaystyle \left.{\overline {u}}\right)d\Gamma +{\int }_{{\Gamma }_{2}}w_{3}\left(B\left(u^{h}\left(x\right)\right)-\right.}$${\displaystyle \left.g\right)d\Gamma =0}$

where the weighting functions ${\displaystyle w_{1}}$ , ${\displaystyle w_{2}}$ , ${\displaystyle w_{3}}$ can have different definitions. Define ${\displaystyle {\delta }_{i}=w_{1}=w_{2}=w_{3}}$ , ${\displaystyle {\delta }_{i}}$ is Dirac ${\displaystyle \delta }$ function, thereby the following collocation type equation can be obtained

 ${\displaystyle L\left(u^{h}\left(x_{i}\right)\right)-f\left(x_{i}\right)=}$${\displaystyle 0,x_{i}\in \Omega ,i=1,2,\cdots ,n_{0}}$
(8)
 ${\displaystyle u^{h}\left(x_{i}\right)-{\overline {u}}=0,x_{i}\in {\Gamma }_{1},i=}$${\displaystyle n_{0}+n_{2}+1,\cdots ,n_{0}+n_{2}+n_{1}}$
(9)
 ${\displaystyle B\left(u^{h}\left(x_{i}\right)\right)-g\left(x_{i}\right)=}$${\displaystyle 0,x_{i}\in {\Gamma }_{2},i=n_{0}+1,\cdots ,n_{0}+n_{2}}$
(10)

where ${\displaystyle n_{1}}$ is the number of nodes on the boundary ${\displaystyle {\Gamma }_{1}}$ , ${\displaystyle n_{2}}$ is the number of nodes on the boundary ${\displaystyle {\Gamma }_{2}}$ , ${\displaystyle n_{0}=N-n_{1}-n_{2}}$ is the number of nodes in the domain ${\displaystyle \Omega }$ , ${\displaystyle N}$ is total number of nodes. Thus, Eq.(8)~Eq.(10) can be expressed in matrix form as

${\displaystyle Ku^{h}\left(x\right)=F}$

## 2.3 The construction of meshless method with ridge basis functions

Let ${\displaystyle C_{A}^{n}=C_{A}^{n}\left(X\right)}$ be the numerical solution of ${\displaystyle C_{A}\left(X,t_{n}\right)}$ , ${\displaystyle C_{B}^{n}=C_{B}^{n}\left(X\right)}$ be the numerical solution of ${\displaystyle C_{B}\left(X,t_{n}\right)}$ , where ${\displaystyle t_{n}=n\Delta t}$ , ${\displaystyle n=1,2,\cdots ,s}$ , ${\displaystyle \Delta t=T/s}$ , then Eq.(1) can be discretized in the time leyer as

${\displaystyle {\frac {C_{A}^{n+1}-C_{A}^{n}}{\Delta t}}=D_{A}({\frac {{\partial }^{2}C_{A}^{n+1}}{\partial x^{2}}}+}$${\displaystyle {\frac {{\partial }^{2}C_{A}^{n+1}}{\partial y^{2}}})-V_{A}({\frac {\partial C_{A}^{n+1}}{\partial x}}+}$${\displaystyle {\frac {\partial C_{A}^{n+1}}{\partial y}})-{\frac {\alpha }{{\theta }_{A}}}\left(C_{A}^{n+1}-\right.}$${\displaystyle \left.C_{B}^{n+1}\right)+f_{a}^{n+1},}$ ${\displaystyle {\frac {C_{B}^{n+1}-C_{B}^{n}}{\Delta t}}=D_{B}({\frac {{\partial }^{2}C_{B}^{n+1}}{\partial x^{2}}}+}$${\displaystyle {\frac {{\partial }^{2}C_{B}^{n+1}}{\partial y^{2}}})-V_{B}({\frac {\partial C_{B}^{n+1}}{\partial x}}+}$${\displaystyle {\frac {\partial C_{B}^{n+1}}{\partial y}})-{\frac {\alpha }{{\theta }_{B}}}\left(C_{B}^{n+1}-\right.}$${\displaystyle \left.C_{A}^{n+1}\right)+f_{b}^{n+1}.}$

If ${\displaystyle \phi \left(\cdot \right)}$ is not a polynomial, then ridge basis functions can approximate to almost all the functions. Here, we respectively denote functions ${\displaystyle C_{A}^{h}\left(X\right)}$ , ${\displaystyle C_{B}^{h}\left(X\right)}$ as approximations of ${\displaystyle C_{A}\left(X\right)}$ , ${\displaystyle C_{B}\left(X\right)}$ :

${\displaystyle C_{A}^{h}\left(X\right)=\sum _{i=1}^{N_{i}}{\alpha }_{1i}\sum _{v=1}^{m}\phi \left(a_{v}\cdot \left(X-\right.\right.}$${\displaystyle \left.\left.I_{i}\right)\right)+\sum _{i=1}^{N_{b}}{\beta }_{1i}\sum _{v=1}^{m}\phi \left(a_{v}\cdot \left(X-\right.\right.}$${\displaystyle \left.\left.B_{i}\right)\right),}$
${\displaystyle C_{B}^{h}\left(X\right)=\sum _{i=1}^{N_{i}}{\alpha }_{2i}\sum _{v=1}^{m}\phi \left(a_{v}\cdot \left(X-\right.\right.}$${\displaystyle \left.\left.I_{i}\right)\right)+\sum _{i=1}^{N_{b}}{\beta }_{2i}\sum _{v=1}^{m}\phi \left(a_{v}\cdot \left(X-\right.\right.}$${\displaystyle \left.\left.B_{i}\right)\right).}$

where ${\displaystyle X=\left(x,y\right)}$ , ${\displaystyle I_{i}\in \Omega }$ , ${\displaystyle B_{i}\in \partial \Omega }$ , ${\displaystyle N_{i}}$ is the number of nodes in the region, ${\displaystyle N_{b}}$ is the number of nodes on the boundary, ${\displaystyle {\alpha }_{1i}}$ , ${\displaystyle {\beta }_{1i}}$ , ${\displaystyle {\alpha }_{2i}}$ , ${\displaystyle {\beta }_{1i}}$ are unknown coefficients, ${\displaystyle a_{v}}$ is a fixed direction, ${\displaystyle m}$ is the total number of the disk directions, ${\displaystyle \phi }$ is a ridge basis function.

Let ${\displaystyle \varphi \left(\Vert X-I\Vert \right)=\sum _{v=1}^{m}\phi \left(a_{v}\cdot \left(X-\right.\right.}$${\displaystyle \left.\left.I\right)\right)}$ , then

 ${\displaystyle C_{A}^{h}\left(X\right)=\sum _{i=1}^{N_{i}}{\alpha }_{1i}\varphi \left(\Vert X-\right.}$${\displaystyle \left.I_{i}\Vert \right)+\sum _{i=1}^{N_{b}}{\beta }_{1i}\varphi \left(\Vert X-\right.}$${\displaystyle \left.B_{i}\Vert \right),}$
(11)
 ${\displaystyle C_{B}^{h}\left(X\right)=\sum _{i=1}^{N_{i}}{\alpha }_{2i}\varphi \left(\Vert X-\right.}$${\displaystyle \left.I_{i}\Vert \right)+\sum _{i=1}^{N_{b}}{\beta }_{2i}\varphi \left(\Vert X-\right.}$${\displaystyle \left.B_{i}\Vert \right).}$
(12)

By adopting ridge basis functions collocation method, it is required that Eq.(11) and Eq.(12) satisfy the partial differential equation in the solution domain, which can be obtained

${\displaystyle {\begin{array}{c}\left[\left(1+{\frac {\alpha \Delta t}{{\theta }_{A}}}\right)-\Delta tD_{A}\Delta +\Delta tV_{A}\nabla \right]\left[\sum _{i=1}^{N_{i}}{\alpha }_{1i}\varphi \left(\Vert X-I_{i}\Vert \right)+\sum _{i=1}^{N_{b}}{\beta }_{1i}\varphi \left(\Vert X-B_{i}\Vert \right)\right]\\=C_{A}^{n}+{\frac {\alpha \Delta t}{{\theta }_{A}}}C_{B}^{n+1}+\Delta tf_{a}^{n+1},\end{array}}}$
${\displaystyle {\begin{array}{c}\left[\left(1+{\frac {\alpha \Delta t}{{\theta }_{B}}}\right)-\Delta tD_{B}\Delta +\Delta tV_{B}\nabla \right]\left[\sum _{i=1}^{N_{i}}{\alpha }_{2i}\varphi \left(\Vert X-I_{i}\Vert \right)+\sum _{i=1}^{N_{b}}{\beta }_{2i}\varphi \left(\Vert X-B_{i}\Vert \right)\right]\\=C_{B}^{n}+{\frac {\alpha \Delta t}{{\theta }_{B}}}C_{A}^{n+1}+\Delta tf_{b}^{n+1}.\end{array}}}$

Thus there is

 ${\displaystyle {\begin{cases}{\begin{array}{l}\sum _{i=1}^{N_{i}}\left[\left(1+{\frac {\alpha \Delta t}{{\theta }_{A}}}\right)\varphi \left(\Vert I_{j}-I_{i}\Vert \right)-\Delta tD_{A}\Delta \varphi \left(\Vert I_{j}-I_{i}\Vert \right)+\Delta tV_{A}\nabla \varphi \left(\Vert I_{j}-I_{i}\Vert \right)\right]{\alpha }_{1}{}_{i}^{n+1}\\+\sum _{i=1}^{N_{b}}\left[\left(1+{\frac {\alpha \Delta t}{{\theta }_{A}}}\right)\varphi \left(\Vert I_{j}-B_{i}\Vert \right)-\Delta tD_{A}\Delta \varphi \left(\Vert I_{j}-B_{i}\Vert \right)+\Delta tV_{A}\nabla \varphi \left(\Vert I_{j}-B_{i}\Vert \right)\right]{\beta }_{1}{}_{i}^{n+1}\\{\mbox{ }}{\mbox{ }}{\mbox{ }}=C_{A}^{n}\left(I_{j}\right)+{\frac {\alpha \Delta t}{{\theta }_{A}}}C_{B}^{n+1}\left(I_{j}\right)+\Delta tf_{a}^{n+1}\left(I_{j}\right),j=1,2,\cdots ,N_{i}\\\sum _{i=1}^{N_{i}}\left[\left(1+{\frac {\alpha \Delta t}{{\theta }_{B}}}\right)\varphi \left(\Vert I_{j}-I_{i}\Vert \right)-\Delta tD_{B}\Delta \varphi \left(\Vert I_{j}-I_{i}\Vert \right)+\Delta tV_{B}\nabla \varphi \left(\Vert I_{j}-I_{i}\Vert \right)\right]{\alpha }_{2}{}_{i}^{n+1}\\+\sum _{i=1}^{N_{b}}\left[\left(1+{\frac {\alpha \Delta t}{{\theta }_{B}}}\right)\varphi \left(\Vert I_{j}-B_{i}\Vert \right)-\Delta tD_{B}\Delta \varphi \left(\Vert I_{j}-B_{i}\Vert \right)+\Delta tV_{B}\nabla \varphi \left(\Vert I_{j}-B_{i}\Vert \right)\right]{\beta }_{2}{}_{i}^{n+1}\\{\mbox{ }}{\mbox{ }}{\mbox{ }}=C_{B}^{n}\left(I_{j}\right)+{\frac {\alpha \Delta t}{{\theta }_{B}}}C_{A}^{n+1}\left(I_{j}\right)+\Delta tf_{b}^{n+1}\left(I_{j}\right),j=1,2,\cdots ,N_{i}\\\sum _{i=1}^{N_{i}}{\alpha }_{1}{}_{i}^{0}\varphi \left(\Vert I_{j}-I_{i}\Vert \right)+\sum _{i=1}^{N_{b}}{\beta }_{1}{}_{i}^{0}\varphi \left(\Vert I_{j}-B_{i}\Vert \right)=C_{A}{}_{j}^{0},j=1,2,\cdots ,N_{i}+N_{b}\\\sum _{i=1}^{N_{i}}{\alpha }_{2}{}_{i}^{0}\varphi \left(\Vert I_{j}-I_{i}\Vert \right)+\sum _{i=1}^{N_{b}}{\beta }_{2}{}_{i}^{0}\varphi \left(\Vert I_{j}-B_{i}\Vert \right)=C_{B}{}_{j}^{0},j=1,2,\cdots ,N_{i}+N_{b}\\\sum _{i=1}^{N_{i}}{\alpha }_{1}{}_{i}^{n+1}\varphi \left(\Vert B_{j}-I_{i}\Vert \right)+\sum _{i=1}^{N_{b}}{\beta }_{1}{}_{i}^{n+1}\varphi \left(\Vert B_{j}-B_{i}\Vert \right)=f_{1}\left(B_{j}\right),j=1,2,\cdots ,N_{b}\\\sum _{i=1}^{N_{i}}{\alpha }_{2}{}_{i}^{n+1}\varphi \left(\Vert B_{j}-I_{i}\Vert \right)+\sum _{i=1}^{N_{b}}{\beta }_{2}{}_{i}^{n+1}\varphi \left(\Vert B_{j}-B_{i}\Vert \right)=f_{2}\left(B_{j}\right),j=1,2,\cdots ,N_{b}\end{array}}\end{cases}}}$
(13)

and let

${\displaystyle {\psi }_{A}\left(\Vert X-I\Vert \right)=\left(1+{\frac {\alpha \Delta t}{{\theta }_{A}}}\right)\varphi \left(\Vert X-\right.}$${\displaystyle \left.I\Vert \right)-\Delta tD_{A}\Delta \varphi \left(\Vert X-\right.}$${\displaystyle \left.I\Vert \right)+\Delta tV_{A}\nabla \varphi \left(\Vert X-\right.}$${\displaystyle \left.I\Vert \right),}$
${\displaystyle {\psi }_{B}\left(\Vert X-I\Vert \right)=\left(1+{\frac {\alpha \Delta t}{{\theta }_{B}}}\right)\varphi \left(\Vert X-\right.}$${\displaystyle \left.I\Vert \right)-\Delta tD_{B}\Delta \varphi \left(\Vert X-\right.}$${\displaystyle \left.I\Vert \right)+\Delta tV_{B}\nabla \varphi \left(\Vert X-\right.}$${\displaystyle \left.I\Vert \right),}$
${\displaystyle F_{A}^{n+1}\left(I_{j}\right)=C_{A}^{n}\left(I_{j}\right)+{\frac {\alpha \Delta t}{{\theta }_{A}}}C_{B}^{n+1}\left(I_{j}\right)+}$${\displaystyle \Delta tf_{a}^{n+1}\left(I_{j}\right),}$
${\displaystyle F_{B}^{n+1}\left(I_{j}\right)=C_{B}^{n}\left(I_{j}\right)+{\frac {\alpha \Delta t}{{\theta }_{B}}}C_{A}^{n+1}\left(I_{j}\right)+}$${\displaystyle \Delta tf_{b}^{n+1}\left(I_{j}\right).}$

Then Eq.(13) can be expressed in the following matrix equcation

 ${\displaystyle \lbrace {\begin{array}{c}H_{1}a=F_{1},\\H_{2}b=F_{2}.\end{array}}}$
(14)

where ${\displaystyle H_{1}}$ , ${\displaystyle H_{2}}$ are ${\displaystyle \left(N_{b}+N_{i}\right)\times \left(N_{b}+N_{i}\right)}$ matrices, ${\displaystyle a}$ , ${\displaystyle b}$ are the undetermined coefficient vectors, ${\displaystyle F_{1}}$ , ${\displaystyle F_{2}}$ are the given right end term vectors, namely

${\displaystyle H_{1}=\left[{\begin{array}{llllll}\varphi \left(\Vert B_{1}-B_{1}\Vert \right)&\cdots &\varphi \left(\Vert B_{1}-B_{N_{b}}\Vert \right)&\varphi \left(\Vert B_{1}-I_{1}\Vert \right)&\cdots &\varphi \left(\Vert B_{1}-I_{N_{i}}\Vert \right)\\{\mbox{ }}{\mbox{ }}\vdots &&{\mbox{ }}{\mbox{ }}\vdots &{\mbox{ }}{\mbox{ }}\vdots &&{\mbox{ }}{\mbox{ }}\vdots \\\varphi \left(\Vert B_{N_{b}}-B_{1}\Vert \right)&\cdots &\varphi \left(\Vert B_{N_{b}}-B_{N_{b}}\Vert \right)&\varphi \left(\Vert B_{N_{b}}-I_{1}\Vert \right)&\cdots &\varphi \left(\Vert B_{N_{b}}-I_{N_{i}}\Vert \right)\\{\psi }_{A}\left(\Vert I_{1}-B_{1}\Vert \right)&\cdots &{\psi }_{A}\left(\Vert I_{1}-B_{N_{b}}\Vert \right)&{\psi }_{A}\left(\Vert I_{1}-I_{1}\Vert \right)&\cdots &{\psi }_{A}\left(\Vert I_{1}-I_{N_{i}}\Vert \right)\\{\mbox{ }}{\mbox{ }}\vdots &&{\mbox{ }}{\mbox{ }}\vdots &{\mbox{ }}{\mbox{ }}\vdots &&{\mbox{ }}{\mbox{ }}\vdots \\{\psi }_{A}\left(\Vert I_{N_{i}}-B_{1}\Vert \right)&\cdots &{\psi }_{A}\left(\Vert I_{N_{i}}-B_{N_{b}}\Vert \right)&{\psi }_{A}\left(\Vert I_{N_{i}}-I_{1}\Vert \right)&\cdots &{\psi }_{A}\left(\Vert I_{N_{i}}-I_{N_{i}}\Vert \right)\end{array}}\right],}$ ${\displaystyle a={\left[{\beta }_{11}^{n+1},\cdots ,{\beta }_{1N_{b}}^{n+1},{\alpha }_{11}^{n+1},\cdots ,{\alpha }_{1N_{i}}^{n+1}\right]}^{T},}$
${\displaystyle F_{1}={\left[f_{1}\left(B_{1}\right),\cdots ,f_{1}\left(B_{N_{b}}\right),F_{A}^{n+1}\left(I_{1}\right),\cdots ,F_{A}^{n+1}\left(I_{N_{i}}\right)\right]}^{T}.}$

${\displaystyle H_{2}=\left[{\begin{array}{llllll}\varphi \left(\Vert B_{1}-B_{1}\Vert \right)&\cdots &\varphi \left(\Vert B_{1}-B_{N_{b}}\Vert \right)&\varphi \left(\Vert B_{1}-I_{1}\Vert \right)&\cdots &\varphi \left(\Vert B_{1}-I_{N_{i}}\Vert \right)\\{\mbox{ }}{\mbox{ }}\vdots &&{\mbox{ }}{\mbox{ }}\vdots &{\mbox{ }}{\mbox{ }}\vdots &&{\mbox{ }}{\mbox{ }}\vdots \\\varphi \left(\Vert B_{N_{b}}-B_{1}\Vert \right)&\cdots &\varphi \left(\Vert B_{N_{b}}-B_{N_{b}}\Vert \right)&\varphi \left(\Vert B_{N_{b}}-I_{1}\Vert \right)&\cdots &\varphi \left(\Vert B_{N_{b}}-I_{N_{i}}\Vert \right)\\{\psi }_{B}\left(\Vert I_{1}-B_{1}\Vert \right)&\cdots &{\psi }_{B}\left(\Vert I_{1}-B_{N_{b}}\Vert \right)&{\psi }_{B}\left(\Vert I_{1}-I_{1}\Vert \right)&\cdots &{\psi }_{B}\left(\Vert I_{1}-I_{N_{i}}\Vert \right)\\{\mbox{ }}{\mbox{ }}\vdots &&{\mbox{ }}{\mbox{ }}\vdots &{\mbox{ }}{\mbox{ }}\vdots &&{\mbox{ }}{\mbox{ }}\vdots \\{\psi }_{B}\left(\Vert I_{N_{i}}-B_{1}\Vert \right)&\cdots &{\psi }_{B}\left(\Vert I_{N_{i}}-B_{N_{b}}\Vert \right)&{\psi }_{B}\left(\Vert I_{N_{i}}-I_{1}\Vert \right)&\cdots &{\psi }_{B}\left(\Vert I_{N_{i}}-I_{N_{i}}\Vert \right)\end{array}}\right],}$

${\displaystyle b={\left[{\beta }_{21}^{n+1},\cdots ,{\beta }_{2N_{b}}^{n+1},{\alpha }_{21}^{n+1},\cdots ,{\alpha }_{2N_{i}}^{n+1}\right]}^{T},}$
${\displaystyle F_{2}={\left[f_{2}\left(B_{1}\right),\cdots ,f_{2}\left(B_{N_{b}}\right),F_{B}^{n+1}\left(I_{1}\right),\cdots ,F_{B}^{n+1}\left(I_{N_{i}}\right)\right]}^{T}.}$

# 3. Existence and uniqueness of solution

In the two-flow domain model, the region with faster flow rate is defined as region A, and the region with slower flow rate is defined as region B. For region A, let ${\displaystyle L=\left(1+\alpha \Delta t/{\theta }_{A}\right)-\Delta tD_{A}\Delta +}$${\displaystyle \Delta tV_{A}\nabla }$ , in order to discuss the existence and uniqueness of the solution of (14), we first perform coordinate transformation on L. We take ${\displaystyle Y^{n+1}=QX^{n+1}}$ ，where ${\displaystyle Q}$ is a non-degenerate matrix，and record as

${\displaystyle v_{A}^{n+1}\left(Y^{n+1}\right)=C_{A}^{n+1}\left(P^{-1}Y^{n+1}\right)}$

then there is

 ${\displaystyle L=\left(1+{\frac {\alpha \Delta t}{{\theta }_{A}}}\right)-}$${\displaystyle \Delta tD_{A}{\epsilon }_{k}\Delta +\Delta tV_{A}\nabla }$
(15)

Let ${\displaystyle \lambda ={\left({\lambda }_{1},{\lambda }_{2}\right)}^{T}}$ , ${\displaystyle v_{A}^{n+1}\left(Y^{n+1}\right)=w_{A}^{n+1}\left(Y^{n+1}\right)e^{\lambda \cdot Y^{n+1}}}$ , we have

 ${\displaystyle Dv_{A}^{n+1}\left(Y^{n+1}\right)=\lambda w_{A}^{n+1}\left(Y^{n+1}\right)e^{\lambda \cdot Y^{n+1}}+}$${\displaystyle e^{\lambda \cdot Y^{n+1}}Dw_{A}^{n+1}\left(Y^{n+1}\right),}$
(16)
 ${\displaystyle D^{2}v_{A}^{n+1}\left(Y^{n+1}\right)=w_{A}^{n+1}\left(Y^{n+1}\right)e^{\lambda \cdot Y^{n+1}}{\lambda }^{T}\lambda +}$${\displaystyle e^{\lambda \cdot Y^{n+1}}\left({\lambda }^{T}Dw_{A}^{n+1}\left(Y^{n+1}\right)+\right.}$${\displaystyle \left.\lambda Dw_{A}^{n+1}\left(Y^{n+1}\right)\right)+}$${\displaystyle e^{\lambda \cdot Y^{n+1}}D^{2}w_{A}^{n+1}\left(Y^{n+1}\right).}$
(17)

By substituting Eq.(16), (17) into Eq.(15), we can get

 ${\displaystyle -\Delta tD_{A}{\epsilon }_{k}\Delta w_{A}^{n+1}\left(Y^{n+1}\right)-}$${\displaystyle \left(2\Delta tD_{A}{\epsilon }_{k}\lambda -\Delta tV_{A}\right)\nabla w_{A}^{n+1}\left(Y^{n+1}\right)+}$${\displaystyle \left(\Delta tV_{A}\lambda -\Delta tD_{A}{\epsilon }_{k}{\lambda }^{T}\lambda +\right.}$${\displaystyle \left.\left(1+{\frac {\alpha \Delta t}{{\theta }_{A}}}\right)\right)w_{A}^{n+1}\left(Y^{n+1}\right)=}$${\displaystyle L\cdot e^{-\lambda Y^{n+1}}.}$
(18)

Substitute Eq.(18) into Eq.(15), and then combine with Eq.(13) to obtain

${\displaystyle \sum _{i=1}^{N_{i}}\varphi \left(\Vert {\overline {B}}_{j}-\right.}$${\displaystyle \left.{\overline {I}}_{i}\Vert \right){\alpha }_{1i}^{n+1}+}$${\displaystyle \sum _{i=1}^{N_{b}}\varphi \left(\Vert {\overline {B}}_{j}-\right.}$${\displaystyle \left.{\overline {B}}_{i}\Vert \right){\beta }_{1i}^{n+1}=}$${\displaystyle {\overline {f}}_{1}\left({\overline {B}}_{j}\right),}$
${\displaystyle \sum _{i=1}^{N_{i}}\left({\overline {c}}-\Delta tD_{A}{\epsilon }_{k}\Delta \right)\varphi \left(\Vert {\overline {I}}_{j}-\right.}$${\displaystyle \left.{\overline {I}}_{i}\Vert \right){\alpha }_{1i}^{n+1}+}$${\displaystyle \sum _{i=1}^{N_{b}}\left({\overline {c}}-\Delta tD_{A}{\epsilon }_{k}\Delta \right)\varphi \left(\Vert {\overline {I}}_{j}-\right.}$${\displaystyle \left.{\overline {B}}_{i}\Vert \right){\beta }_{1i}^{n+1}=}$${\displaystyle {\overline {F}}_{A}^{n+1}\left({\overline {I}}_{j}\right).}$

where let ${\displaystyle {\overline {\psi }}_{A}\left(\Vert {\overline {X}}-{\overline {I}}\Vert \right)=}$${\displaystyle \left({\overline {c}}-\Delta tD_{A}{\epsilon }_{k}\Delta \right)\varphi \left(\Vert {\overline {X}}-\right.}$${\displaystyle \left.{\overline {I}}\Vert \right)}$ ; ${\displaystyle T={\overline {c}}-\Delta tD_{A}{\epsilon }_{k}\Delta \left(k=\right.}$${\displaystyle \left.1,2\right)}$ , we take ${\displaystyle {\epsilon }_{k}=1}$ , ${\displaystyle {\overline {c}}}$ is a constant ${\displaystyle \left({\overline {c}}>0\right)}$ ; ${\displaystyle {\overline {F}}_{A}^{n+1}=F_{A}^{n+1}\left(P^{-1}Y^{n+1}\right)e^{-\lambda Y^{n+1}}}$ , ${\displaystyle \lambda ={\left({\lambda }_{1},{\lambda }_{2}\right)}^{T}}$ , ${\displaystyle {\lambda }_{k}=V_{A}/2D_{A}{\epsilon }_{k}}$ , ${\displaystyle {\overline {f}}_{1}=f_{1}\left(P^{-1}Y^{n+1}\right)}$ ; ${\displaystyle {\overline {B}}}$ , ${\displaystyle {\overline {I}}}$ are the points after ${\displaystyle B}$ , ${\displaystyle I}$ coordinate transformation.

Therefore, defined

${\displaystyle {\Psi }_{1}{\mbox{=}}\left[{\begin{array}{ccc}\varphi \left(\Vert {\overline {B}}_{1}-{\overline {B}}_{1}\Vert \right)&\cdots &\varphi \left(\Vert {\overline {B}}_{1}-{\overline {B}}_{N_{b}}\Vert \right)\\\vdots &\ddots &\vdots \\\varphi \left(\Vert {\overline {B}}_{N_{b}}-{\overline {B}}_{1}\Vert \right)&\cdots &\varphi \left(\Vert {\overline {B}}_{N_{b}}-{\overline {B}}_{N_{b}}\Vert \right)\end{array}}\right],}$ ${\displaystyle {\Psi }_{2}{\mbox{=}}\left[{\begin{array}{ccc}\varphi \left(\Vert {\overline {B}}_{1}-{\overline {I}}_{1}\Vert \right)&\cdots &\varphi \left(\Vert {\overline {B}}_{1}-{\overline {I}}_{N_{i}}\Vert \right)\\\vdots &\ddots &\vdots \\\varphi \left(\Vert {\overline {B}}_{N_{b}}-{\overline {I}}_{1}\Vert \right)&\cdots &\varphi \left(\Vert {\overline {B}}_{N_{b}}-{\overline {I}}_{N_{i}}\Vert \right)\end{array}}\right],}$
${\displaystyle {\Psi }_{3}{\mbox{=}}\left[{\begin{array}{ccc}{\overline {\psi }}_{A}\left(\Vert {\overline {I}}_{1}-{\overline {B}}_{1}\Vert \right)&\cdots &{\overline {\psi }}_{A}\left(\Vert {\overline {I}}_{1}-{\overline {B}}_{N_{b}}\Vert \right)\\\vdots &\ddots &\vdots \\{\overline {\psi }}_{A}\left(\Vert {\overline {I}}_{N_{i}}-{\overline {B}}_{1}\Vert \right)&\cdots &{\overline {\psi }}_{A}\left(\Vert {\overline {I}}_{N_{i}}-{\overline {B}}_{N_{b}}\Vert \right)\end{array}}\right],}$ ${\displaystyle {\Psi }_{4}{\mbox{=}}\left[{\begin{array}{ccc}{\overline {\psi }}_{A}\left(\Vert {\overline {I}}_{1}-{\overline {I}}_{1}\Vert \right)&\cdots &{\overline {\psi }}_{A}\left(\Vert {\overline {I}}_{1}-{\overline {I}}_{N_{i}}\Vert \right)\\\vdots &\ddots &\vdots \\{\overline {\psi }}_{A}\left(\Vert {\overline {I}}_{N_{i}}-{\overline {I}}_{1}\Vert \right)&\cdots &{\overline {\psi }}_{A}\left(\Vert {\overline {I}}_{N_{i}}-{\overline {I}}_{N_{i}}\Vert \right)\end{array}}\right].}$

then there is

${\displaystyle {\Psi }_{A}=\left[{\begin{array}{cc}{\Psi }_{1}&{\Psi }_{2}\\{\Psi }_{3}&{\Psi }_{4}\end{array}}\right],{\overline {F}}_{1}={\left[{\overline {f}}_{1}\left({\overline {B}}_{1}\right),\cdots ,{\overline {f}}_{1}\left({\overline {B}}_{N_{b}}\right),{\overline {F}}_{A}^{n+1}\left({\overline {I}}_{1}\right),\cdots ,{\overline {F}}_{A}^{n+1}\left({\overline {I}}_{N_{i}}\right)\right]}^{T}.}$

Consequently, the matrix equation is

 ${\displaystyle {\Psi }_{A}a={\overline {F}}_{1}}$
(19)

Lemma 1[15] The function ${\displaystyle \Phi \left(x\right)=\phi \left(a{\dot {x}}\right)}$ is a positive definite function if and only if its Fourier tarnsform ${\displaystyle F\left(\Phi \right)}$ is almost everywhere greater than 0.

Lemma 2[11] The Fourier transform ${\displaystyle F\left[\phi \right]\left(w\right)}$ of the ridge basis function ${\displaystyle \phi \left(r\right)}$ is almost everywhere greater than 0, the interpolation matrix ${\displaystyle Q}$ is definite.

Theorem 1 If the ridge basis function ${\displaystyle \phi \left(x\right)}$ is a positive definite function and ${\displaystyle {\left({\Psi }_{3}{\Psi }_{1}^{-1}{\Psi }_{2}-{\Psi }_{4}\right)}^{-1}}$ exists, Eq.(19) has a unique solution.

Proof. Choose the ridge basis function ${\displaystyle \phi \left(x\right)}$ , which is an even function. It is known from Lemma 1 and Lemma 2 that ${\displaystyle {\Psi }_{1}}$ is definite. Hence, ${\displaystyle {\Psi }_{1}^{-1}}$ exists. As ${\displaystyle -{\phi }^{''}\left(x\right)}$ is also an even function, then its Fourier transform is

${\displaystyle F\left[-{\phi }^{''}\left(x\right)\right]=-{\left(iw\right)}^{2}F\left[\phi \left(x\right)\right]=}$${\displaystyle w^{2}F\left[\phi \left(x\right)\right]>0.}$

Similarly, ${\displaystyle {\Psi }_{4}}$ is also definite. Hence, ${\displaystyle {\Psi }_{4}^{-1}}$ also exists.

Elementary transformation of ${\displaystyle \left[{\begin{array}{cc}{\Psi }_{1}&{\Psi }_{2}\\{\Psi }_{3}&{\Psi }_{4}\end{array}}\right]}$ , there is

${\displaystyle \left[{\begin{array}{cc}E&-{\Psi }_{2}{\Psi }_{4}^{-1}\\0&E\end{array}}\right]\left[{\begin{array}{cc}{\Psi }_{1}&{\Psi }_{2}\\{\Psi }_{3}&{\Psi }_{4}\end{array}}\right]{\mbox{=}}\left[{\begin{array}{cc}{\Psi }_{1}-{\Psi }_{2}{\Psi }_{4}^{-1}{\Psi }_{3}&0\\{\Psi }_{3}&{\Psi }_{4}\end{array}}\right]}$

We can get the following determinant

${\displaystyle \vert {\begin{array}{cc}E&-{\Psi }_{2}{\Psi }_{4}^{-1}\\0&E\end{array}}\vert \vert {\begin{array}{cc}{\Psi }_{1}&{\Psi }_{2}\\{\Psi }_{3}&{\Psi }_{4}\end{array}}\vert =\vert {\begin{array}{cc}{\Psi }_{1}-{\Psi }_{2}{\Psi }_{4}^{-1}{\Psi }_{3}&0\\{\Psi }_{3}&{\Psi }_{4}\end{array}}\vert =\vert {\Psi }_{1}-{\Psi }_{2}{\Psi }_{4}^{-1}{\Psi }_{3}\vert \cdot \vert {\Psi }_{4}\vert }$

and because

${\displaystyle \left({\Psi }_{1}-{\Psi }_{2}{\Psi }_{4}^{-1}{\Psi }_{3}\right)\left[{\Psi }_{1}^{-1}-\right.}$${\displaystyle \left.{\Psi }_{1}^{-1}{\Psi }_{2}{\left({\Psi }_{3}{\Psi }_{1}^{-1}{\Psi }_{2}-{\Psi }_{4}\right)}^{-1}{\Psi }_{3}{\Psi }_{1}^{-1}\right]=}$${\displaystyle E}$

we can get that ${\displaystyle {\left({\Psi }_{3}{\Psi }_{1}^{-1}{\Psi }_{2}-{\Psi }_{4}\right)}^{-1}}$ exists. Subsequently, ${\displaystyle \vert {\Psi }_{A}\vert \not =0}$ , the equation has a unique solution.

Similarly, for the area B, it can be proved that there is a unique solution to its matrix equation.

In summary, there is a unique solution to Eq.(14).

4. Numerical examples

In this section, numerical examples of 1D and 2D two-flow domain model are given by using the method described above. In order to verify the validity and accuracy of the proposed method, we compare it with the finite difference method (FDM). The relative error formula given here is

${\displaystyle L_{2}={\sqrt {\sum _{i=1}^{N}{\left[u\left(x_{i}\right)-u_{exact}\left(x_{i}\right)\right]}^{2}/\sum _{i=1}^{N}{\left[u\left(x_{i}\right)\right]}^{2}}}.}$

Example 1 Consider the initial boundary value problem of 1D two-flow domain model

${\displaystyle {\begin{cases}{\begin{array}{l}{\frac {\partial C_{A}}{\partial t}}=D_{A}{\frac {{\partial }^{2}C_{A}}{\partial x^{2}}}-V_{A}{\frac {\partial C_{A}}{\partial x}}-{\frac {\alpha }{{\theta }_{A}}}\left(C_{A}-C_{B}\right)+f_{a},x\in {\Omega }_{1},t\in \left[0,T\right],\\{\frac {\partial C_{B}}{\partial t}}=D_{B}{\frac {{\partial }^{2}C_{B}}{\partial x^{2}}}-V_{B}{\frac {\partial C_{B}}{\partial x}}-{\frac {\alpha }{{\theta }_{B}}}\left(C_{B}-C_{A}\right)+f_{b},x\in {\Omega }_{1},t\in \left[0,T\right],\\C_{A}\left(x,0\right)=C_{A0},C_{B}\left(x,0\right)=C_{B0},x\in {\Omega }_{1},\\C_{A}\left(0,t\right)=f_{a1},C_{A}\left(1,t\right)=f_{a2},C_{B}\left(0,t\right)=f_{b1},C_{B}\left(1,t\right)=f_{b2},t\in \left[0,T\right].\end{array}}\end{cases}}}$

Let ${\displaystyle D_{A}=D_{B}=1}$ , ${\displaystyle V_{A}=V_{B}=1}$ , ${\displaystyle {\theta }_{A}={\theta }_{B}=2}$ , ${\displaystyle \alpha =0.001}$ , ${\displaystyle {\Omega }_{1}=\left[0,1\right]}$ . The exact solutions are ${\displaystyle C_{A}\left(x,t\right)=t^{2}x{\left(x-1\right)}^{2}}$ , ${\displaystyle C_{B}\left(x,t\right)=t^{2}x^{2}{\left(x-1\right)}^{3}}$ , where ${\displaystyle C_{A0}}$ , ${\displaystyle C_{B0}}$ , ${\displaystyle f_{a1}}$ , ${\displaystyle f_{a2}}$ , ${\displaystyle f_{b1}}$ , ${\displaystyle f_{b2}}$ are determined by exact solutions. The space step is ${\displaystyle h=1/(N-1)}$ ( ${\displaystyle N}$ is the number of nodes) , the time step is ${\displaystyle \Delta t=0.001}$ . The Gaussian function ${\displaystyle \phi \left(r\right)=exp\left(-c{\left(a\cdot r\right)}^{2}\right)}$ as basis function, where the direction of calculation is selected from a unit circle. The direction is:

${\displaystyle a_{v}=cos\left({\theta }_{v}\right)}$${\displaystyle {\theta }_{v}=\left(v-1\right)\pi /\left(m-1\right)}$${\displaystyle v=1,2,\ldots ,m}$

When ${\displaystyle m=3}$ , calculation results are shown in Table 1, where ${\displaystyle c_{1}}$ , ${\displaystyle c_{2}}$ are parameters of ${\displaystyle C_{A}}$ , ${\displaystyle C_{B}}$ , ${\displaystyle e_{1}}$ , ${\displaystyle e_{2}}$ are calculation errors of ${\displaystyle C_{A}}$ , ${\displaystyle C_{B}}$ , respectively.

Table 1 The results of the numerical solution of the method in this paper (T=1)
 Inserted nodes (N) Space step (h) Parameter (c) Calculation error (e) 11 1/10 c1=4.55 e1=7.2261e-04 c2=2.51 e2=8.5523e-04 21 1/20 c1=28.0 e1=6.2200e-04 c2=25.0 e2=6.1254e-04 51 1/50 c1=217 e1=4.7024e-04 c2=247 e2=3.8390e-04

Figure 1 Comparison of exact solutions and numerical solutions (N=31)

According to Table 1, it can be seen that the meshless method with ridge basis functions is feasible, and calculation accuracy has gradually improved with the increase of nodes. Fig.1 shows the comparison between the exact solution and the numerical solution of the method in this paper when the number of nodes is 31. It can be seen from Fig.1 that the numerical solution of the proposed method is well matched with the exact solution even when T is large.

Table 2 Comparison of errors between the method in this paper and FDM
 Time Time space Inserted nodes The method of this paper FDM T=1 0.001 31 5.5257e-04 3.1313e-03 4.8538e-04 3.4988e-03 61 4.0158e-04 2.5859e-03 3.4339e-04 2.3371e-03 101 3.6495e-04 2.4729e-03 2.3089e-04 1.9831e-03 0.01 31 4.2783e-03 2.4518e-02 6.8133e-03 2.9222e-02 61 2.7401e-03 2.4005e-02 2.8381e-03 2.2630e-02 101 1.2704e-03 2.3896e-02 1.3564e-03 1.3382e-02

Table3 Comparison of errors of the method in this paper
 Time Inserted nodes Parameter m=4 Calculation error m=5 Calculation error m=6 Calculation error T=1 11 c1=4.550 e1=6.2837e-04 e1=6.1823e-04 e1=5.6218e-04 c2=2.450 e2=1.0022e-03 e2=3.7018e-03 e2=1.2782e-03 41 c1=127.0 e1=5.7576e-04 e1=6.1657e-04 e1=8.4694e-04 c2=134.0 e2=2.6641e-04 e2=1.9289e-04 e2=4.0056e-04 101 c1=905.7 e1=1.5110e-04 e1=9.0536e-04 e1=1.4458e-04 c2=903.7 e2=2.2702e-04 e2=1.5241e-04 e2=2.2188e-04

When T=1, Table 2 gives the comparison of the calculation results of the proposed method and the FDM under different time steps and different node numbers. Table 3 shows the error comparison of the numerical solution of the proposed method after selecting different nodes and directions. As can be seen from Table 2, under different time steps, as the number of spatial nodes increases, the error accuracy and the convergent result of the proposed method are over than the FDM. It can be seen from Table 3 that the number of inserted nodes, the selection of free parameters and the direction of the basis function all affect the numerical accuracy of the meshless method with ridge basis functions.

Example 2 Consider the initial boundary value problem of 2D two-flow domain model

${\displaystyle {\begin{cases}{\begin{array}{l}{\frac {\partial C_{A}}{\partial t}}=D_{A}\Delta C_{A}-V_{A}\left({\frac {\partial C_{A}}{\partial x}}+{\frac {\partial C_{A}}{\partial y}}\right)-{\frac {\alpha }{{\theta }_{A}}}\left(C_{A}-C_{B}\right)+f_{a},\left(x,y\right)\in {\Omega }_{2},t\in \left[0,T\right],\\{\frac {\partial C_{B}}{\partial t}}=D_{B}\Delta C_{B}-V_{B}\left({\frac {\partial C_{B}}{\partial x}}+{\frac {\partial C_{B}}{\partial y}}\right)-{\frac {\alpha }{{\theta }_{B}}}\left(C_{B}-C_{A}\right)+f_{b},\left(x,y\right)\in {\Omega }_{2},t\in \left[0,T\right],\\C_{A}\left(x,y,0\right)=C_{A0},C_{B}\left(x,y,0\right)=C_{B0},\left(x,y\right)\in {\Omega }_{2},\\C_{A}\left(x,y,t\right)=f_{1},C_{B}\left(x,y,t\right)=f_{2},\left(x,y\right)\in {\Omega }_{2},t\in \left[0,T\right].\end{array}}\end{cases}}}$

Let ${\displaystyle D_{A}=D_{B}=1}$ , ${\displaystyle V_{A}=V_{B}=1}$ , ${\displaystyle {\theta }_{A}={\theta }_{B}=1}$ , ${\displaystyle \alpha =0.001}$ , ${\displaystyle {\Omega }_{2}=\left[0,1\right]\times \left[0,1\right]}$ . The exact solutions are ${\displaystyle C_{A}=e^{1/2\left(x+y\right)-t}}$ , ${\displaystyle C_{B}=\left(t+1\right)x^{2}y^{2}}$ , where ${\displaystyle C_{A0}}$ , ${\displaystyle C_{B0}}$ , ${\displaystyle f_{1}}$ , ${\displaystyle f_{2}}$ are determined by exact solutions. Space step size is ${\displaystyle h=1/(N-1)}$, time step size is ${\displaystyle \Delta t=0.001}$ . The MQ type function ${\displaystyle \phi \left(r\right)={\sqrt {c^{2}+{\left(a\cdot r\right)}^{2}}}}$ as basis function, where the direction of calculation is selected from a unit circle. The direction is:

${\displaystyle a_{v}=cos\left({\theta }_{v}\right)}$ , ${\displaystyle {\theta }_{v}=\left(v-1\right)\pi /\left(m-1\right)}$ , ${\displaystyle v=1,2,\cdots ,m}$

When ${\displaystyle m=9}$ , calculation results are shown in Table 4.

Table 4 The calculation results of the method in this paper (T=1)
 Inserted nodes (N) Space step (h) Parameter (c) Calculation error (e) 11×11 1/10 0.4 e1=7.4021e-05 e2=7.2808e-04 21×21 1/20 0.2 e1=3.7560e-05 e2=3.5269e-04 41×41 1/40 0.1 e1=1.4015e-05 e2=2.1705e-04

As can be seen from table 4 that the meshless method with ridge basis functions for solving the 2D two-flow domain model problem can improve computing efficiency and accuracy. Table 5 gives comparison of the calculation results of the proposed method and the FDM under different node numbers at T=1. Table 6 gives the error comparison of the numerical solutions of the proposed method when selecting different node numbers and calculation directions. As can be seen from table 5, with the increase of the number of nodes, the error accuracy of the method in this paper is higher than the FDM. According to table 6, the factors affecting the numerical accuracy of the meshless method with ridge basis functions include not only the number of nodes inserted, the selection of free parameters, but also the selection of the direction of the basis function.

Table 5 Comparison of error between the method in this paper and FDM
 Time Time space Inserted nodes The method of this paper FDM T=1 0.001 5×5 2.5995e-04 2.8608e-02 8.6418e-04 7.5926e-02 8×8 8.6031e-05 2.7613e-02 7.9283e-04 7.0071e-02 16×16 6.4917e-05 2.5862e-02 4.6672e-04 5.8683e-02 31×31 2.3104e-05 2.2896e-02 2.3731e-04 4.6187e-02

Table6 Comparison of errors of the method in this paper
 Time Inserted nodes Parameter m=11 Calculation error m=13 Calculation error m=15 Calculation error T=1 5×5 c=2.0 e1=2.4564e-04 e1=1.1316e-04 e1=6.3281e-04 e2=8.0200e-04 e2=1.0241e-03 e2=1.0646e-03 11×11 c=0.4 e1=7.1922e-05 e1=7.1233e-05 e1=7.0651e-05 e2=7.0068e-04 e2=7.1778e-04 e2=7.0662e-04 41×41 c=0.1 e1=1.3912e-05 e1=2.4231e-05 e1=1.3042e-05 e2=2.3381e-04 e2=3.4299e-04 e2=2.7869e-04

Fig.2 shows the comparison between the exact solution and the numerical solution of the proposed method when the time is T=1 and the number of nodes is 11×11. It can be found that the numerical solution of the meshless method with ridge basis functions is basically consistent with the exact solution.

Figure 2 Comparison of exact solutions and numerical solutions (N=11×11, T=1)

# 5. Conclusions

In this paper, the meshless method with ridge basis functions for 2D two-flow domain model is prorosed, the existence and uniqueness of its numerical solution is analyzed and numerical examples of 1D and 2D two-flow domain model problems are given. We compare the errors between the proposed method and the finite difference method. The numerical results show that the error precision and convergence of the meshless method with ridge basis functions exceed FDM with the increase of the number of nodes. Therefore, the meshless method with ridge basis functions is an effective numerical method that can be used to solve the two-flow domain model problem.

No.

## References

[1] Marco M., Cianci R., Paladino O., Some analytical solutions for two-dimensional convection-dispersion equation in cylindrical geometry. Environ. Model. Softw., 21:681-688, 2016.

[2] Coats K. H., Smith B. D., Dead-end pore volume and dispersion in porous media. Soc. Petro. Eng. J., 4:73-84, 1964.

[3] Skopp J., Gardner W. R., Tyler E. J., Solute movement in structured soil: Two-region model with small interaction. Soil Sci. Soc. Am. Proc., 45:837-842, 1981.

[4] Zhang X., Liu Y., Ma S., Theory and application of meshless method. Prog. Mech., 39:1-36, 2009.

[5] Wu Z. M., Radial basis function, scattered data fitting and numerical solution of meshless partial differential equation. J. Eng. Math., 19:1-12, 2002.

[6] Reutskiy S. Y., A meshless radial basis function method for 2D steady-state heat conduction problems in anisotropic and inhomogeneous media. Eng. Anal. Bound. Elem., 66:1-11, 2016.

[7] Zabihi F., Saffarian M., A meshless method using radial basis functions for the numerical solution of two-dimensional ZK–BBM equation. Int. J. Appl. Comput. Math., 2016.

[8] Ahmad J., Elyas S., Numerical simulation of nonlinear coupled Burgers’ equation through meshless radial point interpolation method. Eng. Anal. Bound. Elem., 95:187-199, 2018.

[9] Gordon Y., Maiorov V., Meyer M., et al, On the best approximation by ridge function in the uniform norm. Constr. Approx., 18:61-85, 2002.

[10] Ismailov V. E., A note on the best ${\displaystyle L_{2}}$ approximation by ridge function. Appl. Math. E-Notes, 7:71-76, 2007.

[11] Zhang L. W., Error estimates for interpolation with ridge basis function. J. Fudan Univ., Nat. Sci. Ed., 44: 301-306, 2005.

[12] Zhang L. W., Approximation to limited linear combined ridge basis function of plane wave. J. Shandong Univ., Nat. Sci. Ed., 41:87-92, 2006.

[13] Shu H. M., Huang C. Q., Li C. W., A novel meshless method based on ridge basis function. J. Chin. Petro. Univ., 32:108-113, 2008.

[14] Wang Z. G., Qin X. Q., Guo W., et al, Meshless method with ridge basis functions. Appl. Math. Comput., 217:1870-1886, 2010.

[15] Tian S. L., Fang B. Y., Wang Z. G., Approximation of two-point boundary value problems for differential games based on ridge basis function. J. Shandong Univ., Nat. Sci. Ed., 46:38-41, 2011.

### Document information

Published on 12/09/19
Submitted on 02/07/19

Volume 3, 2019

### Document Score

0

Views 45
Recommendations 0