離散余弦擬譜方法求解反應(yīng)擴(kuò)散方程

張榮培, 楊程程, 劉 佳

(1. 沈陽(yáng)師范大學(xué) 數(shù)學(xué)與系統(tǒng)科學(xué)學(xué)院, 沈陽(yáng) 110034; 2. 沈陽(yáng)師范大學(xué) 大學(xué)外語(yǔ)教學(xué)部, 沈陽(yáng) 110034)

0 引 言

離散余弦變換(discrete cosine transform,DCT)主要用到余弦函數(shù)逼近給定樣本函數(shù)的最小二乘和插值[1],在信號(hào)處理[2]、圖像壓縮[3]以及偏微分方程求解領(lǐng)域[4-5]得到了廣泛的應(yīng)用。由于DCT能夠?qū)⒖沼虻男盘?hào)轉(zhuǎn)換到頻域上,具有良好的相關(guān)性的性能。

其次,將離散余弦變換應(yīng)用于求解偏微分方程。將要逼近的函數(shù)u的最小二乘余弦插值多項(xiàng)式代入微分方程,余弦擬譜離散化應(yīng)用于空間導(dǎo)數(shù),然后在這些離散點(diǎn)集上應(yīng)用配置方法,結(jié)合離散余弦變換,最終得到關(guān)于頻域上的常微分方程組。通過(guò)向后Euler方法離散時(shí)間導(dǎo)數(shù)使之具有更好的穩(wěn)定性,應(yīng)用Picard迭代求解非線性代數(shù)方程組。在求解過(guò)程中,使用了快速余弦變換,提高了計(jì)算效率。

1 離散余弦變換

首先,將求解區(qū)域劃分成J個(gè)網(wǎng)格,然后取每個(gè)網(wǎng)格的中心點(diǎn)得到離散點(diǎn)集。對(duì)于函數(shù)u,將[a,b]剖分為J個(gè)網(wǎng)格,定義

定理1 當(dāng)-(2J-1)≤k≤2J-1時(shí),有

(1)

故得證。

定理2 函數(shù){φk},k=0,1,…J-1在點(diǎn)集{xj}上正交,即

(2)

當(dāng)k=m時(shí),式(2)顯然成立。當(dāng)k≠m時(shí),即k-m≠0,有-(J-1)≤k-m≤J-1,1≤k+m≤2J-1由定理1可證。

這說(shuō)明了余弦函數(shù){φk},k=0,1,…,J-1在點(diǎn)集{xj}(j=1,2,…,J)上是正交的。

(3)

(4)

如果將離散點(diǎn)集代入最小二乘余弦逼近多項(xiàng)式(3)中,得到函數(shù)值uj等于SN(xj)。下面給出證明,

(5)

2 在反應(yīng)擴(kuò)散方程中的應(yīng)用

反應(yīng)擴(kuò)散方程(組)是一類非線性偏微分方程,在生態(tài)學(xué)[6-8]、燃燒學(xué)[9]和化學(xué)反應(yīng)流[10]等領(lǐng)域有著重要的應(yīng)用。近年來(lái),許多學(xué)者對(duì)該方程提出了許多數(shù)值方法,比如緊致積分因子法[11]、間斷有限元法[12]、 格子Boltzmann方法[13]。下面將上面的離散余弦變換應(yīng)用于求解偏微分方程,得到余弦擬譜方法[14]?？紤]如下一維反應(yīng)擴(kuò)散方程

(6)

由定理3和上述二階導(dǎo)數(shù)離散得

該非線性方程組可由如下Picard迭代法求解

3 數(shù)值算例

下面給出反應(yīng)擴(kuò)散方程的一些算例。在凝膠反應(yīng)中,Gray-Scott方程模擬化學(xué)反應(yīng)。對(duì)應(yīng)的無(wú)量綱反應(yīng)擴(kuò)散方程為

其中:u,v為系統(tǒng)變量,代表反應(yīng)物——基質(zhì)和催化劑的濃度;Du,Dv為擴(kuò)散系數(shù);A,B為系統(tǒng)控制參量——無(wú)尺度的流速度和催化劑的衰變常數(shù)[15]。

圖1 u,v的數(shù)值解圖,當(dāng)A=0.01,B=0.1時(shí)Fig.1 The plot of solution of u,v at A=0.01,B=0.1

在圖2中,取t=4 000,A=0.01,B=0.1,可以看到,2個(gè)脈沖每一個(gè)都分裂成2個(gè)行波脈沖,現(xiàn)在產(chǎn)生了4個(gè)脈沖,又過(guò)了一段時(shí)間,這4個(gè)脈沖再一次分裂,此時(shí)得到了8個(gè)脈沖。這個(gè)過(guò)程在2個(gè)脈沖之外繼續(xù),直到填滿整個(gè)間隔,在這個(gè)自我復(fù)制的過(guò)程中,v在2個(gè)峰之間成幅度減小,在脈沖之間u不趨近于1,u的最大值隨著2個(gè)行波脈沖之間距離的增加而增加。

圖2 u,v的數(shù)值解圖,當(dāng)A=0.01,B=0.1時(shí)Fig.2 The plot of solution of u,v at A=0.01,B=0.1

圖3 u,v的數(shù)值解圖,當(dāng)Fig.3 The plot of solution of u,v at

4 結(jié) 論

本文推導(dǎo)了離散余弦變換是如何得到的,通過(guò)離散余弦變換的方法對(duì)反應(yīng)擴(kuò)散方程進(jìn)行了求解。首先通過(guò)考慮區(qū)域上有具有齊次Neumann邊界的函數(shù)u,將求解區(qū)域劃分為N個(gè)網(wǎng)格,取每個(gè)網(wǎng)格的中心點(diǎn)得到離散點(diǎn)集。然后根據(jù)余弦函數(shù)族在點(diǎn)集上的正交性,得到函數(shù)u的最小二乘余弦逼近多項(xiàng)式SN。將離散點(diǎn)集代入SN中,可以證明SN(xj)=uj。通過(guò)建立映射得到了離散余弦變換,然后應(yīng)用離散余弦變換得到了關(guān)于頻域上的常微分方程組,通過(guò)向后Euler方法離散時(shí)間導(dǎo)數(shù),最后應(yīng)用Picard迭代求解非線性代數(shù)方程組。通過(guò)對(duì)數(shù)值算例的數(shù)值結(jié)果進(jìn)行仿真計(jì)算,證實(shí)了通過(guò)離散余弦變換來(lái)求解反應(yīng)擴(kuò)散方程是有效的,快速的。