擴散方程高精度加權(quán)差分格式的MATLAB實現(xiàn)

胡 敏

（攀枝花學院 數(shù)學與計算機學院，四川 攀枝花617000）

0 引言

MATLAB是一款具有精確數(shù)值計算功能和豐富圖形處理函數(shù)的軟件，[1]偏微分方程的數(shù)值解可直觀地以二維、三維圖形方式顯示在屏幕上.因此，近年來，越來越多的人開始使用MATLAB來求解偏微分方程.[2-7]一維擴散方程是最簡單的偏微分方程之一，其定解問題的數(shù)值方法有差分法、有限元法和邊界元法.本文主要討論一維擴散方程的一個高精度加權(quán)差分格式的MATLAB實現(xiàn).[8]

1 求解擴散方程的基本思想

用有限差分法求解偏微分方程問題必須把連續(xù)問題進行離散化，因此求解擴散方程的基本思想是把連續(xù)的定解區(qū)域用有限個離散點構(gòu)成的網(wǎng)格來代替，這些離散點稱作網(wǎng)格的節(jié)點.[9]

把連續(xù)定解區(qū)域上連續(xù)變量的函數(shù)用定義在網(wǎng)格上的離散變量函數(shù)來近似，于是原微分方程和定解條件就近似地以代數(shù)方程組代之，即有限差分方程組，解此方程組就可以得到原問題在離散點上的近似解.然后再利用插值方法便可以從離散解得到定解問題在整個區(qū)域上的近似解.[10]

2 擴散方程的高精度加權(quán)差分格式

一維擴散方程

其中，擴散系數(shù)a為常數(shù).先進行網(wǎng)格剖分，取空間步長h和時間步長τ，用兩族平行直線x＝xj＝j(luò)h（j＝0，1，…，N）和t＝tk＝kτ（k＝0，1，…，M）（其中N，M 都是正整數(shù)）將矩形域分割成矩形網(wǎng)格，網(wǎng)格節(jié)點為 （xj，tk）.用ukj表示定義在網(wǎng)點 （xj，tk）的函數(shù)，0≤j ≤N，0≤k ≤M.

利用二階微商三次樣條四階逼近公式[11]有

其中 （uxx）j表示二階偏導(dǎo)數(shù)uxx在點 （jh，t）處的近似值.

由于方程（1）在整個求解區(qū)域內(nèi)成立，于是將（2）式代入（1）式得

由于（3）式在任意時刻均成立，則

將（4）式和（5）式分別乘以θ和1-θ（0≤θ≤1），再相加得

又因為

于是得高精度加權(quán)差分格式[8]

3 高精度加權(quán)差分格式的MATLAB編程

利用上述高精度差分格式求擴散方程的近似解，需要求解一個大型稀疏矩陣方程組，這時采用迭代法是最合適的，而且計算出的結(jié)果也比較精確，且其程序設(shè)計簡單，能有效地解決一些高階問題，是解大型稀疏方程組的一種重要方法.本文采用高斯-賽德爾（Gauss-Seidel）迭代法來實現(xiàn)擴散方程的高精度加權(quán)差分格式所構(gòu)造的稀疏矩陣方程組的求解.

%解一維擴散方程 格式 （Ut-aUxx＝0，a＞0）

%用g-s（高斯-賽德爾）迭代法解

%m，n為x，t方向的網(wǎng)格數(shù)

format long

syms temp1；

syms temp2；

syms temp3；

a＝1；

l＝pi；

T＝pi＾2；

N＝10；

M＝200；

h＝l／N；

to＝T／M；

r＝（a＊to）／h＾2；

for j＝1：N＋1

x（j）＝（j-1）＊h；

for k＝1：M＋1

t（k）＝（k-1）＊to；

u0（j，k）＝exp（-t（k））＊sin（x（j））；

end

end

u0 %求解精確解u0

for j＝1：N＋1

x（j）＝（j-1）＊h；

u1（j，1）＝sin（x（j））；

u2（j，1）＝sin（x（j））；

u3（j，1）＝sin（x（j））；

end

for k＝1：M＋1

t（k）＝（k-1）＊to；

u1（1，k）＝0；u1（N＋1，k）＝0；

u2（1，k）＝0；u2（N＋1，k）＝0；

u3（1，k）＝0；u3（N＋1，k）＝0；

end

for k＝1：M

for j＝2：N

temp1＝（（1＋9＊r）＊（u0（j-1，k）

