李友愛

(北京工商大學(xué) 計(jì)算機(jī)與信息工程學(xué)院, 北京 100048)

1 問題的提出

由于具有抗腐蝕能力強(qiáng)、不與酸堿反應(yīng)、制造成本低、防水、質(zhì)輕、容易被塑制成不同形狀等眾多優(yōu)點(diǎn)，塑料已經(jīng)被廣泛地應(yīng)用于食品和藥品包裝. 但是這些塑料包裝材料中所含有的化學(xué)物可能遷移到被包裝的食品和藥品中，因而會(huì)危害人體健康. 如何檢測(cè)和分析塑料包裝材料化學(xué)物的遷移并形成相關(guān)的理論，已經(jīng)成為食品和藥品包裝領(lǐng)域一個(gè)重要的研究課題. 目前，主要檢測(cè)和分析手段是實(shí)驗(yàn)的方法. 然而，化學(xué)物遷移實(shí)驗(yàn)對(duì)實(shí)驗(yàn)儀器要求非常高，且實(shí)驗(yàn)復(fù)雜、費(fèi)用昂貴，很難在普通的實(shí)驗(yàn)室完成. 因此，很多國內(nèi)外的研究者將研究的重點(diǎn)向理論研究轉(zhuǎn)移，其中的一個(gè)研究熱點(diǎn)就是將化學(xué)物的遷移轉(zhuǎn)化為一個(gè)微分方程，然后通過求解這一微分方程模擬化學(xué)物的遷移.

大量科學(xué)實(shí)驗(yàn)表明(Figge[1]，Chang[2]， Till等[3])， 當(dāng)化學(xué)物在塑料內(nèi)的擴(kuò)散系數(shù)D以及它在塑料薄膜—食品間的分配系數(shù)K已知時(shí)，可基于Fick第二定律對(duì)遷移過程進(jìn)行預(yù)測(cè). 通常認(rèn)為遷移僅發(fā)生在包裝材料厚度方向上，可用一維二階偏微分方程式

(1)

描述. 若擴(kuò)散系數(shù)D與濃度無關(guān)，則上式簡(jiǎn)化為

(2)

對(duì)于上述遷移模型，只有在如下特殊的情況下才有可能求得解析解[4]:1)初始時(shí)刻，化學(xué)物均勻分布在包裝薄膜中；2)化學(xué)物從包裝薄膜一側(cè)進(jìn)入食品，交界面處沒有傳質(zhì)阻力(傳質(zhì)系數(shù) 很大)；3)任一時(shí)刻食品中的化學(xué)物均勻分布；4)在整個(gè)遷移過程中，擴(kuò)散系數(shù)D和分配系數(shù)KP/F=CP,e/CF,e為常數(shù)；5)遷移過程任何時(shí)刻，在包裝薄膜和食品的界面上都是平衡的； 6)忽略包裝材料邊界效應(yīng)及其與食品的相互作用. 一般情況下，我們無法獲得解析方法而只能借助于數(shù) 值方法.

對(duì)于模型(2),常用的空間離散方法有有限差分[5]、有限元[6]和有限體積方法. 對(duì)空間離散后，在時(shí)間方向的離散常用的有顯式格式和隱式格式. 用顯式格式計(jì)算簡(jiǎn)單，但是為了保證數(shù)值穩(wěn)定性，時(shí)間步長(zhǎng)必須較小；而用隱式格式計(jì)算時(shí)，可取較長(zhǎng)的時(shí)間步長(zhǎng)，但計(jì)算量要大得多. 并且這些方法有個(gè)共同的缺點(diǎn)，即時(shí)間步長(zhǎng)相對(duì)較小，且精度不是很高. 鐘萬勰等在文獻(xiàn)[7]中提出了一種精細(xì)時(shí)程積分法. 這種方法的一個(gè)顯著特點(diǎn)就是具有高精度. 本文的目的是將精細(xì)時(shí)程積分法推廣應(yīng)用到化學(xué)物遷移模型(2)并給出誤差分.

2 精細(xì)時(shí)程積分法

一般而言，化學(xué)物遷移模型(2)經(jīng)過空間半離散后將得到如下關(guān)于時(shí)間的常微分方程組：

(3)

其中，v為濃度態(tài)變量,H為常系數(shù)矩陣,r為非齊次項(xiàng)，表示系統(tǒng)外部作用量. 方程(3)的解v因非齊次項(xiàng)r(t) 而分為兩部分：v1和v2，即

v(t)=v1(t)+v2(t),

(4)

v1(t)=exp(Ht)·v0,

(5)

(6)

其中v0為初始條件.則該系統(tǒng)的解在時(shí)刻tk為

(7)

根據(jù)方程(7), 解在時(shí)刻tk+1(τ=tk+1-tk, 其中τ為時(shí)間步長(zhǎng)) 為

(8)

通過推導(dǎo)vk和vk+1的關(guān)系，可得逐步積分公式

(9)

式(9)為精細(xì)積分法的基本公式，是基于解析方法導(dǎo)出的，將初始條件代入上式，即可逐步求出系統(tǒng)在各時(shí)刻的解，精細(xì)積分法的關(guān)鍵是計(jì)算上式中的指數(shù)矩陣

(10)

為簡(jiǎn)單起見，應(yīng)用二階Taylor展開計(jì)算上式中括號(hào)內(nèi)的部分，得

(11)

其中

