Modeling of the migration of endothelial cells on bioactive  micropatterned polymers

Thierry Colin; Marie-Christine Durrieu; Julie Joie; Yifeng Lei; Youcef Mammeri; Clair Poignard; Olivier Saut; Thierry Colin; Marie-Christine Durrieu; Julie Joie; Yifeng Lei; Youcef Mammeri; Clair Poignard; Olivier Saut

doi:10.3934/mbe.2013.10.997

Mathematical Biosciences and Engineering

2013, Volume 10, Issue 4: 997-1015. doi: 10.3934/mbe.2013.10.997

Previous Article Next Article

Modeling of the migration of endothelial cells on bioactive micropatterned polymers

1.
Univ. Bordeaux, IMB, UMR 5251, F-33400 Talence
2.
INSERM, IECB, UMR 5248, F-33607 Pessac
3.
Univ. Bordeaux, IECB, UMR 5248, F-33607 Pessac
4.
INRIA, F-33400 Talence
5.
CNRS, IMB, UMR 5251, F-33400 Talence

Received: 01 June 2012 Accepted: 29 June 2018 Published: 01 June 2013
MSC : Primary: 92B05, 92C17.

In this paper, a macroscopic model describing endothelial cellmigration on bioactive micropatterned polymers is presented. It isbased on a system of partial differential equations ofPatlak-Keller-Segel type that describes theevolution of the cell densities. The model is studiedmathematically and numerically. We prove existence and uniquenessresults of the solution to the differential system. We also show thatfundamental physical properties such as mass conservation, positivityand boundedness of the solution are satisfied. The numerical study allows us to show that the modeling results are in good agreement with the experiments.

Keywords:

Citation: Thierry Colin, Marie-Christine Durrieu, Julie Joie, Yifeng Lei, Youcef Mammeri, Clair Poignard, Olivier Saut. Modeling of the migration of endothelial cells on bioactive micropatterned polymers[J]. Mathematical Biosciences and Engineering, 2013, 10(4): 997-1015. doi: 10.3934/mbe.2013.10.997

Related Papers:

[1]	Jin Li, Jinzheng Qu, Xibo Duan, Xiaoning Su . An image encryption algorithm based on heat flow cryptosystems. Networks and Heterogeneous Media, 2023, 18(3): 1260-1287. doi: 10.3934/nhm.2023055
[2]	Xiongfa Mai, Ciwen Zhu, Libin Liu . An adaptive grid method for a singularly perturbed convection-diffusion equation with a discontinuous convection coefficient. Networks and Heterogeneous Media, 2023, 18(4): 1528-1538. doi: 10.3934/nhm.2023067
[3]	Claudio Canuto, Anna Cattani . The derivation of continuum limits of neuronal networks with gap-junction couplings. Networks and Heterogeneous Media, 2014, 9(1): 111-133. doi: 10.3934/nhm.2014.9.111
[4]	Liping Yu, Hans Kleppe, Terje Kaarstad, Svein M. Skjaeveland, Steinar Evje, Ingebret Fjelde . Modelling of wettability alteration processes in carbonate oil reservoirs. Networks and Heterogeneous Media, 2008, 3(1): 149-183. doi: 10.3934/nhm.2008.3.149
[5]	Stefan Berres, Ricardo Ruiz-Baier, Hartmut Schwandt, Elmer M. Tory . An adaptive finite-volume method for a model of two-phase pedestrian flow. Networks and Heterogeneous Media, 2011, 6(3): 401-423. doi: 10.3934/nhm.2011.6.401
[6]	Yaxin Hou, Cao Wen, Yang Liu, Hong Li . A two-grid ADI finite element approximation for a nonlinear distributed-order fractional sub-diffusion equation. Networks and Heterogeneous Media, 2023, 18(2): 855-876. doi: 10.3934/nhm.2023037
[7]	Jérôme Droniou . Remarks on discretizations of convection terms in Hybrid mimetic mixed methods. Networks and Heterogeneous Media, 2010, 5(3): 545-563. doi: 10.3934/nhm.2010.5.545
[8]	Iryna Pankratova, Andrey Piatnitski . Homogenization of convection-diffusion equation in infinite cylinder. Networks and Heterogeneous Media, 2011, 6(1): 111-126. doi: 10.3934/nhm.2011.6.111
[9]	Didier Georges . Infinite-dimensional nonlinear predictive control design for open-channel hydraulic systems. Networks and Heterogeneous Media, 2009, 4(2): 267-285. doi: 10.3934/nhm.2009.4.267
[10]	Jinhu Zhao . Natural convection flow and heat transfer of generalized Maxwell fluid with distributed order time fractional derivatives embedded in the porous medium. Networks and Heterogeneous Media, 2024, 19(2): 753-770. doi: 10.3934/nhm.2024034

Abstract

1. Introduction

Optimal control problems governed by convection-diffusion equations arise in many real-world problems such as wastewater treatment problems ^[1] and air pollution problems ^[2]. Especially it is also an important step toward optimal control problems of the fluid. Such problems have profound application prospects, which simulate some chemical or biological process between different species. It is difficult or impossible to obtain their analytical solutions, and the only way to solve them is to seek approximate solutions by using numerical methods. Thus, the extraction of their numerical algorithms is of great practical significance for scholars. In this paper, we will propose the barycentric interpolation collocation methods for optimal control problems governed by the nonlinear convection-diffusion equation, which can be written as

$\begin{equation} \begin{matrix} \min J\left( v, u \right) = \dfrac{1}{2}\int_{\Omega }{{{\left( v-\hat{{v}} \right)}^{2}}}d\Omega +\dfrac{w }{2}\int_{\Omega }{{{u}^{2}}}d\Omega , \\ s.t\text{ }-\alpha \Delta v+\beta \cdot \nabla v+\gamma v+\phi(v) = f+u, \text{ }v\in \Omega , \\ v = 0, \text{ }v\in \partial \Omega . \\ \end{matrix} \end{equation}$ $

(1.1)

where $ \Delta $ is the Laplace operator, $ \nabla $ is the gradient operator and $ v $ and $ u $ are the state variable and control variable, respectively. $ \hat{v} $ is the desired variable, and $ \phi(v) $ represents the nonlinear term. $ w $ is the so-called regularization parameter. $ J\left(v, u \right) $ stands for the cost function. $ \Omega $ is a bounded domain with the boundary $ \partial \Omega $.

In recent years, a large number of numerical methods have been presented for optimal control problems.

Weng et al. ^[3] presented a stabilized finite element method for optimal control problems governed by a convection dominated diffusion equation. Sandilya and Kumar ^[4] proposed a discontinuous interpolated finite volume approximation of semilinear elliptic optimal control problems. Y$ \rm \ddot{u} $cel et al. ^[5] applied a discontinuous Galerkin method for optimal control problems governed by a system of convection-diffusion PDEs with nonlinear reaction terms. Lu et al. ^[6] presented a priori error estimates of interpolation coefficients mixed finite element methods for semilinear Dirichlet boundary optimal control problems. Lin et al. ^[7] developed Galerkin spectral methods for optimal control problem governed by elliptic equations. Porwal and Shakya ^[8] presented a finite element method for an elliptic optimal control problem with integral state constraints. In ^[9], a bilinear pseudo-spectral method based on Chebyshev polynomials was used for convection-diffusion optimal control problems, and the coupled Sylvester system was solved by some direct and iterative methods. Wang et al. ^[10] investigated a spectral Galerkin approximation of an optimal control problems governed by a fractional advection-diffusion-reaction equation with an integral fractional Laplacian. Casanova et al. ^[11] adopted radial basis function techniques for optimal constrained optimization problems governed by linear convection-diffusion PDEs.

Recently, many researchers have developed a variety of meshless methods ^[12,13,14] for convection-diffusion problems. In particular, the barycentric interpolation collocation is a novel meshless method, which is different from the traditional element-based numerical method. It can effectively avoid the cumulative error caused by the difference scheme. Additionally, it can handle irregular domains. This method has attracted the attention of scholars due to its high precision, fast speed, good stability and convenient program execution. It was extended to solve various PDEs including linear optimal control problems ^[15,16], Sine-Gordon equations ^[17], high-dimensional Fredholm integral equations ^[18], telegraph equations ^[19], viscoelastic wave equations ^[20], the Allen-Cahn equation ^[21], nonlinear parabolic partial differential equations ^[22], etc. Lately, Yi and Yao ^[23] presented error analysis of a barycentric Lagrange interpolation (BLI) collocation scheme.

To the best of our knowledge, there are few studies about using the barycentric interpolation collocation method for the nonlinear optimal control problem described by Eq (1.1). Based on the above work, we proposed two numerical schemes for the nonlinear optimal control in two-dimensions by using BLI collocation or barycentric rational interpolation (BRI) collocation. Moreover, consistency analysis of discrete schemes is presented.

The rest of this paper is organized as follows: First, we obtain the continuous optimality system using Lagrangian multipliers in Section 2. Then, barycentric interpolation collocation methods are briefly introduced in Section 3. Two barycentric interpolation collocation schemes are proposed in Section 4, and we also provide the consistency analysis of discrete schemes. In Section 5, numerical simulations are presented to verify the accuracy and efficiency of the proposed collocation methods. Finally, conclusions are drawn in Section 6.

