基于MATLAB对一维奇异摄动方程高精度格式的研究

DOI：10.19850/j.cnki.2096-4706.2024.02.022

收稿日期：2023-06-24

基金项目：宁夏回族自治区大学生创新项目（S2022-11407-036）

摘  要：为研究更加精确的求解奇异摄动方程数值解的解法，文章研究传统中心差分格式，迎风格式以及高阶精度格式这三种数值求解格式，使用MATLAB软件进行求解一维奇异摄动对流扩散方程的算例并可视化，对不同方法下求解方程的结果进行了对比和研究，通过对数值算例求解结果的分析，发现在任意网格下，当摄动参数趋于0时，高阶精度格式下解的精度和稳定性较传统的求解格式有显著提升，使在求解奇异摄动方程的方法选择上更具有针对性。

关键词：奇异摄动方程；对流扩散方程；MATLAB；高精度格式

中图分类号：TP391；O241.82    文献标识码：A    文章编号：2096-4706（2024）02-0102-07

Research on One-dimensional Singular Perturbation Equation High Precision Scheme Based on MATLAB

PAN Yihang

（School of Mathematics and Information Science， North Minzu University， Yinchuan  750021， China）

Abstract： In order to study a more accurate numerical solution of singular perturbation equation， this paper studies three numerical solutions of traditional central difference scheme， upwind scheme and high order precision scheme. It uses MATLAB software to solve the one-dimensional singularly perturbed convection-diffusion equations for example and visualization， compares and studies the results of solving equations with different methods. Through the analysis of the results of numerical examples， it is found that under any grid， when the perturbation parameter approaches 0， the accuracy and stability of solution with high order precision scheme are significantly improved compared with traditional solution schemes， which makes the method selection of solving singular perturbation equation be more targeted.

Keywords： singular perturbation equation; convection diffusion equation; MATLAB; high precision scheme

0  引  言

奇異摄动方程在数学、物理学和工程技术等领域中均得到了广泛的应用。这类问题通常在高阶导数项包含有一个正的小参数，称为摄动参数。当摄动参数向零无限趋近时，问题的解会在某些很小的区域内发生非常剧烈的非物理震荡，此时解就出现了内点层或边界层。这使得求解奇异摄动问题的解析解工作变得十分困难。所以，我们通常研究它的近似解的数值方法。伴随着计算机的发展，数值解法的研究大大受益于计算机快速的矩阵运算，受到了科研工作者们的重视。

MATLAB软件因具有强大的矩阵运算能力而收到研究者的喜爱，近年来，针对使用MATLAB求解偏微分方程进行了大量的研究[1-3]；同时，对求解对流扩散方程的紧致格式也有大量的研究[4-8]，Lele[9]构造了一类高精度的中心型线性紧致格式；傅德薰和马延文等[10]将迎风机制引入到了紧致差分格式，基于以上研究，本文构造了一种高精度的数值求解格式。最后通过MATLAB软件，使用三种不同的求解格式，实现对一维奇异摄动对流扩散方程两点边值问题的数值模拟，进而研究三种格式对方程的求解情况并在文中给出了详细的数值求解方法与算法步骤。

1  奇异摄动对流扩散方程

1.1   数学模型

本文考虑如下一维奇异摄动对流扩散两点边值问题：

（1）

其中：0＜ε ? 1是摄动参数，a（x）≥a0＞0，a（x）与f （x）均是足够光滑的函数。

1.2   网格剖分

首先将空间区域[0，1]网格剖分为N个等距子区间：0 = x1，x2，…，xN-1，xN，xN+1 = 1，定义空间步长h = 1 / N，xi = 0+（i-1）Δx，i = 1，2，…，N。文中用ui，， 分别表示u在点xi的数值解与一、二阶导数近似值u（xi），u′（xi），u″（xi）分别表示u在点xi处的精确解与一、二阶导数精确值，即ui ≈ u（xi）， ≈ u′（xi）， ≈ u″（xi）。

2  数值方法

2.1  中心差分格式

根据Taylor展开式：

（2）

（3）

可得到各阶导数的近似。

一阶导数的中心差商[11]：

（4）

二阶导数的中心差商：

（5）

扩散项内点使用二阶导数中心差商替代，对流项内点则使用一阶导数中心差商替代，即可得到传统的中心差分格式。

2.2  迎风格式

同样，由式（2）（3）组合可得到向前差商[11]：

（6）

和向后差商：

（7）