＋u0（j＋1，k））＋（10-18＊

r）＊u0（j，k）＋...-（1-9＊

r）＊（u0（j-1，k＋1）＋u0（j

＋1，k＋1）））／（10＋18＊

r）；

u1（j，k＋1）＝temp1；

temp2＝（（1＋6＊r）＊（u0（j-1，k）＋

u0（j＋1，k））＋（10-12＊r）＊

u0（j，k）＋...-（1-6＊r）＊

（u0（j-1，k＋1）＋u0（j＋1，k

＋1）））／（10＋12＊r）；

u2（j，k＋1）＝temp2；

temp3＝（（1＋3＊r）＊（u0（j-1，k）＋

u0（j＋1，k））＋（10-6＊r）＊

u0（j，k）＋...-（1-3＊r）＊

（u0（j-1，k＋1）＋u0（j＋1，k

＋1）））／（10＋6＊r）；

u3（j，k＋1）＝temp3；

end

end

u1 %求解數(shù)值解u1

u2 %求解數(shù)值解u2

u3 %求解數(shù)值解u3

for k＝1：M＋1

for j＝1：N＋1

R1（j，k）＝abs（u0（j，k）-u1（j，k））；

R2（j，k）＝abs（u0（j，k）-u2（j，k））；

R3（j，k）＝abs（u0（j，k）-u3（j，k））；

end

end

R1，R2， R3 %計算誤差R1，R2， R3

Rmax1＝max（R1） %求誤差R1的最大值

Rmax2＝max（R2） %求誤差R2的最大值

Rmax3＝max（R3） %求誤差R3的最大值

%精確解與數(shù)值解的比較：

x＝0：0.1：1；

hold on

plot（x，u0（：，M＋1），＇b＇）；

plot（x，u1（：，M＋1），＇y＇）；

plot（x，u2（：，M＋1），＇r＇）；

plot（x，u3（：，M＋1），＇g＇）；

title（＇t＝1，h＝0.1，τ＝0.005時精確解和數(shù)值解的比較＇）

hold off

4 運行結(jié)果分析

考慮擴散方程混合問題

其解析解為u（x，t）＝e-tsinx.

取r＝0.5，θ＝0.25，θ＝0.5和θ＝0.75代入高精度加權(quán)差分格式（8）計算數(shù)值解，在某時刻的最大誤差分別記作e1，e2，e3，數(shù)據(jù)如表1.

在同一個坐標系中，繪出二維圖作比較，如圖1.

表1 某時刻的最大誤差

圖1 數(shù)值解與精確解的比較

從表1和圖1可以看出，當θ＝0.5時，誤差最小；并且當網(wǎng)格分得越小時誤差更小.建議在解決實際問題時選擇誤差小的格式進行分析求解.

[1]（美）Stephen J，Chapman.MATLAB編程[M].邢樹軍，鄭碧波譯.北京：科學出版社，2008：3.

[2]王 飛，裴永祥.有限差分方法的 MATLAB編程[J].新疆師范大學學報：自然科學版，2003（4）：21-27.

[3]趙德奎，劉 勇.MATLAB在有限差分數(shù)值計算中的應(yīng)用[J].四川理工學院學報，2005（4）：61-64.

[4]田 兵.用 MATLAB解偏微分方程[J].陰山學刊，2006（4）：12-13.

[5]謝煥田，吳 艷.拉普拉斯有限差分法的 MATLAB實現(xiàn)[J].四川理工學院學報，2008（3）：1-2.

[6]史 策.熱傳導(dǎo)方程有限差分法的 MATLAB實現(xiàn)[J].咸陽師范學院學報，2009（4）：27-29.

[7]薛 瓊，肖小峰.二維熱傳導(dǎo)方程有限容積法的 MATLAB實現(xiàn)[J].計算機工程與應(yīng)用，2009（4）：27-29.

[8]田振夫，張艷萍.擴散方程的高精度加權(quán)差分格式[J].中國科學技術(shù)大學學報，1999（2）：237-241.

[9]陳金甫，關(guān) 治.偏微分數(shù)值解法：第二版[M].北京：清華大學出版社，2004：274.

[10]李瑞遐，何志慶.微分方程數(shù)值方法[M].北京：華東理工大學出版社，2005：97.

[11]田振夫.泊松方程的高精度三次樣條差分方法[J].西北師范大學學報：自然科學版，1996（2）：13-17.