基于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—），男，漢族，福建南平人，本科在讀，研究方向：偏微分方程數值解。