在式(11)中，高階的Taylor展開值得推薦，會(huì)提高結(jié)果的精度. 注意到 Δt很小，故Tα相對(duì)于I也很小. 如果直接將他們累加會(huì)導(dǎo)致有效數(shù)位的很大損失. 有鑒于此，精確積分算法先將值較小的矩陣Tα累加，然后將之加到矩陣I. 注意到

(I+Tα)2=I+2Tα+Tα×Tα,

(12)

算法設(shè)計(jì)為

for (iter=0;iter

(13)

當(dāng)循環(huán)結(jié)束后，

T=I+Tα.

(14)

方程(13)和(14)是PTI計(jì)算指數(shù)矩陣(10)的關(guān)鍵步驟.

3 精細(xì)時(shí)程積分法的誤差估計(jì)

若r(t)=0，則方程(3)簡(jiǎn)化為

(15)

H為常系數(shù)矩陣,其通解可形式地寫成

v=exp(H·τ)v0,

(16)

其中v0為初始時(shí)刻的值. 令T=exp(H·τ)=[exp(H·τ/m)]m，又令Δt=τ/m，m=2N(在文獻(xiàn)[7]中，取N=20)，則

T=[exp(H·Δt)]m.

令A(yù)=exp(H·Δt)≈I+HΔt+(HΔt)2/2!. 下面考察用精細(xì)時(shí)程積分法來近似精確解的誤差. 不妨假定H對(duì)稱負(fù)定，λ1,…,λn為H的特征值，且λi≤0. 又設(shè)|λ1|≤|λ2|≤…≤|λn|， 則

H=Qdiag(λ1,…,λn)QT,

(17)

且QQT=I. 令vh=Amv0，于是誤差為

v-vh=(exp(H·τ)-Am)v0=
Qdiag(Δ1,Δ2,…,Δn)QTv0,

(18)

exp(λ1Δt)=1+λ1Δt+(λ1Δt)2/2!+
exp(λ1Δξ1)(λ1Δt)3/3!,
0 <Δξ1<Δt,

(19)

即

1+λ1Δt+(λ1Δt)2/2! =exp(λ1Δt)-
exp(λ1Δξ1)(λ1Δt)3/3!.

(20)

由方程(20)及(18)可知

Δ1=exp(λ1τ)-(1+λ1Δt+(λ1Δt)2/2!)m=
exp(λ1τ) -(exp(λ1Δt)-
exp(λ1Δξ1)(λ1Δt)3/3!)m=
exp(λ1τ) -exp(λ1Δtm)(1-
exp(λ1(Δξ1-Δt))(λ1Δt)3/3!)m=
exp(λ1τ)-exp(λ1τ)(1+x)m=
exp(λ1τ)[1-(1+x)m],

(21)

其中x=-exp(λ1(Δξ1-Δt))(λ1Δt)3/3!. 利用一階Taylor展開可得

(1+x)m=1+m(1+y1)m-1x, 0

(22)

從式(22)可知和(21)式, 可得

Δ1=-exp(λ1τ)m(1+y1)m-1x=
-exp(λ1τ)m(1+y1)m-1(-exp(λ1(Δξ1-
Δt))(λ1Δt)3/3!)=
exp(λ1τ)exp(λ1(Δξ1-Δt))m(1+
y1)m-1(λ1Δt)3/3!=
exp(λ1(Δξ1-Δt))mexp(λ1τ)(1+
y1)m-1(λ1Δt)3/3!.

(23)

下面依次考察式(23)中各項(xiàng). 由于 0

exp(λ1τ)(1+y1)m-1=
exp(λ1Δt)(exp(λ1Δt)+exp(λ1Δt)y1)m-1≤
exp(λ1Δt)(exp(λ1Δt)-
exp(λ1Δξ1)(λ1Δt)3/3!)m-1=
exp(λ1Δt)(1+λ1Δt+(λ1Δt)2/2!)m-1.

(24)

因?yàn)棣=τ/m,λ1<0， 假設(shè)|λ1|≤|λ2|≤…≤|λn|

由上知,(24)式為

(25)

將(25)式代入(23)式得

Δ1≤exp(λ1Δξ1)|(λ1Δt)3/3!|.

(26)

定理1 設(shè)v為方程(15)精確解，vh為用精細(xì)時(shí)程積分法得到的近似解， 則v-vh的l2誤差為

‖v-vh‖l2≤max1≤i≤ncexp(λiΔξi)|λi/m|3τ3‖v0‖l2,

(27)

其中λi,i=1,2,…,n為矩陣A的特征值，m=2N，τ為時(shí)間步長(zhǎng)，v0為初始時(shí)刻的值.

證明：在式(18)的兩邊左乘QT，令u=QT(v-vh)，并注意到Q為正交矩陣，則

u=QT(v-vh)=diag(Δ1,Δ2, …,Δn)QTv0.

(28)

令

y=QTv0=(y1,…,yn)T,

(29)

則

u=diag(Δ1,Δ2, …,Δn)y.

(30)

于是u的離散l2范數(shù)為

(31)

因Q為正交矩陣，故有

(32)

‖v-vh‖l2=‖u‖l2,

(33)

即

‖v-vh‖l2≤max1≤i≤n(Δi)‖v0‖l2.

(34)

于是有

‖v-vh‖l2≤
max1≤i≤ncexp(λiΔξi)|λi/m|3τ3‖v0‖l2,

(35)

其中i=1,…,n，c表示不同的常數(shù).