2. Optimality conditions

In this section, we consider the optimize-then-discretize approach based on the Lagrangian multiplier method for solving Eq (1.1). First, we define the Lagrangian multiplier $ p $ on $ \Omega $ as follows

$\begin{equation} \mathcal{L} = \dfrac{1}{2}{\int_{\Omega }{\left( v-\hat{v}\right)}^{2}}d\Omega \text{+}\dfrac{w }{2}\int_{\Omega }{{{u}^{2}}}d\Omega +\int_{\Omega }{\left( -\alpha \Delta v+\beta \cdot \nabla v+\gamma v+\phi(v)-f-u \right)p}d\Omega +\int_{\partial \Omega }{vp}ds\ . \end{equation}$ $

(2.1)

Then, the state equations are obtained by differentiating the derivative of $ \mathcal{L} $ with respect to $ p $

$\begin{equation} \left\{ \begin{array}{lr} -\alpha \Delta v+\beta \cdot \nabla v+\gamma v+\phi(v) = f+u, \text{ }v\in \Omega, &\\ v = 0, \text{ }v\in \partial \Omega . \\ \end{array} \right. \end{equation}$ $

(2.2)

Similarly, the adjoint equations are gained by calculating the $ \rm Fr\acute{e}chet $ derivative of $ \mathcal{L} $ with respect to $ v $

$\begin{equation} \left\{ \begin{array}{lr} -\alpha \Delta p-\beta \cdot \nabla p+\gamma p+\phi'(v)p = v-\hat{v}, \text{ }p\in \Omega, &\\ p = 0, \text{ }p\in \partial \Omega . \\ \end{array} \right. \end{equation}$ $

(2.3)

Finally, by differentiating with respect to $ u $, the optimality equation is written as follows

$\begin{equation} wu-p = 0, \text{ }on \; \Omega. \end{equation}$ $

(2.4)

Combining Eqs (2.2)–(2.4), Eq (1.1) can be transformed into a continuous optimality system

$\begin{equation} \left\{ \begin{array}{lr} -\alpha \Delta v+\beta \cdot \nabla v+\gamma v+\phi(v) = f+u, \; \; \text{ }v\in \Omega, &\\ v = 0, \text{ }v\in \partial \Omega , \\ -\alpha \Delta p-\beta \cdot \nabla p+\gamma p+\phi'(v)p = v-\hat{v}, \text{ }p\in \Omega, &\\ p = 0, \text{ }p\in \partial \Omega , \\ wu-p = 0, on\; \Omega.\\ \end{array} \right. \end{equation}$ $

(2.5)

3. Barycentric interpolation collocation method

In this section, we mainly introduce BLI and BRI collocation.

3.1. Barycentric Lagrange interpolation collocation method

This part presents a novel meshless method named BLI to solve nonlinear optimal control problems. Suppose that $ m\text{+}1 $ distinct nodes $ {{x}_{k}} $($ k = 0, 1, \cdots, m $) are given together with a set of functional values $ {{h}_{k}} $ at the discrete nodes $ {{x}_{k}} $ correspondingly. Let $ q\left(x \right) $ denote the approximate interpolation polynomial of $ h(x) $, satisfying $ q\left({{x}_{k}} \right) = {{h}_{k}} $ ($ k = 0, 1, \cdots, m $).

As $ q\left(x \right) $ is unique, we can rewrite it in the Lagrange interpolation polynomial form

$\begin{equation} h(x)\approx q(x) = \sum\limits_{k = 0}^{m}{{{L}_{k}}\left( x \right)}{{h}_{k}}, \quad \text{ }k = 0, 1, \cdots , m, \\ \\ \end{equation}$ $

(3.1)

where $ L_{k}(x) $ is the Lagrange basis function, and $ L_{k}(x) = \dfrac{{\prod\limits_{i\ne k}^{m}{(x-{{x}_{i}})}}}{{\prod\limits_{i\ne k}^{m}{({{x}_{k}}-{{x}_{i}})}}}, \text{ }i = 0, 1, \cdots, m $.

Define $ \tilde{L}\left(x \right) = \left(x-{{x}_{0}} \right)\left(x-{{x}_{1}} \right)\cdots \left(x-{{x}_{m}} \right) $; thus, the basis function $ {{L}_{k}}\left(x \right) $ can be expressed as

$\begin{equation} {{L}_{k}}\left( x \right) = \tilde{L}(x)\dfrac{{\omega_{k}}}{x-{{x}_{k}}}, \quad \text{ }k = 0, 1, \cdots , m, \\ \end{equation}$ $

(3.2)

where $ {\omega_{k}} = \dfrac{1}{\prod\limits_{k\ne i}{\left({{x}_{k}}-{{x}_{i}} \right)}} $ is the barycentric weight.

Supposing that $ q(x) = 1 $ in Eq (3.1) and inserting Eq (3.2) into Eq (3.1), the results are shown as

$ \left\{ {

$\begin{array}{*{20}{l}} {1 = \tilde{L}(x)\sum\limits_{k = 0}^m{\dfrac{{\omega_{k}}}{x-{{x}_{k}}}}, }&{}&{}&{}&{{\rm{(3.3a)}}}\\ {q(x) = \tilde{L}(x)\sum\limits_{k = 0}^m{\dfrac{{\omega_{k}}}{x-{{x}_{k}}}}h_{k}.}&{}&{}&{}&{{\rm{(3.3b)}}} \end{array}$ } \right. $

By dividing Eq (3.3b) by Eq (3.3a) and canceling the common polynomial $ \tilde{L}(x) $, the BLI formula for $ v(x) $ can be expressed as

$\begin{equation} q\left( x \right) = \dfrac{\sum\limits_{k = 0}^m \dfrac{\omega_{k}}{x-x_{k}}h_{k}}{\sum \limits_{k = 0}^{M} \dfrac{\omega_{k}}{x-x_{k}}}: = \sum\limits_{k = 0}^m {{\psi _k}(x){h_k}}.\\ \end{equation}$ $

(3.4)

As a matter of fact, the BLI formula has excellent numerical stability when Chebyshev points are selected, so the nodes selected in this article and the corresponding barycentric weights are

$\begin{equation} x_k = cos(\dfrac{k}{m}\pi), \; \; \; k = 0, 1, \cdots, m. \end{equation}$ $

(3.5)

$\begin{equation} \omega_k = (-1)^{k}\delta^{k}, \; \; \; \; {{{\delta }^{k}} = \left\{ {\begin{array}{*{20}{l}} {\dfrac{1}{2}, \; \; \; k = 0\; \rm or \rm \; m, }\\ {1, \; \; \; \mbox{else}.}\\ \end{array}} \right.} \end{equation}$ $

(3.6)

3.2. Barycentric rational interpolation collocation method

In this subsection, we consider the derivation process of the BRI collocation method. Suppose that $ (x_{0}, h_{0}), (x_{1}, h_{1}), \cdots(x_{m}, h_{m}) $ are given $ m+1 $ interpolation nodes of the interval $ [a, b] $, and that the function $ h(x) $ satisfies $ h(x_{k}) = h_{k}, k = 0, 1, 2\cdots m $. Additionally, $ d (0\leq d \leq m) $ is a chosen integer; We define the following BRI formula:

$\begin{equation} r\left( x \right) = \dfrac{\sum\limits_{k = 0}^{m}{\dfrac{{{w }_{k}}}{x-{{x}_{k}}}{{h}_{k}}}}{\sum\limits_{k = 0}^{m}{\dfrac{{{w }_{k}}}{x-{{x}_{k}}}}} = \sum\limits_{k = 0}^{m}{{r_{k}}\left( x \right){{h}_{k}}}, \\ \end{equation}$ $

(3.7)

where $ {{r }_{k}}\left(x \right) = \dfrac{\dfrac{{{w }_{k}}}{x-{{x}_{k}}}}{{\sum\limits_{k = 0}^{m}{\dfrac{{{w }_{k}}}{x-{{x}_{k}}}}}} $ is the basic function and the interpolation weight $ {{w }_{k}} = \sum\limits_{i\in {{J}_{k}}}{{{\left(-1 \right)}^{i}}\prod\limits_{l = i, l\ne k}^{i+d}{\dfrac{1}{{{x}_{k}}-{{x}_{l}}}}} $, $ {{J}_{k}} = \left\{ i\in J|k-d\le i\le k \right\} $, $ J = \left\{ 0, 1, \cdots, m-d \right\} $. More details about the BRI formula are available ^[24,25].

Similarly, the BRI formula for the binary function $ h(x, y) $ with $ r_{M, N}(x, y) $ can be obtained. Let us partition the domain $ \Omega = [a, b]\times [c, d] $ into $ \Omega_{h} = {(x_{k}, y_{j}), k = 0, 1, \cdots M, j = 0, 1, 2, \cdots N} $. Then $ d_{1}\; and \; d_{2} $ are two integers, where $ 0\leq d_{1}\leq M $, $ 0\leq d_{2} \leq N $, $ h(x_{k}, y_{j}) = h_{kj} $. Thus, $ r_{M, N}(x, y) $ can be expressed as

$\begin{equation} h(x, y):\approx r_{M, N}(x, y) = \sum\limits_{k = 0}^{M}\sum\limits_{j = 0}^{N} r_{k}(x)r_{j}(y)h_{kj}, \end{equation}$ $

(3.8)

where $ r_{k}(x) $ and $ r_{j}(y) $ are basis functions, which can be expressed as

$\begin{equation} {{r }_{k}}\left( x \right) = \dfrac{\dfrac{{{\bar{w} }_{k}}}{x-{{x}_{k}}}}{{\sum\limits_{s = 0}^{M}{\dfrac{{{\bar{w} }_{s}}}{x-{{x}_{s}}}}}}, \; {{\bar{w} }_{k}} = \sum\limits_{i\in {{J}_{k}}}{{{\left( -1 \right)}^{i}}\prod\limits_{l = i, l\ne k}^{i+d_{1}}{\dfrac{1}{{{x}_{k}}-{{x}_{l}}}}}, \; J_{k} = \left\{ k\in{0, 1, \cdots , M-d_{1}: k-d_{1}}\leq i\leq k \right\}, \end{equation}$ $

(3.9)

$\begin{equation} {{r }_{j}}\left( y \right) = \dfrac{\dfrac{{{\bar{v} }_{j}}}{y-{{x}_{j}}}}{{\sum\limits_{s = 0}^{N}{\dfrac{{{\bar{v} }_{s}}}{y-{{y}_{s}}}}}}, \; {{\bar{v} }_{j}} = \sum\limits_{i\in {{J}_{j}}}{{{\left( -1 \right)}^{i}}\prod\limits_{l = i, l\ne j}^{i+d_{2}}{\dfrac{1}{{{y}_{j}}-{{y}_{l}}}}}, \; J_{j} = \left\{ i\in{0, 1, \cdots , N-d_{2}: j-d_{2}}\leq i\leq j \right\}. \end{equation}$ $

(3.10)

3.3. Differential matrix

In order to solve the optimal control problem, derivatives expressions of the barycentric interpolation formula should be derived. Obviously, the $ \mu $-order derivative of $ q(x) $ at the different nodes $ x_{i} $ ($ i = 0, 1, \cdots, M $) can be obtained as follows:

$\begin{equation} {q^{\left(\mu\right)}}\left( {{x}_{i}} \right) = \sum\limits_{k = 0}^{M}{\psi _{_{k}}^{\left( \mu \right)}\left( {{x}_{i}} \right){{h}_{k}}}: = \sum\limits_{k = 0}^{M}{D_{ik}^{\left( \mu \right)}{{h}_{k}}}, \; \; \; \text{ }\mu = 1, 2, \cdots, \end{equation}$ $

(3.11)

where $ D_{ik}^{(\mu)} $ denotes the element of the differential matrix $ D^{(\mu)} $, the recursive formula of which ^[26,27] is as follows according to the principles of mathematical induction

$\begin{equation} \left\{ {\begin{array}{*{20}{l}} {D_{ik}^{(1)} = \dfrac{w_{i}}{w_{k}}\dfrac{1}{(x_{i}-{x_{k}})}, \; \; \; i \ne k}, \\ {D_{ii}^{(1)} = -\sum \limits_{k = 0, k \ne i}^M{D_{ik}^{(1)}} }. \end{array}} \right. \end{equation}$ $

(3.12)

$\begin{equation} \left\{ {\begin{array}{*{20}{l}} {D_{ik}^{(\mu)} = \mu\left(D_{ik}^{(1)}D_{ii}^{(\mu-1)}-\dfrac{D_{ik}^{(\mu-1)}}{x_{i}-{x_{k}}}\right), \; \; \; i \ne k}, \\ {D_{ii}^{(\mu)} = -\sum \limits_{k = 0, k \ne i}^M{D_{ik}^{(\mu)}} }. \end{array}} \right. \end{equation}$ $

(3.13)

4. Discrete system and consistency analysis

4.1. Discrete scheme based on barycentric interpolation collocation method

In this part, we will establish numerical schemes for Eq (1.1) by discretizing $ x $ and $ y $ directions with barycentric interpolation formulae. We consider the domain $ \Omega = [0, 1]\times[0, 1] $, and that the $ x $ direction is discretized by $ M+1 $ nodes with $ N+1 $ nodes along the $ y $ direction. Supposing that $ \beta = (\beta_{1}, \beta_{2}) $, the optimality conditions of Eq (1.1) are as follows

$\begin{equation} \left\{ \begin{array}{lr} -\alpha \dfrac{{{\partial }^{2}}v}{\partial {{x}^{2}}}-\alpha \dfrac{{{\partial }^{2}}v}{\partial {{y}^{2}}}+\beta_{1}\dfrac{\partial v}{\partial x}+\beta_{2}\dfrac{\partial v}{\partial y}+\gamma v+\phi(v)-\dfrac{1}{w}p = f, \; \; \text{ }v\in \Omega , \\ v(0, y) = 0, \; \; \; \; \; \; \; v(1, y) = 0, \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \text{ }y\in [0, 1] , \\ v(x, 0) = 0, \; \; \; \; \; \; \; v(x, 1) = 0, \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \text{ }x\in [0, 1] .\\ \end{array} \right. \end{equation}$ $

(4.1)

$\begin{equation} \left\{ \begin{array}{lr} -\alpha \dfrac{{{\partial }^{2}}p}{\partial {{x}^{2}}}-\alpha \dfrac{{{\partial }^{2}}p}{\partial {{y}^{2}}}-\beta_{1}\dfrac{\partial p}{\partial x}-\beta_{2}\dfrac{\partial p}{\partial y}+\gamma p+\phi'(v)p-v = -\hat{v}, \; \; \text{ }p\in \Omega , \\ p(0, y) = 0, \; \; \; \; \; \; \; p(1, y) = 0, \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \text{ }y\in [0, 1] , \\ p(x, 0) = 0, \; \; \; \; \; \; \; p(x, 1) = 0, \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \text{ }x\in [0, 1] .\\ \end{array} \right. \end{equation}$ $

(4.2)

According to the barycentric interpolation formula, $ v(x, y) $ and $ p(x, y) $ can be written as follows

$\begin{equation} v\left( x, y \right) = \sum\limits_{k = 0}^{M}{\sum\limits_{j = 0}^{N}{{C_{k}}\left( x \right){{D }_{j}}\left( y \right)}}{{v}_{kj}}, \text{ }p\left( x, y \right) = \sum\limits_{k = 0}^{M}{\sum\limits_{j = 0}^{N}{{C_{k}}\left( x \right){{D }_{j}}\left( y \right)}}{{p}_{kj}}, \ \end{equation}$ $

(4.3)

where $ C_{k}(x) $ and $ D_{j}(y) $ represent basis functions of the barycentric interpolation method along the $ x $ and $ y $ directions. And $ {{v}_{kj}} = v({{x}_{k}}, {{y}_{j}})\; and\; {{p}_{kj}} = p({{x}_{k}}, {{y}_{j}}). $

Consider the $ \lambda+\mu $ order partial derivative of $ v\left(x, y \right) $ with respect to the spatial variables $ x $ and $ y $; we have

$\begin{equation} \dfrac{{{\partial }^{\lambda +\mu }}v(x, y)}{\partial {{x}^{\lambda }}\partial {{y}^{\mu }}}: = {{v}^{\left( \lambda , \mu \right)}}\left( x, y \right) = \sum\limits_{k = 0}^{M}{\sum\limits_{j = 0}^{N}{C _{k}^{\left( \lambda \right)}}}\left( x \right)D _{j}^{\left( \mu \right)}\left( y \right){{v}_{kj}}, \text{ }\lambda , \mu = 0, 1, 2, \cdots .\\ \end{equation}$ $

(4.4)

Thus, the approximate values of the partial derivatives of $ v(x, y) $ at the node $ \left({{x}_{q}}, {{y}_{s}} \right) $ can be written as

$\begin{equation} {{v}^{\left( \lambda , \mu \right)}}\left( {{x}_{q}}, {{y}_{s}} \right) = \sum\limits_{k = 0}^{M}{\sum\limits_{j = 0}^{N}{C _{k}^{\left( \lambda \right)}}}\left( x_{q} \right)D _{j}^{\left( \mu \right)}\left( y_{s} \right){{v}_{kj}}, \text{ }q = 0, 1, 2, \cdots, M, s = 0, 1, 2, \cdots, N. \\ \end{equation}$ $

(4.5)

Similarly, we consider the $ \lambda+\mu $ order partial derivative of $ p(x, y) $ with respect to the variables $ x $ and $ y $, and its approximate value at the node $ \left({{x}_{q}}, {{y}_{s}} \right) $ is

$\begin{equation} \left\{ \begin{array}{lr} {{p}^{\left( \lambda , \mu \right)}}\left( {x, y} \right) = \sum\limits_{k = 0}^{M}{\sum\limits_{j = 0}^{N}{C _{k}^{\left( \lambda \right)}}}\left( x \right)D _{j}^{\left( \mu \right)}\left( y \right){{p}_{kj}}, \text{ }\lambda , \mu = 0, 1, 2, \cdots, \\ {{p}^{\left( \lambda , \mu \right)}}\left( {{x}_{q}}, {{y}_{s}} \right) = \sum\limits_{k = 0}^{M}{\sum\limits_{j = 0}^{N}{C _{k}^{\left( \lambda \right)}}}\left( {{x}_{q}} \right)D _{j}^{\left( \mu \right)}\left( {{y}_{s}} \right){{p}_{kj}}, \text{ }q = 0, 1, \cdots , M, s = 0, 1, \cdots , N.\\ \end{array} \right. \end{equation}$ $

(4.6)

Selecting a suitable $ \lambda $ and $ \mu $ and then substituting Eqs (4.4)–(4.6) into Eqs (4.1) and (4.2), is implied that

$\begin{equation} \left\{ \begin{array}{lr} -\alpha \sum\limits_{k = 0}^{M}{\sum\limits_{j = 0}^{N}{ C_{k}^{(2)}({{x}_{q}}}})D_{j}^{(0)}({{y}_{s}}){{v}_{kj}}-\alpha \sum\limits_{k = 0}^{M}{\sum\limits_{j = 0}^{N}{C _{i}^{(0)}({{x}_{q}}}})D _{j}^{(2)}({{y}_{s}}){{v}_{kj}}+{{\beta }_{1}}\sum\limits_{k = 0}^{M}{\sum\limits_{j = 0}^{N}{C _{k}^{(1)}({{x}_{q}}}})D _{j}^{(0)}({{y}_{s}}){{v}_{kj}} \\ +{{\beta}_{2}}\sum\limits_{k = 0}^{M}{\sum\limits_{j = 0}^{N}{C _{k}^{(0)}({{x}_{q}}}})D _{j}^{(1)}({{y}_{s}}){{v}_{kj}}+\gamma{{v}_{kj}}+\phi(v_{kj})-\dfrac{1}{w }{{p}_{kj}} = {{f}_{kj}} , \\ -\alpha \sum\limits_{k = 0}^{M}{\sum\limits_{j = 0}^{N}{\phi C_{k}^{(2)}({{x}_{q}}}})D_{j}^{(0)}({{y}_{s}}){{p}_{kj}}-\alpha \sum\limits_{k = 0}^{M}{\sum\limits_{j = 0}^{N}{C _{i}^{(0)}({{x}_{q}}}})D _{j}^{(2)}({{y}_{s}}){{p}_{kj}}-{{\beta }_{1}}\sum\limits_{k = 0}^{M}{\sum\limits_{j = 0}^{N}{C _{k}^{(1)}({{x}_{q}}}})D _{j}^{(0)}({{y}_{s}}){{p}_{kj}} \\ -{{\beta}_{2}}\sum\limits_{k = 0}^{M}{\sum\limits_{j = 0}^{N}{C _{k}^{(0)}({{x}_{q}}}})D _{j}^{(1)}({{y}_{s}}){{p}_{kj}}+\gamma {{p}_{kj}}+\phi'(v_{kj})p_{kj}- {{v}_{kj}} = {{-\hat v}_{kj}} , \\ \end{array} \right. \end{equation}$ $

(4.7)

with the following boundary conditions

$\begin{equation} \left\{ \begin{array}{lr} v_{0j} = v(x_{0}, y_{j}) = 0, \; \; \; \; \; \; \; \; \; v_{Mj} = v(x_{M}, y_{j}) = 0, \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \text{ }j = 1, \cdots N-1, \\ v_{k0} = v(x_{k}, y_{0}) = 0, \; \; \; \; \; \; \; \; \; v_{kN} = v(x_{k}, y_{N}) = 0, \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \text{ }k = 1, \cdots M-1, \\ p_{0j} = p(x_{0}, y_{j}) = 0, \; \; \; \; \; \; \; \; p_{Mj} = p(x_{M}, y_{j}) = 0, \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \text{ }j = 1, \cdots N-1, \\ p_{k0} = p(x_{k}, y_{0}) = 0, \; \; \; \; \; \; \; \; p_{kN} = p(x_{k}, y_{N}) = 0, \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \text{ }k = 1, \cdots M-1.\\ \end{array} \right. \end{equation}$ $

(4.8)

This system can be derived by applying the matrix form

$\begin{equation} \left\{ \begin{array}{lr} diag(-\alpha)\cdot(C^{(2)}\otimes I_{N}+I_{M}\otimes D^{(2)}){\bf V}+diag(\beta_{1})\cdot(C^{(1)}\otimes I_{N}){\bf V}+diag(\beta_{2})\cdot(I_{M}\otimes D^{(1)}){\bf V}\\ +diag(\gamma)\cdot{\bf V}+\phi({\bf V})-\dfrac{1}{w}{\bf P} = {\bf F}, \\ diag(-\alpha)\cdot(C^{(1)}\otimes I_{N}+I_{M}\otimes D^{(2)}){\bf P}-diag(\beta_{1})\cdot(C^{(1)}\otimes I_{N}){\bf P}-diag(\beta_{2})\cdot(I_{M}\otimes D^{(1)}){\bf P}\\ +diag(\gamma)\cdot{\bf P}-{\bf V}+\phi'({\bf V}){\bf P} = -\hat{{\bf V}}, \\ \end{array} \right. \end{equation}$ $

(4.9)

where $ \otimes $ stands for the Kronecker product of the matrix, $ C^{(i)}, D^{(i)} (i = 1, 2) $ are differentiating matrices of order $ i $ and $ {\bf V}, {\bf P}, {\bf F}, and\; \hat{{\bf V}} $ are column vectors, which can be respectively expressed as follows

$\begin{equation} \left\{ \begin{array}{lr} {\bf V} = [v_{00}, \cdots, v_{0N}, v_{10}, \cdots, v_{1N} \cdots, v_{M0}, \cdots, v_{MN}]^{T}, \\ {\bf P} = [p_{00}, \cdots, p_{0N}, p_{10}, \cdots, p_{1N} \cdots, p_{M0}, \cdots, p_{MN}]^{T}, \\ {\bf F} = [f_{00}, \cdots, f_{0N}, f_{10}, \cdots, f_{1N}, \cdots, f_{M0}, \cdots, f_{MN}]^{T}, \\ \hat{{\bf V}} = [\hat{v}_{00}, \cdots, \hat{v}_{0N}, \hat{v}_{10}, \cdots, \hat{v}_{1N}, \cdots, \hat{v}_{M0}, \cdots, \hat{v}_{MN}]^{T}.\\ \end{array} \right. \end{equation}$ $

(4.10)

4.2. Consistency analysis

In this section, we present consistency estimates of proposed schemes based on the approximation properties of the BLI or BRI. For the unknown function $ v(x, y) $, the corresponding barycentric Lagrange interpolation is $ v_{h}(x, y) $, so the error $ \zeta(x, y) $ can be defined as

$\begin{equation} \zeta(x, y) = v(x, y)-v_{h}(x, y). \end{equation}$ $

(4.11)

The approximation properties of BLI have been presented ^[23] as follows

Lemma 1. If $ \zeta\left(x, y\right)\in {{C}^{\left(n+1 \right)}}\left(\Omega \right) $, $ \Omega = [a, b]\times[c, d] $ is a smooth and bounded domain, it holds that

$\begin{equation} \left\{ \begin{array}{lr} |\zeta(x, y)|\leq \parallel v^{\left( n+1 \right)}\parallel _{\infty} ({c_{x}(\dfrac{e {L}_{x}}{2M})^{M}}+{c_{y}(\dfrac{e {L}_{y}}{2N})^{N}}), \\ | \dfrac{\partial \zeta\left( x, y \right)}{\partial x} |\leq \parallel v^{\left( n+1 \right)} \parallel _{\infty} ({\dot{c_{x}}(\dfrac{e {L}_{x}}{2(M-1)})^{M-1}}+{c_{y}(\dfrac{e {L}_{y}}{2N})^{N}}), \\ | \dfrac{\partial \zeta\left( x, y \right)}{\partial y}|\leq \parallel v^{\left( n+1 \right)} \parallel _{\infty} ({c_{x}(\dfrac{e {L}_{x}}{2M})^{M}}+{\dot{c_{y}}(\dfrac{e {L}_{y}}{2(N-1)})^{N-1}}), \\ | \dfrac{\partial^{2} \zeta\left( x, y \right)}{\partial x^2} |\leq \parallel v^{\left( n+1 \right)} \parallel _{\infty} ({\ddot{c_{x}}(\dfrac{e {L}_{x}}{2(M-2)})^{M-2}}+{c_{y}(\dfrac{e {L}_{y}}{2N})^{N}}), \\ | \dfrac{\partial^{2} \zeta\left( x, y \right)}{\partial y^2}|\leq \parallel v^{\left( n+1 \right)} \parallel _{\infty} ({c_{x}(\dfrac{e {L}_{x}}{2M})^{M}}+{\ddot{c_{y}}(\dfrac{e {L}_{y}}{2(N-2)})^{N-2}}), \\ \end{array} \right. \end{equation}$ $

(4.12)

where e is the natural logarithm, $ c_{x}, c_{y}, \dot{c_{x}}, \dot{c_{y}}, \ddot{c_{x}}\; and\; \ddot{c_{y}} $ are positive constants, $ L_{x} = \dfrac{b-a}{2} $ and $ L_{y} = \dfrac{d-c}{2} $.

Suppose that $ v(x_{q}, y_{s}) $ and $ p(x_{q}, y_{s}) $ are the corresponding BLIs of $ v(x, y) $ and $ p(x, y) $, and we define the differential operators $ \mathcal{G}_{1} $ and $ \mathcal{G}_{2} $ as follows

$\begin{equation} \begin{aligned} &\mathcal{G}_{1} v(x, y) = f(x, y), \\ &\mathcal{G}_{2} p(x, y) = -\hat{v}(x, y){, }\\ \end{aligned} \end{equation}$ $

(4.13)

and

$\begin{equation} \begin{aligned} &\lim\limits_{(M, N) \to +\infty} \mathcal{G}_{1} v(x, y) = f(x, y), \\ &\lim\limits_{(M, N) \to +\infty} \mathcal{G}_{2}p(x, y) = -\hat{v}(x, y).\\ \end{aligned} \end{equation}$ $

(4.14)

Next, the consistency analysis of the BLI scheme is presented.

Theorem 1. If $ v(x, y) $, $ p(x, y) $ and $ u(x, y) $ $ \in {{C}^{(n+1) }}\left(\Omega \right) $, $ \Omega = [a, b]\times[c, d] $ and the nonlinear term $ \phi(v) \in W^{2, \infty}(\Omega) $ satisfies the Lipschitz condition, it holds that

$\begin{equation} \left| v\left( x, y \right) -v\left( x_{q}, y_{s} \right) \right|\leq C \parallel v^{(n+1)} \parallel _{\infty} ((\dfrac{e {L}_{x}}{2(M-2)})^{M-2}+ (\dfrac{e {L}_{y}}{2(N-2)})^{N-2}), \\ \end{equation}$ $

(4.15)

$\begin{equation} \left| p\left( x, y \right) -p\left( x_{q}, y_{s} \right) \right|\leq \hat{C} \parallel p^{(n+1)} \parallel _{\infty} ((\dfrac{e {L}_{x}}{2(M-2)})^{M-2}+ (\dfrac{e {L}_{y}}{2(N-2)})^{N-2}), \\ \end{equation}$ $

(4.16)

$\begin{equation} \left| u\left( x, y \right) -u\left( x_{q}, y_{s} \right) \right|\leq \tilde{C} \parallel p^{(n+1)} \parallel _{\infty} ((\dfrac{e {L}_{x}}{2(M-2)})^{M-2}+ (\dfrac{e {L}_{y}}{2(N-2)})^{N-2}).\\ \end{equation}$ $

(4.17)

where $ C, \hat{C}\; and\; \tilde{C} $ are positive constants and $ W^{2, \infty}(\Omega) = \{v\in L^{\infty}(\Omega) \mid \partial^{\alpha}v \in L^{\infty}(\Omega), |\alpha|\leq 2\} $ is the Sobolev space.

Proof. The definitions of the operators $ \mathcal{G}_{1} $ and $ \mathcal{G}_{2} $ lead to the following result

$\begin{equation} \begin{aligned} \mathcal{G}_{1} v(x, y)+\mathcal{G}_{2} p(x, y)& = -\alpha[v_{xx}(x, y)+v_{yy}(x, y)]+\beta_{1}v_{x}(x, y)+\beta_{2}v_{y}(x, y) +\phi(v(x, y))\\ &+(\gamma-1) v(x, y)-\dfrac{1}{w}p(x, y) -\alpha[p_{xx}(x, y)+p_{yy}(x, y)] -\beta_{1}p_{x}(x, y)\\ &-\beta_{2}p_{y}(x, y)+\gamma p(x, y)+\phi'(v(x, y))p(x, y)\\ & = f(x, y)-\hat{v}(x, y).\\ \end{aligned} \end{equation}$ $

(4.18)

The discrete equation corresponding to Eq (4.18) is as follows

$\begin{equation} \begin{aligned} \mathcal{G}_{1} v(x_q, y_s)+\mathcal{G}_{2} p(x_q, y_s)& = -\alpha[v_{xx}(x_q, y_s)+v_{yy}(x_q, y_s)]+\beta_{1}v_{x}(x_q, y_s)+\beta_{2}v_{y}(x_q, y_s) +\phi(v(x_q, y_s))\\&+(\gamma-1) v(x_q, y_s)-\dfrac{1}{w}p(x_q, y_s)-\alpha[p_{xx}(x_q, y_s)+p_{yy}(x_q, y_s)]\\ &-\beta_{1}p_{x}(x_q, y_s)-\beta_{2}p_{y}(x_q, y_s)+\gamma p(x_q, y_s)+\phi'(v(x_q, y_s))p(x_q, y_s)\\ & = f(x_q, y_s)-\hat{v}(x_q, y_s).\\ \end{aligned} \end{equation}$ $

(4.19)

To simplify the analysis process, we first consider the terms related to the state variable $ v(x, y) $. Combining Eq (4.18) with Eq (4.19), we obtain

$\begin{equation} \begin{aligned} \mathcal{G}_{1} [v(x, y)- v(x_q, y_s)]& = -\alpha[v_{xx}(x, y)-v_{xx}(x_q, y_s)]-\alpha[v_{yy}(x, y)-v_{yy}(x_q, y_s)]+\beta_{1}[v_{x}(x, y)-v_{x}(x_q, y_s)]\\ &+\beta_{2}[v_{y}(x, y)-v_{y}(x_q, y_s)]+\phi(v(x, y))-\phi(v(x_q, y_s))+(\gamma-1)[v(x, y)- v(x_q, y_s)]\\ :& = \zeta_{1}+\zeta_{2}+\zeta_{3}+\zeta_{4}+\zeta_{5}+\zeta_{6}, \end{aligned} \end{equation}$ $

(4.20)

where $ \zeta_{1}, \zeta_{2}, \zeta_{3}, \zeta_{4}, \zeta_{5}\; and\; \zeta_{6} $ represent

$\begin{equation} \begin{aligned} &\zeta_{1}: = -\alpha[v_{xx}(x, y)-v_{xx}(x_q, y_s)], \; \; \zeta_{2}: = -\alpha[v_{yy}(x, y)-v_{yy}(x_q, y_s)], \\ &\zeta_{3}: = \beta_{1}[v_{x}(x, y)-v_{x}(x_q, y_s)], \ \; \; \; \; \; \zeta_{4}: = \beta_{2}[v_{y}(x, y)-v_{y}(x_q, y_s)], \\ &\zeta_{5}: = \phi(v(x, y))-\phi(v(x_q, y_s)), \; \; \; \; \; \zeta_{6}: = (\gamma-1)[v(x, y)- v(x_q, y_s)].\\ \end{aligned} \end{equation}$ $

(4.21)

Similarly, we have

$\begin{equation} \begin{aligned} \mathcal{G}_{2} [p(x, y)- p(x_q, y_s)]& = -\alpha[p_{xx}(x, y)-p_{xx}(x_q, y_s)]-\alpha[p_{yy}(x, y)-p_{yy}(x_q, y_s)]-\beta_{1}[p_{x}(x, y)-p_{x}(x_q, y_s)]\\ &-\beta_{2}[p_{y}(x, y)-p_{y}(x_q, y_s)]+(\gamma-\dfrac{1}{w})[p(x, y)-p(x_q, y_s)]\\ &+\phi'(v(x, y))p(x, y)-\phi'(v(x_q, y_s))p(x_q, y_s)\\ :& = R_{1}+R_{2}+R_{3}+R_{4}+R_{5}+R_{6}, \end{aligned} \end{equation}$ $

(4.22)

where $ R_{1}, R_{2}, R_{3}, R_{4}, R_{5}\; and\; R_{6} $ represent

$\begin{equation} \begin{aligned} &R_{1}: = -\alpha[p_{xx}(x, y)-p_{xx}(x_q, y_s)], \; \; \; \; \; \; \; R_{2}: = -\alpha[p_{yy}(x, y)-p_{yy}(x_q, y_s)], \\ &R_{3}: = -\beta_{2}[p_{x}(x, y)-p_{x}(x_q, y_s)], \; \; \; \; \; \; \; \; R_{4}: = -\beta_{2}[p_{y}(x, y)-p_{y}(x_q, y_s)], \\ &R_{5}: = (\gamma-\dfrac{1}{w})[p(x, y)-p(x_q, y_s)], \; \; \; \; R_{6}: = \phi'(v(x, y))p(x, y)-\phi'(v(x_q, y_s))p(x_q, y_s).\\ \end{aligned} \end{equation}$ $

(4.23)

Applying Lemma 1, it holds that

$\begin{equation} \begin{aligned} |\zeta_{1}|& = \alpha\left|v_{xx}(x, y)-v_{xx}(x_{q}, y)+v_{xx}(x_{q}, y)-v_{xx}(x_{q}, y_{s})\right|\\ &\leq\alpha\left[\left|v_{xx}(x, y)-v_{xx}(x_{q}, y)\right|+\left| v_{xx}(x_{q}, y)-v_{xx}(x_{q}, y_{s}) \right| \right]\\ &\leq\tilde{C}_{1}\alpha\parallel v^{(n+1)} \parallel _{\infty} ({\ddot{c_{x}}(\dfrac{e {L}_{x}}{2(M-2)})^{M-2}}+{c_{y}(\dfrac{e {L}_{y}}{2N})^{N}})\\ &\leq {C}_{1}\parallel v^{(n+1)} \parallel _{\infty} ({(\dfrac{e {L}_{x}}{2(M-2)})^{M-2}}+{(\dfrac{e {L}_{y}}{2N})^{N}}), \end{aligned} \end{equation}$ $

(4.24)

$\begin{equation} \begin{aligned} |\zeta_{2}|& = \alpha\left|v_{yy}(x, y)-v_{yy}(x, y_{s})+v_{yy}(x, y_{s})-v_{yy}(x_{q}, y_{s})\right|\\ &\leq C_{2}\parallel v^{(n+1)} \parallel _{\infty} ((\dfrac{e {L}_{x}}{2M})^{M}+(\dfrac{e {L}_{y}}{2(N-2)})^{N-2}), \\ \end{aligned} \end{equation}$ $

(4.25)

$\begin{equation} \begin{aligned} |\zeta_{3}|& = \beta_{1}\left|v_{x}(x, y)-v_{x}(x_{q}, y)+v_{x}(x_q, y)-v_{x}(x_q, y_s)\right|\\ &\leq C_{3}\parallel v^{(n+1)} \parallel _{\infty} ((\dfrac{e {L}_{x}}{2(M-1)})^{M-1}+(\dfrac{e {L}_{y}}{2N})^{N}), \\ \end{aligned} \end{equation}$ $

(4.26)

$\begin{equation} \begin{aligned} |\zeta_{4}|& = \beta_{2}\left|v_{x}(x, y)-v_{x}(x, y_{s})+v_{x}(x, y_{s})-v_{x}(x_q, y_s)\right|\\ &\leq C_{4}\parallel v^{(n+1)} \parallel _{\infty} ((\dfrac{e {L}_{x}}{2M})^{M}+(\dfrac{e {L}_{y}}{2(N-1)})^{N-1}).\\ \end{aligned} \end{equation}$ $

(4.27)

Since $ \phi(v) $ satisfies the Lipschitz condition, assuming a positive constant $ K $, we obtain

$\begin{equation} \begin{aligned} |\zeta_{5}| \leq K \mid v(x, y)-v(x_{q}, y_{s})\mid &\leq C_{5}\parallel v^{(n+1)} \parallel _{\infty} ((\dfrac{e {L}_{x}}{2M})^{M}+(\frac{e {L}_{y}}{2N})^{N}).\\ \end{aligned} \end{equation}$ $

(4.28)

As for the $ \zeta_{6} $, we have

$\begin{equation} \begin{aligned} |\zeta_{6}|&\leq \mid \gamma-1\mid \mid v(x, y)- v(x_q, y_s)\mid &\leq C_{6}\parallel v^{(n+1)} \parallel _{\infty} ((\dfrac{e {L}_{x}}{2M})^{M}+(\dfrac{e {L}_{y}}{2N})^{N}).\\ \end{aligned} \end{equation}$ $

(4.29)

Taking $ C = max\{C_{1}, C_{2}, C_{3}, C_{4}, C_{5}, C_{6}\}, $ and substituting Eqs (4.24)–(4.29) into Eq (4.20), we complete the proof of Eq (4.15).

Similarly, the approximation errors of $ R_{1}, R_{2}, R_{3}\; and\; R_{4} $ can be straightforwardly proved. Now we consider the derivation process of $ R_{5}\; and\; R_{6} $. Applying the Lemma 1, we derive

$\begin{equation} \begin{aligned} |R_{5}| = \left|(\gamma-\dfrac{1}{w})[p(x, y)-p(x_q, y_s)]\right| \leq \tilde{C}_{5}\parallel p^{(n+1)} \parallel _{\infty} ((\dfrac{e {L}_{x}}{2M})^{M}+(\frac{e {L}_{y}}{2N})^{N}).\\ \end{aligned} \end{equation}$ $

(4.30)

As $ \phi(v(x, y)) \in W^{2, \infty}(\Omega) $, we have

$\begin{equation} \begin{aligned} |R_{6}|& = |\phi'(v(x, y))p(x, y)-\phi'(v(x_q, y_s))p(x_q, y_s)|\\ & = |\phi'(v(x, y))p(x, y)-\phi'(v(x_q, y_s))p(x, y)+\phi'(v(x_q, y_s))p(x, y)-\phi'(v(x_q, y_s))p(x_q, y_s)|\\ &\leq|[\phi'(v(x, y))-\phi'(v(x_q, y_s))]p(x, y)|+|\phi'(v(x_q, y_s))[p(x, y)-p(x_q, y_s)]|\\ &\leq|[\phi''(v(x_{\gamma}, y_{\xi}))(v(x, y)-v(x_q, y_s))]||p(x, y)|+\dot{C}_{6}[p(x, y)-p(x_q, y_s)]|\\ &\leq \tilde{C}_{6}\parallel v^{(n+1)} \parallel _{\infty} ((\dfrac{e {L}_{x}}{2M})^{M}+(\frac{e {L}_{y}}{2N})^{N}). \end{aligned} \end{equation}$ $

(4.31)

Finally, substituting the error results of $ R_{1}, R_{2}, R_{3}, R_{4}, R_{5}\; and\; R_{6} $ into Eq (4.22), the proof of Eq (4.16) can be completed.

As for the variable $ u(x, y) $, it can be proved that

$\begin{equation} \begin{aligned} \left| u\left( x, y \right) -u( x_{q}, y_{s} ) \right|& = \dfrac{1}{w}\left| p\left( x, y \right)-p( x_{q}, y_{s} ) \right|\\ &\leq \tilde{C} \parallel p^{(n+1)} \parallel _{\infty} ((\dfrac{e {L}_{x}}{2(M-2)})^{M-2}+ (\dfrac{e {L}_{y}}{2(N-2)})^{N-2}).\\ \end{aligned} \end{equation}$ $

(4.32)

This completes the proof of Eq (4.17).

In the following part, we present the error analysis of Eq (1.1), as solved by the BRI scheme. The compatibility errors of two-dimensional nonlinear optimal control problems solved by the BRI method are as follows based on Theorem 1 in ^[19]

Theorem 2. If $ v(x, y) $, $ p(x, y) $ and $ u(x, y) $ $ \in {{C}^{d_{1}+4}[a, b]}\times{{C}^{d_{1}+4}[c, d]} $, $ v(x_{q}, y_{s}) $ and $ p(x_{q}, y_{s}) $ are corresponding numerical solutions solved by BRI. The nonlinear term $ \phi(v) \in W^{2, \infty}(\Omega) $ satisfies the Lipschitz condition, so we have

$\begin{equation} \left| v\left( x, y \right) -v\left( x_{q}, y_{s} \right) \right|\leq C (h_{x}^{d_{1}-1}+h_{y}^{d_{1}-1}), \\ \end{equation}$ $

(4.33)

$\begin{equation} \left| p\left( x, y \right) -p\left( x_{q}, y_{s} \right) \right|\leq C (h_{x}^{d_{1}-1}+h_{y}^{d_{1}-1}), \\ \end{equation}$ $

(4.34)

$\begin{equation} \left| u\left( x, y \right) -u\left( x_{q}, y_{s} \right) \right|\leq C (h_{x}^{d_{1}-1}+h_{y}^{d_{1}-1}).\\ \end{equation}$ $

(4.35)

where $ C $ is the positive constant independent of the grid size $ h_x $ and $ h_y $; $ h_{x} and h_{y} $ are step sizes of the discretized grids in the $ x $ and $ y $ directions, respectively.

4.3. Fully discrete scheme based on Newton's iteration

The nonlinear term $ \phi(V^{k+1}) $ is expanded by using Taylor's formula at $ V^{k} $, and $ \phi'(V^{k+1})P^{k+1} $ is treated with

$\begin{equation} \begin{aligned} \phi(V^{k+1})&\approx\phi(V^{k})+\phi'(V^{^{k}})(V-V^{K}), \\ \phi'(V^{k+1})&P^{k+1}\approx\phi'(V^{k})P^{k+1}. \end{aligned} \end{equation}$ $

(4.36)

Thus, we derive the numerical scheme of Eq (1.1) based on Eq (4.9) as follows

$\begin{equation} \left\{ \begin{array}{lr} \{diag(-\alpha)\cdot(C^{(2)}\otimes I_{N}+I_{M}\otimes D^{(2)})+diag(\beta_{1})\cdot(C^{(1)}\otimes I_{N})+ diag(\beta_{2})\cdot(I_{M}\otimes D^{(1)})\}{\bf V}^{k+1}\\ +diag(\gamma+\phi'({\bf V}^{k})){\bf V}^{k+1}-\dfrac{1}{w}{\bf P}^{k+1} = -\phi({\bf V}^{k})+\phi'({\bf V}^{k}){\bf V}^{k} +{\bf F}^{k+1}, \\ \{diag(-\alpha)\cdot(C^{(2)}\otimes I_{N}+I_{M}\otimes D^{(2)})-diag(\beta_{1})\cdot(C^{(1)}\otimes I_{N})-diag(\beta_{2})\cdot(I_{M}\otimes D^{(1)})\}{\bf P}^{k+1}\\ +diag(\gamma+\phi'({\bf V}_{k})){\bf P}^{k+1}-{\bf V}_{k+1} = -\hat{{\bf V}}_{k+1}. \\ \end{array} \right. \end{equation}$ $

(4.37)

Rewrite Eq (4.37) in matrix form

$\begin{equation} \left[ \begin{matrix} {\bf A+B+C} & diag{(-\dfrac {1}{w})} \\ -{\bf I} & {\bf A-B+C } \\ \end{matrix} \right]\left[ \begin{matrix} {\bf V}^{k+1} \\ {\bf P}^{k+1} \\ \end{matrix} \right] = \left[ \begin{matrix} -\phi({\bf V}^{k})+\phi'({\bf V}^{k}){\bf V}^{k}+{\bf F}^{k+1} \\ -\hat{{\bf V}}_{k+1} \\ \end{matrix} \right], \end{equation}$ $

(4.38)

where $ \textbf{A} = diag(-\alpha)\cdot(C^{(2)}\otimes I_{N}+I_{M}\otimes D^{(2)} $, $ \textbf{B} = diag(\beta_{1})\cdot(C^{(1)}\otimes I_{N})+diag(\beta_{2})\cdot(I_{M}\otimes D^{(1)})\} $ and $ \textbf{C} = diag(\gamma+\phi'(\textbf{V}^{k})) $.

5. Numerical experiments

We will present two numerical examples to illustrate the accuracy and stability of our proposed methods in this section.

5.1. Problem 1

Consider $ \Omega = [0, 1]\times[0, 1] $; the exact solutions of the state and adjoint functions are as follows

$\begin{equation} \left\{ \begin{array}{lr} v(x, y) = \sin (\pi x)\sin (\pi y), \\ p(x, y) = \pi^{2}\sin (\pi x)\sin (\pi y), \\ \end{array} \right. \end{equation}$ $

(5.1)

with

$ f = (1+2\alpha \pi^{2}-\dfrac{\pi^{2}}{w})\sin (\pi x)\sin (\pi y)+\beta_{1}\pi\cos (\pi x)\sin (\pi y)+\beta_{2}\pi\sin (\pi x)\cos (\pi y) +(\sin (\pi x)\sin (\pi y))^{3}, $

$ \hat{v} = (1-2\alpha \pi^{4}-\pi^{2})\sin (\pi x)\sin (\pi y)+\beta_{1}\pi^{3}\cos (\pi x)\sin (\pi y)+\beta_{2}\pi^{3}\sin (\pi x)\cos (\pi y) -3\pi^{2}(\sin (\pi x)\sin (\pi y))^{3}. $

Take the simulation parameters as follows: $ \alpha = 0.3 $, $ \beta = [1, 1] $, $ \gamma = 1 $, $ w = 0.5 $, $ d_{1} = [\dfrac{3M}{4}] $, $ d_{2} = [\dfrac{3N}{4}] $, and vary the numbers of mesh nodes $ M, N $, $ \phi(v(x, y)) = (v(x, y))^{3}. $ The $ L^\infty $ norm errors of different variables $ v, p, u $ solved by a BLI, BRI or five-point finite difference (FD) scheme are shown in Tables 1–3. The computational results in Tables 1, 2 show that. Obviously, both $ \rm BLI $ and $ \rm BRI $ schemes for Eq (1.1) produce the higher-order accurate solutions when selecting fewer nodes. Furthermore, it is observed that the accuracy of the BLI method is slightly higher than that of the BRI method when taking the same mesh nodes. In Table 3, the obtained results illustrate that the convergence rates of FD scheme are almost second order, and thus are in line with the theoretical analysis. It can be seen that the barycentric interpolation collocation method achieves a high accuracy on the order of $ 10^{-7} $ with $ 8\times 8 $ mesh nodes, while the five-point FD approach can only achieve an accuracy on the order of $ 10^{-4} $ using $ 64\times 64 $ mesh nodes. This shows that the barycentric interpolation collocation method possesses higher accuracy than the FD approach.

Table 1. The $ L^\infty $ errors of different variables $ v, p, u $ solved by the BLI scheme.

$ (N, M) $	$ E _{{\infty}} $($ v $)	$ E _{{\infty}} $($ p $)	$ E _{{\infty}} $($ u $)
(4, 4)	4.3856 $ \times 10^{-3} $	9.0961$ \times 10^{-4} $	9.0961$ \times 10^{-4} $
(8, 8)	1.7191$ \times 10^{-7} $	1.0134$ \times 10^{-7} $	1.0134$ \times 10^{-7} $
(12, 12)	1.6189$ \times 10^{-12} $	1.1602$ \times 10^{-12} $	1.1602$ \times 10^{-12} $
(16, 16)	9.0836$ \times 10^{-14} $	4.8873$ \times 10^{-14} $	4.8873$ \times 10^{-14} $

| Show Table

DownLoad: CSV

Table 2. The $ L^\infty $ errors of different variables $ v, p, u $ solved by the BRI scheme.

$ (N, M) $	$ E _{{\infty}} $($ v $)	$ E _{{\infty}} $($ p $)	$ E _{{\infty}} $($ u $)
(4, 4)	4.3856 $ \times 10^{-3} $	9.0691$ \times 10^{-4} $	9.0691$ \times 10^{-4} $
(8, 8)	2.1043$ \times 10^{-6} $	1.1426$ \times 10^{-6} $	1.1426$ \times 10^{-6} $
(12, 12)	1.4114$ \times 10^{-10} $	9.5459$ \times 10^{-11} $	9.5459$ \times 10^{-11} $
(16, 16)	2.0288$ \times 10^{-13} $	1.0036$ \times 10^{-13} $	1.0036$ \times 10^{-13} $

| Show Table

DownLoad: CSV

Table 3. The $ L^\infty $ errors, convergence rates of different variables $ v, p, u $ solved by FD scheme.

$ (N, M) $	$ E _{{\infty}} $($ v $)	$ \rm Order $	$ E _{{\infty}} $($ p $)	$ \rm Order $	$ E _{{\infty}} $($ u $)	$ \rm Order $
(8, 8)	1.6718 $ \times 10^{-2} $	–	6.6866$ \times 10^{-3} $	–	6.6866$ \times 10^{-3} $	–
(16, 16)	4.2523$ \times 10^{-3} $	1.9751	1.7503$ \times 10^{-3} $	1.9336	1.7503$ \times 10^{-3} $	1.9336
(32, 32)	1.0727$ \times 10^{-3} $	1.9869	4.3889$ \times 10^{-4} $	1.9957	4.3889$ \times 10^{-4} $	1.9957
(64, 64)	2.6842$ \times 10^{-4} $	1.9988	1.0996$ \times 10^{-4} $	1.9969	1.0996$ \times 10^{-4} $	1.9969

| Show Table

DownLoad: CSV

Select the nodes $ M = N = 16 $ to obtain the exact solution images, numerical solution images and error images of $ v $, $ p $, as shown in Figures 1 and 2. We can conclude that numerical solution images solved by the BLI and BRI schemes are all consistent with analytical solution images, indicating that collocation schemes are stable. Figure 3 shows that BLI and BRI schemes have exponential convergence effect when choosing different nodes. In addition, the five-point FD scheme has a second-order convergence rate.

Figure 1. Exact solution images, numerical solution images and error images of $ v $, $ p $ solved by BLI scheme for Example 1 ($ M = N = 16 $).

DownLoad: Full-Size Img PowerPoint

Figure 2. Exact solution images, numerical solution images and error images of $ v $, $ p $ solved by BRI scheme for Example 1 ($ M = N = 16 $).

DownLoad: Full-Size Img PowerPoint

Figure 3. Comparison results of convergence rates in space solved by different schemes for Problem 1.

DownLoad: Full-Size Img PowerPoint

5.2. Problem 2

In this simulation, we consider the following exact solutions

$ v(x, y) = (x^{2}-x)(y^{2}-y)e^{x+y}, $

$ p(x, y) = 4\pi(x^{2}-x)(y^{2}-y)e^{x+y}, $

with

$ \phi(v(x, y)) = (v(x, y))^{3}. $

We consider Eq (1.1) on the computational domain $ \Omega = [0, 1]\times [0, 1] $ with homogeneous Dirichlet boundary condition and set the parameters $ \alpha = 0.1, \beta = [1, 1], \gamma = 1, w = 0.3 $ respectively, $ d_{1} = [\dfrac{3M}{4}] $, $ d_{2} = [\dfrac{3N}{4}] $. The simulation has confirmed the high accuracy and efficiency of the proposed schemes through the calculated errors in Tables 4, 5. Furthermore, the convergence rate of the FD method is almost second order in Table 6.

Table 4. The $ L^\infty $ errors of different variables $ v, p, u $ solved by BLI scheme.

$ (N, M) $	$ E _{{\infty}} $($ v $)	$ E _{{\infty}} $($ p $)	$ E _{{\infty}} $($ u $)
(4, 4)	3.9671 $ \times 10^{-2} $	3.7965$ \times 10^{-3} $	3.7965$ \times 10^{-3} $
(8, 8)	2.6812$ \times 10^{-8} $	7.9556$ \times 10^{-9} $	7.9556$ \times 10^{-9} $
(12, 12)	1.3786$ \times 10^{-13} $	3.9897$ \times 10^{-14} $	3.9897$ \times 10^{-14} $
(16, 16)	8.8305$ \times 10^{-14} $	1.4051$ \times 10^{-14} $	1.4051$ \times 10^{-14} $

| Show Table

DownLoad: CSV

Table 5. The $ L^\infty $ errors of different variables $ v, p, u $ solved by BRI scheme.

$ (N, M) $	$ E _{{\infty}} $($ v $)	$ E _{{\infty}} $($ p $)	$ E _{{\infty}} $($ u $)
(4, 4)	3.9671 $ \times 10^{-2} $	3.7965$ \times 10^{-3} $	3.7965$ \times 10^{-3} $
(8, 8)	1.8400$ \times 10^{-6} $	5.3447$ \times 10^{-7} $	5.3447$ \times 10^{-7} $
(12, 12)	4.3919$ \times 10^{-12} $	1.3874$ \times 10^{-12} $	1.3874$ \times 10^{-12} $
(16, 16)	1.5009$ \times 10^{-13} $	2.2794$ \times 10^{-14} $	2.2794$ \times 10^{-14} $

| Show Table

DownLoad: CSV

Table 6. The $ L^\infty $ errors, convergence rates of different variables $ v, p, u $ solved by FD scheme.

$ (N, M) $	$ E _{{\infty}} $($ v $)	$ \rm Order $	$ E _{{\infty}} $($ p $)	$ \rm Order $	$ E _{{\infty}} $($ u $)	$ \rm Order $
(8, 8)	4.3345 $ \times 10^{-1} $	–	2.4250 $ \times 10^{-2} $	–	2.4250 $ \times 10^{-2} $	–
(16, 16)	1.0308$ \times 10^{-1} $	2.0721	5.7603$ \times 10^{-3} $	2.8271	5.7603$ \times 10^{-3} $	2.8271
(32, 32)	2.5497$ \times 10^{-2} $	2.0154	1.4259$ \times 10^{-3} $	2.0738	1.4259$ \times 10^{-3} $	2.0738
(64, 64)	6.3546$ \times 10^{-3} $	2.0045	3.5574$ \times 10^{-4} $	2.0030	3.5574$ \times 10^{-4} $	2.0030

| Show Table

DownLoad: CSV

Figures 4, 5 show the exact solutions, numerical solutions and error estimates for the state and adjoint functions when selecting Chebyshev nodes $ M = N = 48 $. From these images, it is obvious that two barycentric interpolation collocation schemes are effective and stable with higher accuracy compared with the FD method.

Figure 4. Exact solution images, numerical solution images and error images of $ v $, $ p $ solved by BLI scheme for Problem 2 ($ M = N = 48 $).

DownLoad: Full-Size Img PowerPoint

Figure 5. Exact solution images, numerical solution images and error images of $ v $, $ p $ solved by BRI scheme for Problem 2 ($ M = N = 48 $).

DownLoad: Full-Size Img PowerPoint

Figure 6 illustrates the convergence rates of different schemes and shows that both the BLI and BRI schemes can achieve exponential convergence.

Figure 6. Comparing results for space accuracy as solved by different schemes for Problem 2.

DownLoad: Full-Size Img PowerPoint

6. Conclusion

In this paper, we proposed two barycentric interpolation collocation schemes for optimal control problems governed by the nonlinear convection-diffusion equation. The consistency analyses of the two collocation schemes have been derived. The numerical examples show that the efficiency of our proposed schemes which can achieve the higher-oder accurate solutions with fewer nodes compared with the FD method. Our future work will focus on other kinds of PDE-constrained optimal control problems such as parabolic equations, fluid equations and so on.

Acknowledgments

This work is was supported in part by the NSF of China (No. 11701197) and the Natural Science Foundation of Fujian Province (No. 2022J01308).

Conflict of interest

The authors declare that there is no conflict of interest.

References

[1]	Bulletin of Mathematical Biology, 60 (1998), 857-899.
[2]	Acta Biomaterialia, 6 (2010), 3824-3846.
[3]	I. Colloq. Math., 66 (1993), 319-334.
[4]	Electron. J. Differential Equations, 2006, No. 44, 32 pp. (electronic).
[5]	Nature, 436 (2005), 193-200.
[6]	Science, 276 (1997), 1425-1428.
[7]	in Vitro Cell. Dev. Biol., 35 (1999), 441-448.
[8]	C. R. Math. Acad. Sci. Paris, 339 (2004), 611-616.
[9]	Handbook of Numerical Analysis, (eds. P. G Ciarlet and J. L Lions), 2007.
[10]	Annu. Rev. Biomed. Eng., 2 (2000), 227-256.
[11]	Nature, 288 (1980), 551-556.
[12]	Math. Nachr., 195 (1998), 77-114.
[13]	J. Math. Biol., 58 (2008), 183-217.
[14]	Nonlinear Differ. Equ. Appl., 8 (2001), 399-423.
[15]	Springer Series in Comput. Math., 33, Springer, 2003.
[16]	Biomaterials, 20 (1999), 2333-2342.
[17]	Nat. Med., 9 (2003), 685-593.
[18]	Nat Biotechnol, 23 (2005), 821-823.
[19]	J. of Computational Physics, 126 (1996), 202-228.
[20]	in vivo, Nature, 442 (2006), 453-456.
[21]	Journal of Theo. Biol., 30 (1971), 235-248.
[22]	Small, In Revision.
[23]	PLoS ONE, 7 (2012), e41163.
[24]	Journal of Computational Physics, 115 (1994), 200-212.
[25]	Cell, 112 (2003), 19-28.
[26]	Tissue Eng., 12 (2006), 1143-50.
[27]	Biosens. Bioelectron, 14 (1999), 317-325.
[28]	Blood Cells Mol. Dis., 39 (2007), 212-220.
[29]	London Math. Soc. Monographs Series, Princeton University Press. 31, 2005.
[30]	Bull. Math. Biophys., 15 (1953), 311-338.
[31]	Curr. Opin. Biotechnol, 21 (2010), 704-709.
[32]	Macromol Biosci., 10 (2010), 12-27.
[33]	Adv. Differential Equations, 6 (2001), 21-50.
[34]	Arch. Rational Mech. Anal., 153 (2000), 91-151.
[35]	To Appear in AMS. Proc., 141 (2013), 1067-1081. (Avalaible on http://arxiv.org/abs/1009.1965v2).

This article has been cited by:

Mengli Yao, Zhifeng Weng, A Numerical Method Based on Operator Splitting Collocation Scheme for Nonlinear Schrödinger Equation, 2024, 29, 2297-8747, 6, 10.3390/mca29010006

Reader Comments

Your name:*

Email:*
© 2013 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)

通讯作者: 陈斌, bchen63@163.com

1.
沈阳化工大学材料科学与工程学院沈阳 110142

Mathematical Biosciences and Engineering

4.4

Metrics

Article views(3037) PDF downloads(461) Cited by(7)

Preview PDF

Download XML

Export Citation

Article outline

Show full outline

Mathematical Biosciences and Engineering

Modeling of the migration of endothelial cells on bioactive micropatterned polymers

Related Papers:

Abstract

1. Introduction

2. Optimality conditions

3. Barycentric interpolation collocation method

3.1. Barycentric Lagrange interpolation collocation method

3.2. Barycentric rational interpolation collocation method

3.3. Differential matrix

4. Discrete system and consistency analysis

4.1. Discrete scheme based on barycentric interpolation collocation method

4.2. Consistency analysis

4.3. Fully discrete scheme based on Newton's iteration

5. Numerical experiments

5.1. Problem 1

5.2. Problem 2

6. Conclusion

Acknowledgments

Conflict of interest

References

This article has been cited by:

Reader Comments

通讯作者: 陈斌, bchen63@163.com

Metrics

Other Articles By Authors

Catalog

Mathematical Biosciences and Engineering

Modeling of the migration of endothelial cells on bioactive micropatterned polymers

Related Papers:

Abstract

1. Introduction

2. Optimality conditions

3. Barycentric interpolation collocation method

3.1. Barycentric Lagrange interpolation collocation method

3.2. Barycentric rational interpolation collocation method

3.3. Differential matrix

4. Discrete system and consistency analysis

4.1. Discrete scheme based on barycentric interpolation collocation method

4.2. Consistency analysis

4.3. Fully discrete scheme based on Newton's iteration

5. Numerical experiments

5.1. Problem 1

5.2. Problem 2

6. Conclusion

Acknowledgments

Conflict of interest

References

This article has been cited by:

Reader Comments

通讯作者: 陈斌, bchen63@163.com

Metrics

Other Articles By Authors

Related pages

Tools

Export File

Citation

Format

Content

Catalog