ABSTRACT

The singularly perturbed problem’s high accuracy numerical method is always needed. In this paper, barycentric Lagrange interpolation collocation method (BLICM) is proposed for solving a class of singularly perturbed parabolic reaction-diffusion problem. Compared with other methods, the numerical experiment shows BLICM is a high precision method to solve this class of problems.

Keywords:Singular Perturbation, Barycentric Lagrange Interpolation Collocation Method, Parabolic Reaction-Diffusion Problem

一类奇异摄动抛物型反应扩散问题的数值解

赵辰辉，王玉兰

内蒙古工业大学，内蒙古 呼和浩特

收稿日期：2019年6月11日；录用日期：2019年6月21日；发布日期：2019年6月28日

摘 要

奇异摄动问题一直需要高精度的数值方法，本文提出了重心拉格朗日插值配点法来求解一类奇异摄动抛物型反应扩散问题，与其他方法相比，数值实验表明该方法是解决这类问题的一种高精度方法。

关键词 :奇异摄动，重心拉格朗日插值配点法，抛物型反应扩散问题

Copyright © 2019 by author(s) and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY).

http://creativecommons.org/licenses/by/4.0/

1. 引言

我们考虑下列的奇异摄动抛物型反应扩散问题：

$\{\begin{array}{l}\frac{\partial u}{\partial t}-{\epsilon }_1\frac{{\partial }^{2}u}{\partial x^2}+{\alpha }_{1}u+{\alpha }_{2}v+{\alpha }_{3}r=f_1\left(x,t\right),\\ \frac{\partial v}{\partial t}-{\epsilon }_2\frac{{\partial }^{2}v}{\partial x^2}+{\beta }_{1}u+{\beta }_{2}v+{\beta }_{3}r=f_2\left(x,t\right),\left(x,t\right)\in \left[0,1\right]\times \left[0,1\right],\\ \frac{\partial r}{\partial t}-{\epsilon }_3\frac{{\partial }^{2}r}{\partial x^2}+{\gamma }_{1}u+{\gamma }_{2}v+{\gamma }_{3}r=f_3\left(x,t\right),\end{array}$ (1)

具有初始边界条件

$\begin{array}{l}u\left(x,0\right)=f_0\left(x\right),v\left(x,0\right)=g_0\left(x\right),r\left(x,0\right)=s_0\left(x\right),x\in \left[0,1\right]\\ u\left(0,t\right)=f_1\left(t\right),u\left(1,t\right)=f_2\left(t\right),\\ v\left(0,t\right)=g_1\left(t\right),v\left(1,t\right)=g_2\left(t\right),t\in \left[0,1\right]\\ r\left(0,t\right)=s_1\left(t\right),r\left(1,t\right)=s_2\left(t\right),\end{array}$ (2)

其中 $\epsilon$ 是任意小参数， ${\alpha }_i,{\beta }_i,{\gamma }_i$ 为给定常数， $f_1\left(x,t\right),f_2\left(x,t\right),f_3\left(x,t\right),f_1\left(t\right),f_2\left(t\right),g_1\left(t\right),g_2\left(t\right),s_1\left(t\right),s_2\left(t\right)$ 是已知函数， $u\left(x,t\right),v\left(x,t\right),r\left(x,t\right)$ 分别是未知函数。

奇异摄动抛物型反应扩散问题 [1] [2] [3] [4] [5] 在热科学和力学中有广泛的应用，有一些有价值的努力专注于找到解决问题的分析和数值方法 [6] - [11] 。在 [12] [13] [14] [15] [16] 中，J.P. Berrut引入了重心拉格朗日插值，并研究了它的数值稳定性和收敛性。重心拉格朗日插值在切比雪夫节点处无条件稳定。在 [17] [18] 中，Li和Wang给出了BLICM的算法和程序，用于普通的微分方程，积分方程和差分积分方程以及偏微分方程等。在 [19] - [24] 中，Wang运用了BLICM解决耦合粘性Burgers方程，奇异扰动问题和混沌系统。在 [25] 中，BLICM被提出用于求解第二类线性和非线性高维Fredholm积分方程。参考文献 [17] - [26] 的结果表明BLICM具有高精度，良好的稳定性和收敛性。在本文中，BLICM被提出用于求解一类奇异摄动抛物型反应扩散问题。

2. BLICM的描述

我们考虑规则区域 $\Omega =\left[0,1\right]\times \left[0,1\right]$ ，区间 $\left[0,1\right]$ 被分成M个不同的节点，区间 $\left[0,1\right]$ 也被分成N个不同

的节点。

为了提高计算精度，我们选则Chebyshev节点，切比雪夫节点计算公式如下：

$\begin{array}{l}{x}_i=-\cos \left(\frac{i-1}{M-1}\right)\text{π},i=1,2,\cdots ,M,\\ t_j=-\cos \left(\frac{j-1}{N-1}\right)\text{π},j=1,2,\cdots ,N.\end{array}$ (3)

节点 $\left(x_i,t_j\right)$ 处的 $u\left(x,t\right),v\left(x,t\right)$ 和 $r\left(x,t\right)$ 的重心插值可写为：

$\begin{array}{l}u\left(x,t\right)={\displaystyle \underset{i=1}{\overset{M}{\sum }}{\displaystyle \underset{j=1}{\overset{N}{\sum }}{\xi }_i}\left(x\right){\eta }_j\left(t\right)u\left(x_i,t_j\right)},\\ v\left(x,t\right)={\displaystyle \underset{i=1}{\overset{M}{\sum }}{\displaystyle \underset{j=1}{\overset{N}{\sum }}{\xi }_i}}\left(x\right){\eta }_j\left(t\right)v\left(x_i,t_j\right),\\ r\left(x,t\right)={\displaystyle \underset{i=1}{\overset{M}{\sum }}{\displaystyle \underset{j=1}{\overset{N}{\sum }}{\xi }_i\left(x\right){\eta }_j\left(t\right)r\left(x_i,t_j\right)}},\end{array}$ (4)

其中， ${\xi }_i\left(x\right)=\frac{{\displaystyle \underset{k=1,k\ne i}{\overset{M}{\prod }}\left(x-x_k\right)}}{{\displaystyle \underset{k=1,k\ne i}{\overset{M}{\prod }}\left(x_i-x_k\right)}},i=1,2,\cdots ,M$ ， ${\eta }_j\left(t\right)=\frac{{\displaystyle \underset{k=1,k\ne j}{\overset{N}{\prod }}\left(t-t_k\right)}}{{\displaystyle \underset{k=1,k\ne j}{\overset{N}{\prod }}\left(t_j-t_k\right)}},j=1,2,\cdots ,N$ 。

使用公式(4)，节点上函数 $u\left(x,t\right)$ 的l + k阶偏导数可表示为：

$\begin{array}{l}\frac{{\partial }^{l+k}u}{\partial x^l\partial t^k}={\displaystyle \underset{i=1}{\overset{M}{\sum }}{\displaystyle \underset{j=1}{\overset{N}{\sum }}{\xi }_i^{\left(l\right)}\left(x\right){\eta }_j^{\left(k\right)}\left(t\right)u\left(x_i,t_j\right)}},l,k=1,2,\cdots ,\\ \frac{{\partial }^{l+k}v}{\partial x^l\partial t^k}={\displaystyle \underset{i=1}{\overset{M}{\sum }}{\displaystyle \underset{j=1}{\overset{N}{\sum }}{\xi }_i^{\left(l\right)}\left(x\right){\eta }_j^{\left(k\right)}\left(t\right)v\left(x_i,t_j\right)}},l,k=1,2,\cdots ,\\ \frac{{\partial }^{l+k}r}{\partial x^l\partial t^k}={\displaystyle \underset{i=1}{\overset{M}{\sum }}{\displaystyle \underset{j=1}{\overset{N}{\sum }}{\xi }_i^{\left(l\right)}\left(x\right){\eta }_j^{\left(k\right)}\left(t\right)r\left(x_i,t_j\right)}},l,k=1,2,\cdots ,\end{array}$ (5)

在节点 $\left(x_p,t_q\right)$ 处，偏导数的函数值定义为：

$\begin{array}{l}{u}^{\left(l,k\right)}\left(x_p,t_q\right)=\frac{{\partial }^{l+k}u\left(x_p,t_q\right)}{\partial x^l\partial t^k}={\displaystyle \underset{i=1}{\overset{M}{\sum }}{\displaystyle \underset{j=1}{\overset{N}{\sum }}{\xi }_i^{\left(l\right)}\left(x_p\right){\eta }_j^{\left(k\right)}\left(t_q\right)u\left(x_p,t_q\right),p=1,2,\cdots ,M;q=1,2,\cdots ,N,}}\\ v^{\left(l,k\right)}\left(x_p,t_q\right)=\frac{{\partial }^{l+k}v\left(x_p,t_q\right)}{\partial x^l\partial t^k}={\displaystyle \underset{i=1}{\overset{M}{\sum }}{\displaystyle \underset{j=1}{\overset{N}{\sum }}{\xi }_i^{\left(l\right)}\left(x_p\right){\eta }_j^{\left(k\right)}}}\left(t_q\right)v\left(x_p,t_q\right),p=1,2,\cdots ,M;q=1,2,\cdots ,N,\\ r^{\left(l,k\right)}\left(x_p,t{}_q\right)=\frac{{\partial }^{l+k}r\left(x_p,t_q\right)}{\partial x^l\partial t^k}={\displaystyle \underset{i=1}{\overset{M}{\sum }}{\displaystyle \underset{j=1}{\overset{N}{\sum }}{\xi }_i^{\left(l\right)}}}\left(x_p\right){\eta }_j^{\left(k\right)}\left(t_q\right)r\left(x_p,t_q\right),p=1,2,\cdots ,M;q=1,2,\cdots ,N.\end{array}$ (6)

令 $x^0={\left[x_1,x_2,\cdots ,x_M\right]}^{\text{T}}$ ， $t^0={\left[t_1,t_2,\cdots ,t_N\right]}^{\text{T}}$ 分别定义为x轴和t轴，由张量型节点坐标组成的矩阵X和T

分别定义为：

$X={\left[{\left(x^0\right)}^{\text{T}},{\left(x^0\right)}^{\text{T}},\cdots ,{\left(x^0\right)}^{\text{T}}\right]}^{\text{T}},T=\left[t^0,t^0,\cdots ,t^0\right].$ (7)

$x,t$ 是由矩阵 $X,T$ 拉伸的 $\left(N\times M\right)$ 维列向量：

$x={\left[X_1,X_2,\cdots ,X_{M\times N}\right]}^{\text{T}},t={\left[T_1,T_2,\cdots ,T_{M\times N}\right]}^{\text{T}}.$ (8)

向量 $x,t$ 和向量 $x^0,t^0$ 的分量之间的关系如下：

$X_k=X_{\left(i-1\right)N+j}=x_i,T_k=T_{\left(i-1\right)N+j}=t_j,i=1,2,\cdots ,M;J=1,2,\cdots ,N;k=1,2,\cdots ,M\times N.$ (9)

假设u和 $u^{\left(l,k\right)}$ 是由 $u\left(x_i,t_j\right)$ 和 $u^{\left(l,k\right)}\left(x_i,t_j\right)\left(i=1,2,\cdots ,M;j=1,2,\cdots ,N\right)$ 组成的列向量，即：

$u={\left[u_1,u_2,\cdots ,u_{M\times N}\right]}^{\text{T}},u^{\left(l,k\right)}={\left[u_1^{\left(l,k\right)},u_2^{\left(l,k\right)},\cdots ,u_{M\times N}^{\left(l,k\right)}\right]}^{\text{T}},$

$u_p=u\left(X_p,T_p\right),u_p^{\left(l,k\right)}=u^{\left(l,k\right)}\left(X_p,T_p\right),p=1,2,\cdots ,M\times N.$

$v={\left[v_1,v_2,\cdots ,v_{M\times N}\right]}^{\text{T}},v^{\left(l,k\right)}={\left[v_1^{\left(l,k\right)},v_2^{\left(l,k\right)},\cdots ,v_{M\times N}^{\left(l,k\right)}\right]}^{\text{T}},$

$v_p=v\left(X_p,T_p\right),v_p^{\left(l,k\right)}=v^{\left(l,k\right)}\left(X_p,T_p\right),p=1,2,\cdots ,M\times N.$

$r={\left[r_1,r_2,\cdots ,r_{M\times N}\right]}^{\text{T}},r_p^{\left(l,k\right)}={\left[r_1^{\left(l,k\right)},r_2^{\left(l,k\right)},\cdots ,r_{M\times N}^{\left(l,k\right)}\right]}^{\text{T}},$

$r_p=r\left(X_p,T_p\right),r_p^{\left(l,k\right)}=r^{\left(l,k\right)}\left(X_p,T_p\right),p=1,2,\cdots ,M\times N.$

因此，(4)可写成下列的矩阵形式：

$u^{\left(l,k\right)}=D^{\left(l,k\right)}u,v^{\left(l,k\right)}=D^{\left(l,k\right)}v,r^{\left(l,k\right)}=D^{\left(l,k\right)}r.$ (10)

在上面的公式中， $D^{\left(l,k\right)}=C^{\left(l\right)}\otimes D^{\left(k\right)}$ 是矩阵 $C^{\left(l\right)}$ 和 $D^{\left(k\right)}$ 的Kronecker乘积，也被称为在节点 $\left\{\left(x_i,t_j\right),i=1,2,\cdots ,M;j=1,2,\cdots ,N\right\}$ 的 $\left(l,k\right)$ 阶部分差分矩阵。 $C^{\left(l\right)}$ 是x方向节点上的l阶差分矩阵， $D^{\left(k\right)}$ 是t方向节点上的k阶差分矩阵。因此，(1)可写成以下矩阵形式：

$\left[\begin{array}{ccc}{D}^{\left(0,1\right)}-{\epsilon }_1{D}^{\left(2,0\right)}+{\alpha }_{1}I& {\alpha }_{2}I& {\alpha }_{3}I\\ {\beta }_{1}I& D^{\left(0,1\right)}-{\epsilon }_2{D}^{\left(2,0\right)}+{\beta }_{2}I& {\beta }_{3}I\\ {\gamma }_{1}I& {\gamma }_{2}I& D^{\left(0,1\right)}-{\epsilon }_3{D}^{\left(2,0\right)}+{\gamma }_{3}I\end{array}\right]\left[\begin{array}{l}u\\ v\\ r\end{array}\right]=\left[\begin{array}{l}{f}_1\\ f_2\\ f_3\end{array}\right]$ (11)

将公式(3)代入初始边界条件(2)，我们可以得到关于 $u\left(x,t\right),v\left(x,t\right),r\left(x,t\right)$ 在 $q=1,2,\cdots ,N$ 的初始边界条件的离散方程：

$\begin{array}{l}{\displaystyle \underset{i=1}{\overset{M}{\sum }}{\displaystyle \underset{j=1}{\overset{N}{\sum }}{\xi }_i\left(x_p\right){\eta }_j\left(0\right)u_{ij}}}=f_0\left(x_p\right),\\ {\displaystyle \underset{i=1}{\overset{M}{\sum }}{\displaystyle \underset{j=1}{\overset{N}{\sum }}{{\xi }^{\prime }}_i}}\left(x_p\right){{\eta }^{\prime }}_j\left(0\right)v_{ij}=g_0\left(x_p\right),\\ {\displaystyle \underset{i=1}{\overset{M}{\sum }}{\displaystyle \underset{j=1}{\overset{N}{\sum }}{{\xi }^{″}}_i}}\left(x_p\right){{\eta }^{″}}_j\left(0\right)r_{ij}=s_0\left(x_p\right),\end{array}$ $\begin{array}{l}{\displaystyle \underset{i=1}{\overset{M}{\sum }}{\displaystyle \underset{j=1}{\overset{N}{\sum }}{\xi }_i}}\left(0\right){\eta }_j\left(t_q\right)u_{ij}=f_1\left(t_q\right),\\ {\displaystyle \underset{i=1}{\overset{M}{\sum }}{\displaystyle \underset{j=1}{\overset{N}{\sum }}{{\xi }^{\prime }}_i}}\left(0\right){{\eta }^{\prime }}_j\left(t_q\right)v_{ij}=g_1\left(t_q\right),\\ {\displaystyle \underset{i=1}{\overset{M}{\sum }}{\displaystyle \underset{j=1}{\overset{N}{\sum }}{{\xi }^{″}}_i}\left(0\right){{\eta }^{″}}_j\left(t_q\right)r_{ij}}=s_1\left(t_q\right),\end{array}$ $\begin{array}{l}{\displaystyle \underset{i=1}{\overset{M}{\sum }}{\displaystyle \underset{j=1}{\overset{N}{\sum }}{\xi }_i}}\left(1\right){\eta }_j\left(t_q\right)u_{ij}=f_2\left(t_q\right),\\ {\displaystyle \underset{i=1}{\overset{M}{\sum }}{\displaystyle \underset{j=1}{\overset{N}{\sum }}{\xi }_i{}^{\prime }}}\left(1\right){{\eta }^{\prime }}_j\left(t_q\right)v_{ij}=g_2\left(t_q\right),\\ {\displaystyle \underset{i=1}{\overset{M}{\sum }}{\displaystyle \underset{j=1}{\overset{N}{\sum }}{{\xi }^{″}}_i}\left(1\right){{\eta }^{″}}_j\left(t_q\right)r_{ij}}=s_2\left(t_q\right),\end{array}$ (12)

在本文中，我们用置换法来施加初始边界条件。

3. 数值算例

在此部分中，我们展示了用此数值方法得到的数值结果，从而求出下面例子的近似解。

例1 我们考虑以下具有不连续源项的奇异摄动抛物型反应扩散问题 [8] ：

$\{\begin{array}{l}\frac{\partial u}{\partial t}-{\epsilon }_1\frac{{\partial }^{2}u}{\partial x^2}+5u-tv=f_1\left(x,t\right),\\ \frac{\partial v}{\partial t}-{\epsilon }_2\frac{{\partial }^{2}v}{\partial x^2}-tu+5\left(1+xt\right)v-r=f_2\left(x,t\right),\\ \frac{\partial r}{\partial t}-{\epsilon }_3\frac{{\partial }^{2}r}{\partial x^2}-2u-\left(1+x^2\right)v+\left(7+xt\right)r=f_3\left(x,t\right),\end{array}$ (13)

具有初始边界条件：

$u\left(0,t\right)=v\left(0,t\right)=r\left(0,t\right)=0,u\left(1,t\right)=v\left(1,t\right)=r\left(1,t\right)=0,u\left(x,0\right)=v\left(x,0\right)=r\left(x,0\right)=0,$ (14)

其中

$f_1\left(x,t\right)=\{\begin{array}{l}2+x+t,0\le x\le 0.5,\\ 1,0.5<x\le 1,\end{array}$

$f_2\left(x,t\right)=\{\begin{array}{l}20\le x\le 0.5,\\ t,0.5<x\le 1,\end{array}$ (15)

和

$f_3\left(x,t\right)=\{\begin{array}{l}1+x,0\le x\le 0.5,\\ x^2+t,0.5<x\le 1.\end{array}$

图1和图2分别给出了选取不同的节点M和N时得到的数值解。表1给出了参数N变为2倍时的数值结果，在表2中，得到了参数N除以2，参数 $\Delta t$ 除以2时的数值结果。在这种情况下，观察到一阶均匀收敛。 $\epsilon$ 的值越大，结果越小，当 $\epsilon$ 的值达到某一值时，数值结果是稳定的。从表1和表2我们可以看出，此方法比以前的方法具有更高的精度，我们推断，时间离散化的误差控制了数值格式的全局误差。

图1. 例1的数值解图当 ${\epsilon }_1={10}^{-9},{\epsilon }_2={10}^{-5},{\epsilon }_3={10}^{-3}$ 和 $N=256,M=8$

图2. 例1的数值解图当 ${\epsilon }_1={10}^{-9},{\epsilon }_2={10}^{-5},{\epsilon }_3={10}^{-3}$ 和 $N=128,M=2$

表1. 例1最大点积误差与 ${\epsilon }_1,{\epsilon }_2,{\epsilon }_3$ 均匀收敛速度的比较

表2. 例1最大点积误差与 ${\epsilon }_1,{\epsilon }_2,{\epsilon }_3$ 均匀收敛速度的比较

表3例举了当 $\epsilon$ ，N和t选取不同的值时，实例1的算法所需的迭代次数。可以观察到，当 $\epsilon$ 较大时，迭代次数随着N略微减小；但对于取值小的 $\epsilon$ ，只需要一次迭代。

表3. 例1算法所需的迭代次数

例2：我们考虑以下的时间延迟奇异摄动反应扩散问题 [9] ：

$\{\begin{array}{l}\frac{\partial u}{\partial t}-{\epsilon }_1\frac{{\partial }^{2}u}{\partial x^2}+2.1u-\left(1-x\right)v-\left(1+x\right)r=16x^2{\left(1-x\right)}^2,\\ \frac{\partial v}{\partial t}-{\epsilon }_2\frac{{\partial }^{2}v}{\partial x^2}-xu+\left(1.1+x\right)v-xr=t^2,\\ \frac{\partial r}{\partial t}-{\epsilon }_3\frac{{\partial }^{2}r}{\partial x^2}-\left(2+x\right)u-\left(1-x\right)v-\left(3.1+x\right)r=-16x^2{\left(1-x\right)}^2,\end{array}$ (16)

具有初始边界条件

$u\left(0,t\right)=v\left(0,t\right)=r\left(0,t\right)=0,u\left(1,t\right)=v\left(1,t\right)=r\left(1,t\right)=0.$

我们计算了对不同 $\epsilon$ 值的最大误差，图3和表4展现了相应的结果。与Shishkin的方法相比，通过本方法获得的最大误差更小。这表明此方法具有更高的精确度，可以看到，随着 $\epsilon$ 值的增加，误差变小。

图3. 对 ${\epsilon }_1={10}^{-1},{\epsilon }_2={10}^{-2},{\epsilon }_3={10}^{-3}$ 和 $N=64,M=8$ 时例2的数值解

表4. 比较例2的最大点误差和收敛速度

例3：

$\{\begin{array}{l}\frac{\partial u}{\partial t}-{\epsilon }_1\frac{{\partial }^{2}u}{\partial x^2}+\left(5+xt\right)u-\left(x^2+t^2\right)v-\left(2+xt\right)r=2\left(1-{\text{e}}^{-t}\right)+5t\cos \left(xt\right),\\ \frac{\partial v}{\partial t}-{\epsilon }_2\frac{{\partial }^{2}v}{\partial x^2}-5xt{\text{e}}^{-x}u+\left(6+\cos \left(xt\right)\right)v-\cos \left(xt\right)r=\left(10x+1\right)t\left(x^2+t^2-5\sin \left(xt\right)\right),\\ \frac{\partial r}{\partial t}-{\epsilon }_3\frac{{\partial }^{2}r}{\partial x^2}-3x{\text{e}}^{-t}u-\left(xt+\sin \left(x+t\right)\right)v+\left(5+\left(x+1\right){\text{e}}^{-xt}\right)r=-3t\cos \left(x+t\right).\\ u\left(0,t\right)=u\left(1,t\right)=0,t\in \left[0,1\right],u\left(x,0\right)=0,x\in \left[0,1\right].\end{array}$ (17)

我们选取参数 ${\epsilon }_1={10}^{-8},{\epsilon }_2={10}^{-6},{\epsilon }_3={10}^{-4}$ 和 $N=60,M=64$ ，初始条件 $u\left(x,0\right)=v\left(x,0\right)=r\left(x,0\right)=0$ ，边界条件 $u\left(0,y\right)=v\left(0,y\right)=r\left(0,y\right)=0$ 。图4显示了方程组的数值近似值，为了将我们的算法与另一种经典的方法比较，表5显示了与之前相同的 ${\epsilon }_1,{\epsilon }_2,{\epsilon }_3$ 值的最大误差。从中可以看出，在这种情况下，最大误差的大小与表中相似，但小于表中所示的误差。根据所考虑的问题，使用重心插值配点法得到的误差有时会更好，有时会更差点。

图4. 当 ${\epsilon }_1={10}^{-8},{\epsilon }_2={10}^{-6},{\epsilon }_3={10}^{-4}$ 时例3的数值解

表5. 比较例3的最大误差和收敛阶数

为了比较两种算法的计算成本，我们在表6中显示了几个N和M值所需的CPU时间，选择固定值 $\left({\epsilon }_1={10}^{-8},{\epsilon }_2={10}^{-6},{\epsilon }_3={10}^{-4}\right)$ 作为扩散参数。相对于与经典的方法比较，当n值增加时，我们的方法产生了一个加速。

表6. 例3的CPU时间比较

例4：我们考虑以下二维时间延迟奇异摄动反应扩散问题

$\{\begin{array}{l}\frac{\partial u}{\partial t}-{\epsilon }_1\frac{{\partial }^{2}u}{\partial x^2}+\left(1+0.5xy\right)\frac{\partial u}{\partial x}+{\text{e}}^{x^{2}y}\frac{\partial u}{\partial y}+{\text{e}}^{x+y}\left(1+t\right)u-t\left(x+y\right)v-txr=10t^2\sin \left(x+y\right),\\ \frac{\partial v}{\partial t}-{\epsilon }_2\frac{{\partial }^{2}v}{\partial x^2}+\left(5+x^{2}y\right)\frac{\partial u}{\partial x}+\left(3+\sin \left(x+y\right)\right)\frac{\partial u}{\partial y}-\left(x+y\right)u+\left(1+t\right)\left(3+x+y\right)v-t\sin \left(y\right)r=-5\left(1-{\text{e}}^{-t}\right)\left(x^2+y^2\right),\\ \frac{\partial r}{\partial t}-{\epsilon }_3\frac{{\partial }^{2}r}{\partial x^2}+\left(3-xy\right)\frac{\partial u}{\partial x}+\left(1+x+y\right)\frac{\partial u}{\partial y}-x{y}^{2}u-t\left(\sin x+\sin y\right)v+{\text{e}}^t\left(2+\cos \left(x+y\right)\right)r=-4t{\text{e}}^t\cos \left(xy\right),\end{array}$

具有初始边界条件

$\begin{array}{l}u\left(0,y,t\right)=4y\sin \left(t\right),\\ u\left(x,0,t\right)=4x\sin \left(t\right),\end{array}$ $\begin{array}{l}v\left(0,y,t\right)=0,\\ v\left(x,0,t\right)=0,\end{array}$ $\begin{array}{l}r\left(0,y,t\right)=3\left(1-{\text{e}}^t\right),\\ r\left(x,0,t\right)=3\left(1-{\text{e}}^t\right),\end{array}$

$u\left(x,y,0\right)=v\left(x,y,0\right)=r\left(x,y,0\right)=0,$

$\begin{array}{l}u\left(1,y,t\right)=4\left(1+y\right)\sin \left(t\right),\\ u\left(x,1,t\right)=4\left(x+1\right)\sin \left(t\right),\end{array}$ $\begin{array}{l}v\left(1,y,t\right)=y{t}^2,\\ v\left(x,1,t\right)=x{t}^2,\end{array}$ $\begin{array}{l}r\left(1,y,t\right)=3{\text{e}}^y\left(1-{\text{e}}^t\right),\\ r\left(x,1,t\right)=3{\text{e}}^x\left(1-{\text{e}}^t\right).\end{array}$

为了获得数值解，我们运用相同的思想来分割具有更多元素的系统。

图5显示了当 $\epsilon ={10}^{-4}$ 时例3的数值解，表7分别为第一，第二，第三分量 $\epsilon$ 值的结果。我们可以看到， $\epsilon$ 的值越小，最大误差也越小。由此可见，本文方法具有较好的精度。

图5. 在T = 1， $N=32,M=32$ 时例4的数值解

表7. 例4中对v的最大误差和收敛阶数的比较

4. 结论

在本文中，我们首次使用重心插值配点法来解决奇异摄动问题。该方法的一些计算结果与其他方法的比较表明，该方法具有较高的精度和收敛性。从这篇文章中，我们可以发现我们的方法可以应用于解决这样的人口统计模型。因此，我们可以将此方法扩展到更广泛的领域。本文的所有数值计算结果都是通过数学软件MatlabR2007b给出。

致谢

感谢王玉兰老师的支持与帮助。

文章引用

赵辰辉,王玉兰. 一类奇异摄动抛物型反应扩散问题的数值解
Numerical Solution of a Class of Singularly Perturbed Parabolic Reaction-Diffusion Problem[J]. 应用数学进展, 2019, 08(06): 1181-1191. https://doi.org/10.12677/AAM.2019.86136

参考文献