扩散项内点仍由二阶导数中心差商替代，对流项内点则使用一阶迎风格式[1]，即将（6）（7）组合：

（8）

2.3  高精度格式

扩散项内点处采用四阶的Padé格式离散[9]，格式如下：

（9）

对流项内点处采用三阶迎风紧致格式离散[10]，格式如下：

（10）

2.4  求解方法

首先初始化数组[u0，u1，u2，…，uN，uN+1]T。在边界条件已知的情况下，该问题可转化为通过矩阵运算求解向量U = [u1，u2，…，uN]T的问题。

下面本文将上文三种格式的计算方法以矩阵运算的方式给出，以便读者更加直觀的理解三种格式的计算方式。

2.4.1  二阶导数项的替换

1）将中心差分格式（迎风格式）化为矩阵格式： ，其中：

2）将高精度格式化为矩阵格式：，其中：

2.4.2  一阶导数项的替换

不失一般性，我们讨论当  时：

1）将中心差分格式化为矩阵格式为

2）将迎风格式化为矩阵格式为

3）将高精度格式化为矩阵格式为

2.5  算法实现

2.5.1  数值格式函数

文中根据以上三种格式，分别编写了五种数值格式函数（分别是中心差分格式，a（x）＞0时的迎风格式，a（x）＜0时的迎风格式，a（x）＞0时的高精度格式和a（x）＜0时的高精度格式），方便随时根据所需情况进行调用。下面给出高精度格式当a（x）＞0时的数值求解格式函数示例：

function [U1，U2] = zgjd（h，n）

%a（x）>0时的数值求解格式

%% 二阶差分

A=zeros（n-1，n-1）;%正文中的B2矩阵

B=zeros（n-1，n-1）;%正文中的B3矩阵

A（1，1）=5/6;A（1，2）=1/12;

A（n-1，n-2）=1/12;A（n-1，n-1）=5/6;

for i=1：n-3

A（i+1，i）=1/12;

end

for i=1：n-3

A（i+1，i+1）=5/6;

end

for i=1：n-3

A（i+1，i+2）=1/12;

end

B（1，1）=-2;B（1，2）=1;

B（n-1，n-2）=1;B（n-1，n-1）=-2;

for i=1：n-3

B（i+1，i）=1;

end

for i=1：n-3

B（i+1，i+1）=-2;

end

for i=1：n-3

B（i+1，i+2）=1;

end

U2=inv（A）*（（1/（h^2））.*B）;

%% 一阶差分

C=zeros（n-1，n-1）;%正文中的A3矩阵

D=zeros（n-1，n-1）;%正文中的A4矩阵

C（1，1）=2;C（1，2）=0;

C（n-1，n-2）=1;C（n-1，n-1）=2;

for i=1：n-3

C（i+1，i）=1;

end

for i=1：n-3

C（i+1，i+1）=2;

end

for i=1：n-3

C（i+1，i+2）=0;

end

D（1，1）=4/2;D（1，2）=1/2;

D（n-1，n-2）=-5/2;D（n-1，n-1）=4/2;

for i=1：n-3

D（i+1，i）=-5/2;

end

for i=1：n-3

D（i+1，i+1）=4/2;

end

for i=1：n-3

D（i+1，i+2）=1/2;

end

U1=（1/h）.*inv（C）*D;

end

2.5.2  数值实验与可视化

文中给出高精度格式求解算例1的数值实验代码示例：

%% 初始赋值

clear;clc

n=100;%设置步长

h=1/n;

R=zeros（n-1，1）;

eps=input（'输入epsilon的值： '）;

E=eye（n-1）;%算例1中u的系数矩阵

R=zeros（n-1，1）;%算例1中的f（x）

x=h;

for i=1：n-1

R（i，1）=x+1-（x-1）*exp（（x-1）/eps）;

x=x+h;

end

[U1，U2]=zgjd（h，n）;% 调用函数

%% 数值计算

T=-eps*U2+U1+E;

u=inv（T）*R;

U（1）=0; U（n+1）=0;% 边界值

for i=2：n;

U（i）=u（i-1）;% 内点值

end

%% 图像绘制

x=0：h：1;

M=x.*（1-exp（（x-1）/eps））;% 精确解

plot（x，U，'*：'，x，M，'k'）;

xlabel（'X'）;ylabel（'Solution'）;

legend（'本文格式'，'精确解'）;

set（gca，'fontweight'，'bold'，'linewidth'，2）

3  数值实验

