1. 引言
在本文中,我们考虑如下奇异摄动两点边值问题:
(1)
其中,
是一个小的正参数(
),称为摄动参数。这里
、
和
是给定的关于
的函数,且在
上具有光滑性。上述边值问题在给定的边界条件下存在唯一的解。随着
趋于0,这个解
会产生边界层现象[1],即解在过渡层和边界层出现大幅度变化。该问题的解具有多尺度特性[2],即解在边界层内变化迅速,在远离边界层的区域变化相对缓慢,这使得难以用统一的解析表达式或数值方法准确求解整个区域的解。在数值求解时,由于边界层内解的快速变化,容易导致数值方法的不稳定,出现数值振荡或发散等问题。因此,研究解在边界层上的变化显得尤为重要。如何有效地解决边界层现象成为了数值计算的主要困难之一。
许多研究学者构造了求解此类奇异摄动问题的数值方法。Shishkin [3]采用局部加密的方法首先构建了分段均匀网格,利用该网格可以得到拟最优一致误差估计。分段均匀网格简单易于实现,但在处理奇异摄动问题时,可能导致数值解不稳定和解的精度低的问题。而Chebyshev点作为配置点在求解奇异摄动问题时,往往能提供更高的精度,可以有效捕捉到边界层的变化。Kumar和Srinivasan [4]提出了一种新的自适应网格策略,他们使用二阶中心差分格式求解对流扩散奇异摄动问题。自适应网格策略具有极大的灵活性和适应性,能够根据解的特性动态优化网格,适合处理奇异摄动问题,但自适应网格的生成和调整需要复杂的算法和大量的计算开销。与自适应网格策略相比,Chebyshev点作为配置点的方法算法较为简单,计算开销较小。Doolan等[5]和Marusic [6]提出了各种指数拟合有限差分格式。指数拟合有限差分格式是一种通过使用指数类型的插值或逼近来改善在解存在奇异点或边界层时的数值性能方法,但该方法的实现通常较为复杂,在奇异性较强时,可能无法提供足够的精度。相较于复杂且实现较难的指数拟合有限差分格式,Chebyshev点的简单实现和高效稳定的特性使其在处理有奇异摄动的问题时更具吸引力。Schatz和Wahlbin [7]使用有限元法来解决这类问题。赵辰辉等[8]提出了重心拉格朗日插值配点法来求解奇异摄动抛物型反应扩散问题。该方法选用Chebyshev点为配置点,与其他方法相比,该方法具有高精度的特点。
B样条配置法是一种在数值计算和函数逼近等领域广泛应用的方法。该方法以B样条函数作为基函数来构造逼近函数。而配置法则在求解区域内选择一系列特定的点作为配置点,要求在这些点上满足微分方程。目前,已有较多学者采用B样条配置法来求解微分方程。如Khuri和Sayfy [9]提出了一种三次B样条配置法,用于求解线性和非线性二阶边值问题的扩展系统的数值解。韩丹夫等[10]采用改进的三次B样条配置法和Crank-Nicolson格式求解BBMB方程。王仕良[11]利用三次B样条配置法求解Burgers-BBM方程的初边值问题。刘兴霞等[12]运用五次B样条配置法研究了广义KdV方程的孤立波解。以上数值结果表明:B样条配置法是求解微分方程的有效工具。
由于Chebyshev点在区间内是非均匀分布的,具有“两端密、中间疏”的特点。这种分布特性与奇异摄动问题解的特性非常吻合,即:在边界层内,解变化剧烈,它可由Chebyshev点中“密”的部分进行捕捉;另一方面,在远离边界层的区域,解变化相对缓慢,它可由Chebyshev点中“疏”的部分进行计算。鉴于以上原因,本文构造以Chebyshev点作为配置点的三次B样条配置法,用以求解奇异摄动两点边值问题,并探讨该方法的有效性。
2. 数值算法
2.1. Chebyshev点
Chebyshev点的通式[13]为:
(2)
本文将区间
划分为n份,划分节点依次为
,进而可以得出Chebyshev点的步长:
,其中
,
。同时假设:
,
,
,
。
2.2. 非均匀网格下的三次B样条函数
非均匀B样条的递推形式[14]为:
(3)
其中,
;
。
由此可推导出非均匀网格下的三次B样条函数形式[14]为(4)式:
(4)
其中
,进一步得,
(5)
(6)
2.3. B样条配置法求解两点边值问题
考虑如下(7)式的两点边值问题:
(7)
下面介绍B样条配置法。
将(7)式的解记为
,则其可写成(8)式[15]:
(8)
其中,节点
为节点
处的B样条基函数。这里,
的取值为−1到
,以确保端点处使用相邻节点处的B样条值。由此可得[15]:
(9)
(10)
(11)
配置法要求在配置点上满足方程,故有
(12)
上式中,
、
、
由式(9)~(11)代入。
选用(2)式所示的
个Chebyshev点作为配置点,也即满足如下(13)式的线性方程组:
(13)
此外,由边界条件可知:
(14)
(15)
即:
,
。
由此,式(13)可表示为:
(16)
其中:
,
,
。这里
,
,
。
求解式(16),可得系数
,进一步地,将其代入(8)式,即可得到近似函数
。
3. 数值求解
本文用最大误差来衡量算法的准确性。最大误差的定义如下:
(17)
这里,
代表
处的真实解;
代表
处的数值解。
为了说明本文构造的基于Chebyshev点的B样条配置法的有效性,该方法的数值结果与基于等距节点下的B样条配置法的数值结果进行了比较。本文的算例在单机上运行,机器配置为:处理器型号是12th Gen Intel(R) Core(TM) i5-12450H 2.00 GHz,机带RAM 16.0 GB (15.7 GB可用)。
3.1. 奇异摄动扩散反应问题
考虑如下的奇异摄动扩散反应方程:
(18)
其中,
是摄动因子。该问题的解析解由Matlab的dsolve命令求得。
图1为摄动参数
和剖分份数为1000时,基于Chebyshev点和等距节点的B样条配置法所得解与解析解的比较。图2是对图1左端边界层附近的局部放大。由图2可知,等距节点下的数值解在边界层附近表现出数值振荡;而基于Chebyshev点的数值解能够很好地抑制这种现象。图3给出了摄动参数
和剖分份数为1000时,Chebyshev点和等距节点下的误差以及最大误差。等距节点下数值解的最大误差为0.2414,而Chebyshev点下数值解的最大误差为0.0017127。由此可知,基于Chebyshev点的B样条配置法所得的数值解更为精确。
Figure 1. Comparison of three solutions to Example 1 when n = 1000
图1. 当n = 1000时算例1三种解的比较
Figure 2. Local amplification at the left end of Figure 1
图2. 图1左端的局部放大
Figure 3. Error diagram between the isometric node and the Chebyshev node in Example 1
图3. 算例1中等距节点与Chebyshev节点下的误差图
表1和表2分别给出了摄动参数为
和
时,在给定的最大误差限下,基于等距与Chebyshev点的B样条配置法所对应的剖分份数与运行时间。需要注意的是,当摄动参数
和最大误差小于等于10−3时,等距情况要求剖分的份数远远大于14,000份。限于单机的运行能力,本文没有给出具体数值。但从已有的数据来看,在相同最大误差下,基于Chebyshev点的B样条配置法需要剖分的份数远少于等距节点的剖分份数,其运行时间也更短。因此,以Chebyshev点作为配置点的方法克服了由等距节点造成的计算量的猛增和不必要的计算成本的投入。
Table 1. Section number and running time of
in Example 1
表1. 算例1中
的剖分份数和运行时间
最大误差 |
|
|
|
|
等距 |
Chebyshev |
等距 |
Chebyshev |
等距 |
Chebyshev |
剖分份数 |
1270 |
140 |
3940 |
430 |
12,390 |
1330 |
运行时间 |
0.030132 s |
0.005108 s |
0.223570 s |
0.010497 s |
1.351258 s |
0.033031 s |
Table 2. Section number and running time of
in Example 1
表2. 算例1中
的剖分份数和运行时间
最大误差 |
|
|
|
|
等距 |
Chebyshev |
等距 |
Chebyshev |
等距 |
Chebyshev |
剖分份数 |
14,000 |
420 |
|
1320 |
|
4160 |
运行时间 |
1.823096 s |
0.009472 s |
|
0.035149 s |
|
0.238262 s |
3.2. 奇异摄动对流扩散反应问题
考虑如下的奇异摄动对流扩散反应方程:
(19)
其中,
是摄动因子。该问题的解析解由Matlab的dsolve命令求得。
Figure 4. Comparison of the three solutions to Example 2 when n = 1000
图4. 当n = 1000时算例2三种解的比较
Figure 5. Local amplification at the left end of Figure 4
图5. 图4左端的局部放大
Figure 6. Error diagram between the isometric node and the Chebyshev node in Example 2
图6. 算例2中等距节点与Chebyshev节点下的误差图
图4为摄动参数
和剖分份数为1000时,基于Chebyshev点和等距节点的B样条配置法所得解与解析解的比较。图5是对图4左端边界层附近的局部放大。由图5可知,等距节点下的数值解在边界层附近有很强的数值振荡;而基于Chebyshev点的数值解则表现得很稳定。图6给出了摄动参数
和剖分份数为1000时,Chebyshev点和等距节点下的误差以及最大误差。等距节点下数值解的最大误差为1.81109,而Chebyshev点下数值解的最大误差为0.038134。由此可知,基于Chebyshev点的B样条配置法所得的数值解更为精确。
Table 3. Section number and running time of
in Example 2
表3. 算例2中
的剖分份数和运行时间
最大误差 |
|
|
|
|
等距 |
Chebyshev |
等距 |
Chebyshev |
等距 |
Chebyshev |
剖分份数 |
1090 |
220 |
3980 |
2010 |
20000 |
3200 |
运行时间 |
0.019321 s |
0.003833 s |
0.170831 s |
0.072157 s |
3.796157 s |
0.188011 s |
Table 4. Section number and running time of
in Example 2
表4. 算例2中
的剖分份数和运行时间
最大误差 |
|
|
|
|
等距 |
Chebyshev |
等距 |
Chebyshev |
等距 |
Chebyshev |
剖分份数 |
9990 |
630 |
|
2040 |
|
6660 |
运行时间 |
0.780684 s |
0.008838 s |
|
0.059052 s |
|
0.566314 s |
表3和表4分别给出了摄动参数为
和
时,在给定的最大误差限下,基于等距与Chebyshev点的B样条配置法所对应的剖分份数与运行时间。这里,当摄动参数
和最大误差小于等于10−3时,等距情况要求剖分的份数远大于15,000份。限于单机的运行能力,本文没有给出具体数值。由表3和表4可知,在相同最大误差下,基于Chebyshev点的B样条配置法需要剖分的份数远少于等距节点的剖分份数,其运行时间也更短。因此,在该例中,以Chebyshev点为配置点的方法其优越性也高于基于等距节点的方法。
4. 理论解释
基于Chebyshev点的配置法在数值计算中具有优越性,特别是在处理插值、逼近和求解微分方程等问题时。这种优越性可以从最小化插值误差、收敛性、减少数值振荡、全局特性与高效性和数值稳定性与算法五个方面解释。
4.1. 最小化插值误差
Chebyshev点是基于Chebyshev多项式生成的,这些多项式具有特定的正交性质,能有效减小插值时的误差。在Chebyshev点上进行插值时,插值多项式的误差能够均匀分布在整个区间,尤其是在区间的边缘区域。与等距节点相比,Chebyshev点的分布可以有效减小插值中的“龙格现象”。
Chebyshev多项式的正交性意味着它们在某种程度上能够避免相互干扰,从而提供更精确的插值结果。
4.2. 收敛性
基于Chebyshev点的配置法具有良好的收敛性,使用Chebyshev点进行多项式插值或逼近时,其收敛速度通常优于使用均匀分布点。在某些情况下,误差能够以更高的速率趋于零。这是因为Chebyshev点使得高阶多项式在插值过程中的符号变化最小化。
在多项式插值中,由于Chebyshev点的特殊位置,接近真实函数的能力大大增强,使得相对误差而言,在较少的节点数下也可实现较高的精度。
4.3. 减少数值振荡
在处理高阶多项式时,特别是在插值的过程中,Chebyshev点能够有效地减少高频振荡的产生。这在数值计算中是至关重要的,因为高频振荡通常会导致不稳定的数值结果。
Chebyshev点的选择使得在求解过程中不容易受到算法不稳定性的影响,如条件数大的问题,这有助于保持算法的稳定性。
4.4. 全局特性与高效性
使用Chebyshev点,当节点数增加时,相对错误减小的速度非常快,使得其在追求高精度的同时,不会造成过多的计算量。
Chebyshev点的使用使得在多维插值和多项式逼近中,可以通过较少的节点实现较高的精度,从而降低了计算复杂性。
4.5. 数值稳定性与算法
Chebyshev点配置法提高了数值计算过程的稳定性。对于求解偏微分方程或高阶导数问题,增加节点的同时不至于导致数值错误显著增加,保持了算法的鲁棒性。
5. 结论
本文推导了非等距节点下的B样条函数,详细给出了B样条配置法求解两点边值问题的步骤和实施细节。通过奇异摄动扩散反应问题、奇异摄动对流扩散反应问题的研究,表明该算法的有效性、高精度和高效率。
本文提出的方法很好地克服了由均匀网格造成的计算资源的浪费,并且更高分辨率地在边界层捕捉到稳定的数值解。该方法也有望用于其他微分方程的求解。
基金项目
河南科技大学大学生创新创业训练计划项目(项目编号:2024235)。
NOTES
*通讯作者。