为了验证数值格式的准确性、稳定性和有效性，我们考虑下面两个数值算例展开数值实验。本文所有的结果都是用双精度表示；所有程序与实验均在“MATLAB 2016a”软件上运行。

注1：最大值误差为 ；均方根误差为 。

注2：文章表中 。

3.1  算例1

考虑如下常系数对流扩散方程：

已知该方程的精确解为：

其中：

由图1，我们可以看出本文格式在ε较大的情况下，中心差分格式和高精度格式得到的数值解和精确解相吻合，但可以明显看出迎风格式的数值解与精确解存在偏离；但在ε较小的情况下，迎风格式和高精度格式的数值解和精确解吻合度较高，而中心差分格式的数值结果偏差则要大很多。即本文高精度格式相对于中心差分格式减少了数值振荡提高了ε较小时的数值计算精度，相对于迎风格式提高了ε较大时的数值计算精度。

分析表1结果，可知文中高精度格式相对于迎风格式在10-3＜ε＜10-1范围内精度提高了10倍，但在10-6＜ε＜10-4范围内精度有所降低，这是因为高精度格式和中心差分格式的数值解在边界x = 1附近出现数值振荡，但相对于中心差分格式精度提高了近102倍，可见本文高精度格式有效减少了数值震荡，这与图1的数值结果是相匹配的。

3.2  算例2

考虑如下变系数对流扩散方程：

已知该方程的精确解为：

其中：

由图2，我们可以看出本文格式在ε较大的的情况下，中心差分格式和高精度格式得到的数值解和精确解相吻合，但可以明显看出迎风格式的数值解与精确解存在偏离；但在ε较小的情况下，迎风格式和高精度格式的数值解和精确解吻合度较高，而中心差分格式的数值结果偏差则要大很多。

分析表2结果，可再次验证前文分析，文中高精度格式相对于迎风格式在10-3＜ε＜10-1范围内精度提高了10倍，但在10-6＜ε＜10-4范围内精度有所降低，与算例1不同的是，算例2种高精度格式和中心差分格式的数值解在边界x = 0附近出现了数值振荡，但相对于中心差分格式精度也提高了近102倍，减小了数值震荡。

算例2通过对变系数对流扩散方程的研究，进一步验证了由算例1所得出的结论。

4  结  论

本文讨论了三种奇异摄动对流扩散方程的求解格式，并使用MATLAB软件编写了求解相关问题的程序，对兩个算例进行了数值实验，主要从先验误差估计的角度入手，得出本文高精度格式减少了中心差分格式数值振荡提高了ε较小时的数值计算精度，提高了迎风格式ε较大时的数值计算精度的结论。通过与传统的中心差分格式和迎风格式的数值计算结果的对比，进一步验证了高精度格式的有效性和稳定性，为计算数学的发展提供了新的数据。

参考文献：

[1] 龙亮，张振华，姜林，等.放射性气体扩散方程有限差分法的MATLAB实现 [J].咸阳师范学院学报，2017，32（6）：40-42.

[2] 王忆锋，唐利斌.通过有限差分和MATLAB矩阵运算直接求解一维薛定谔方程 [J].红外，2010，31（3）：42-46.

[3] 冯立伟，徐涛，屈福志.泊松方程的有限差分法的MATLAB实现 [J].电脑知识与技术，2017，13（13）：233-235.

[4] 王涛，刘铁钢.求解对流扩散方程的一致四阶紧致格式 [J].计算数学，2016，38（4）：391-404.

[5] 杨录峰.几类奇异摄动问题的高精度数值方法研究 [D].兰州：兰州大学，2021.

[6] 岑仲迪.奇异摄动对流扩散问题的有限差分法 [D].杭州：浙江大学，2005.

[7] 王小妹，陈豫眉.求解一维对流方程的四阶紧致差分格式 [J].安阳师范学院学报，2022（2）：4-11.

[8] 何学飞.对流扩散方程的新型差分格式 [D].重庆：重庆大学，2016.

[9] Lele S K. Compact finite difference schemes with spectral-like resolutio [J].Journal of Computational Physics，1992，103（1）：16-42.

[10] 刘宏，傅德薰，马延文.迎风紧致格式与驱动方腔流动问题的直接数值模拟 [J].中国科学（A辑 数学 物理学 天文学 技术科学），1993（6）：657-665.

[11] 陆金甫，关治，偏微分方程数值解[M].北京：清华大学出版社，2004.

作者简介：潘轶航（2002—），男，汉族，福建南平人，本科在读，研究方向：偏微分方程数值解。