自由邊界問題的重心插值迭代配點(diǎn)法研究

趙岳月，王兆清，*，李金

（1.山東建筑大學(xué) 工程力學(xué)研究所，山東 濟(jì)南250101；2.山東建筑大學(xué) 理學(xué)院，山東 濟(jì)南250101）

0 引言

許多工程問題都可以轉(zhuǎn)化為自由邊界問題求解，如自由邊界隨時(shí)間變化時(shí)水的融化和凝固問題［1］（典型的 Stefan問題［2］）、氣體在多孔介質(zhì)中的吸收和擴(kuò)散（氧氣通過生物薄膜的吸收和擴(kuò)散問題［3］）；以及自由邊界值不隨時(shí)間變化的巖土體的無壓滲流面問題［4-6］等。自由邊界問題一般采用微分方程描述，考慮一類控制方程為二階線性微分方程的自由邊界問題（未知邊界為固定值），其形式由式（1）表示為

式中：p（x）、q（x）、f（x）在［a，bf］上連續(xù)可導(dǎo)；bf和y（x）未知。求解自由邊界問題的難點(diǎn)在于微分方程的求解區(qū)域的部分邊界是待定的，其待定邊界bf要作為解的一部分來求。

為使自由邊界問題有解，在自由邊界上一般給定過約束邊界條件。其中一個(gè)邊界條件用于求解微分方程，另一個(gè)邊界條件用于確定自由邊界位置。自由邊界位置的計(jì)算精度依賴于所選擇的數(shù)值計(jì)算方法。有限差分法［7-9］和有限元方法［10-11］是2種經(jīng)典的求解微分方程自由邊值問題的數(shù)值方法，其計(jì)算精度依賴于差分步長 ／單元尺寸的大小，要得到問題的高精度解，需要細(xì)化計(jì)算網(wǎng)格。自適應(yīng)投影方法將自由邊界問題用有限差分格式離散為一個(gè)線性互補(bǔ)問題，在投影迭代過程中自動(dòng)調(diào)整參數(shù)［12］；改進(jìn)的插值型邊界無單元法不僅要循環(huán)積分，還需要計(jì)算邊界節(jié)點(diǎn)上的插值型最小二乘法的形函數(shù)矩陣［13］，其計(jì)算相對(duì)比較復(fù)雜且計(jì)算節(jié)點(diǎn)較多。譜Tau方法是一種高精度數(shù)值計(jì)算自由邊界問題的數(shù)值方法［14］，選取較少的節(jié)點(diǎn)計(jì)算得到了高精度的數(shù)值解。重心插值配點(diǎn)法則采用重心型插值函數(shù)作為未知函數(shù)的近似函數(shù)，利用重心插值微分矩陣可以直接得到微分方程的配點(diǎn)法離散計(jì)算矩陣［15］。重心插值配點(diǎn)法具有類似譜方法的計(jì)算精度。

文章給定一個(gè)自由邊界初始假設(shè)值，采用重心插值配點(diǎn)法求解微分方程，利用自由邊界上的任意一個(gè)定界條件構(gòu)造出確定自由邊界位置的2種迭代格式，即Newton法和弦截法，提出了數(shù)值求解自由邊值問題的重心插值迭代配點(diǎn)法，并以數(shù)值算例進(jìn)行分析，驗(yàn)證了迭代配點(diǎn)法對(duì)于自由邊界問題求解的可行性和精確性。

1 自由邊值問題的重心插值迭代配點(diǎn)法建立

1.1 邊值問題的重心插值配點(diǎn)法

將微分方程（1）的未知函數(shù)近似為重心插值型多項(xiàng)式 xi，為區(qū)間［a，b］上按第二類切比雪夫（Chebyshev）節(jié)點(diǎn)分布離散的N個(gè)節(jié)點(diǎn)，在此節(jié)點(diǎn)上的重心插值由式（2）［15］表示為

式中：wj為插值權(quán)其中j＝1，2，…，N。

近似插值函數(shù)利用微分矩陣離散的形式由式（3）表示為

式中：P ＝ diag［p（xi）］；Q ＝ diag［q（xi）］；y ＝［y（xi）］T；f＝［f（xi）］T；P、Q為對(duì)角矩陣；y、f為列向量；I為n階單位矩陣；D（m）為未知函數(shù)的m階插值微分矩陣［10］。

將邊界條件利用微分矩陣離散后施加到式（3）中，求解近似函數(shù)y（x）。邊界條件的施加方法有消去法、置換法、附加法。為了方便處理邊界條件，選用了置換法。

1.2 自由邊值問題的重心插值迭代配點(diǎn)法

求邊值問題需要先假設(shè)一個(gè)bf的值，以式（1）的自由邊界問題為例，y（bf）和 y′（bf）任意與 y（a）構(gòu)成邊值問題的邊界條件，bf由第三個(gè)條件確定。假設(shè)值一般不是真值，于是采用迭代法逼近準(zhǔn)確值，文章引進(jìn)2種迭代格式Newton法和弦截法以滿足計(jì)算方法的普遍適用性。

1.2.1 Newton迭代格式

（1）當(dāng)施加的邊界條件為第一類邊界條件，判定條件為自由邊界的導(dǎo)數(shù)值時(shí)，將近似函數(shù)的導(dǎo)函數(shù)在bf處（設(shè)bf是準(zhǔn)確值b*f的近似值）按泰勒公式展開，由式（4）表示為

當(dāng) y″（bf）≠0時(shí)，有于是得到了由y′判定的Newton迭代公式，由式（5）表示為

式中為第k+1次迭代得到的自由邊界值；為k次得到的自由邊界值；y′（）為利用的值插值得到的近似函數(shù)的導(dǎo)函數(shù)在處的值；y″（）為近似函數(shù)在bkf處的二階導(dǎo)數(shù)值。

（2）當(dāng)施加的邊界條件為混合邊界條件，判定條件為自由邊界的函數(shù)值時(shí)，類似于（1）的推導(dǎo)，將近似函數(shù)在bf處泰勒展開為bf）+…，得到由y判定的Newton迭代格式，由式（6）表示為

由于計(jì)算近似函數(shù)時(shí)已經(jīng)將邊界條件的一階導(dǎo)數(shù)值施加到了微分方程中，那么y′（）為已知的定值，可能會(huì)造成計(jì)算中迭代不收斂，且當(dāng)y′（）＝0時(shí)，迭代格式無意義。于是引進(jìn)了另一種迭代格式弦截法。

1.2.2 弦截法迭代格式

將式（6）中近似函數(shù) y（x）在點(diǎn)處的切線斜率y′（）改用曲線上2點(diǎn)弦的斜率代替，即得到由y判定的弦截法迭代格式，由式（7）表示為

式中分別為第 k+2、k+1、k次迭代得到的自由邊界值是利用的值插值得到的近似函數(shù)在處的值。

弦截法的迭代格式中要得到的值需要有第k+1次的的值第k次的的值，所以在計(jì)算過程中要假設(shè)2個(gè)初始值。在計(jì)算中直接令，則那么在此迭代格式中只需要假設(shè)一個(gè)初始值。

同樣地，將式（4）中y″（）也改用曲線上2點(diǎn)弦的斜率代替，得到y(tǒng)′判定的弦截法迭代格式，由式（8）表示為

1.2.3 重心插值配點(diǎn)法的迭代步驟

（1）令 k＝1，且 k≤100，假設(shè)的 b′f值，（弦截法令＝a，假設(shè)的值）；

（2）將求解區(qū)域［a］（弦截法離散成N個(gè)節(jié)點(diǎn)；用式（2）、（3）計(jì)算第k次近似函數(shù)yk，得到的值，其中。（在弦截法中

（3）用迭代公式（7）計(jì)算第k+1次的自由邊界值（弦截法用式（5）計(jì)算第k+2次的自由邊界

值）和的值；

（4）Newton法用 y′（）作為判定條件，若則重復(fù)過程（2）和（3）；若k＞100，解不收斂，跳出循環(huán)；若 y′（）＜ε，迭代結(jié)束，輸出、yk。為所求自由邊界的近似值；yk為在邊界［a］上的近似函數(shù)。

（5）弦截法用作為判定條件，若，重復(fù)過程（2）和（3）；若k＞100，解不收斂，跳出循環(huán)；若，迭代結(jié)束，輸出為所求自由邊界的近似值，yk+1為在邊界上的近似函數(shù)。

2 數(shù)值算例驗(yàn)證

數(shù)值算例的計(jì)算程序利用MATLAB編寫，其中m為自由邊值的迭代次數(shù)，N為插值節(jié)點(diǎn)的數(shù)量。

2.1 算例1

二階線性微分方程的自由邊值問題［14］。

邊界條件為

問題的解析解為

在問題求解中將y（0）＝1和y（xs）＝0作為邊界條件，y′（xs）＝0作為判定條件，采用2種迭代格式來計(jì)算，比較其迭代次數(shù)和誤差精度，誤差的控制精度為10-10。假設(shè)初始邊界為xs＝1，其結(jié)果見表1和2。

表1 Newton法迭代配點(diǎn)表

表2 弦截法迭代配點(diǎn)表

由表1和2可知，2種方法求解微分方程的自由邊界問題時(shí)，均得到了10-11～10-13量級(jí)的高精度解。在邊界條件和節(jié)點(diǎn)數(shù)量相同的情況下，Newton迭代法要比弦截法的收斂速度快，只需要迭代3、4次就可以得到滿足控制誤差精度的解。而對(duì)于xs和近似函數(shù)y（x）的誤差弦截法要稍稍優(yōu)于Newton法約10-2，但是在實(shí)際計(jì)算中Newton法的計(jì)算精度足以滿足微分方程的計(jì)算要求。

2.2 算例2

線性微分方程自由邊值問題。

邊界條件為

問題的解析解為

在問題求解中首先將y（0）＝1和y′（xs）＝0作為邊界條件，y（xs）＝0作為判定條件，假設(shè)初始邊界為xs＝1。由于y′（xs）＝0，采用Newton迭代格式時(shí)無法計(jì)算，于是選用弦截法迭代格式。誤差的控制精度為10-10，計(jì)算結(jié)果見表3。而當(dāng)選取y（0）＝1和y（xs）＝0作為邊界條件，y′（xs）＝0作為判別條件時(shí)，Newton迭代法和弦截法都可以求解。

表3 弦截法迭代配點(diǎn)表

3 結(jié)論

通過上述算例分析表明，2種迭代格式都可以解線性微分方程的自由邊值問題，而且計(jì)算精度都很高。

（1）Newton法的迭代速度較弦截法快，迭代3、4次就可以得到高精度的解，且計(jì)算公式簡單。

（2）由于邊界條件以及控制方程自身的原因，可能會(huì)造成Newton迭代配點(diǎn)法無法計(jì)算，而對(duì)弦截法計(jì)算無影響，其計(jì)算精度與Newton法基本相當(dāng)，只是在初始值的處理上比較復(fù)雜。

（3）重心插值迭代配點(diǎn)法是一種高精度的迭代配點(diǎn)法，在求解線性微分方程自由邊界問題時(shí)，可以根據(jù)問題的特點(diǎn)選取適當(dāng)?shù)牡袷?。?jì)算誤差精度隨節(jié)點(diǎn)的增加成量級(jí)提高，可以達(dá)到10-11～10-13。在計(jì)算自由邊值問題時(shí)具有計(jì)算公式簡單，數(shù)值程序?qū)嵤┓奖?，?jì)算精度高等優(yōu)點(diǎn)。

參考文獻(xiàn)：

［1］Mitchell S L，Vynnycky M，Gusev I G，et al.An accurate numerical solution for the transient heating of an evaporating spherical droplet［J］.Applied Mathematics＆ Computation，2011，217（22）：9219-9233.

［2］Borodin M A.Stefan problem［J］.Journal of Mathematical Sciences，2011，178（1）：13-40.

［3］Bougoffa L.The adomian decomposition method for solving a moving boundary problem arising from the diffusion of oxygen in absorbing tissue［J］.The Scientific World Journal，2014，2014：1-7.

［4］王兆清，李術(shù)才，李樹忱.無壓滲流問題分析的多節(jié)點(diǎn)有限元方法［J］.巖土力學(xué)，2008，29（10）：2647-2650.

［5］李樹忱，王兆清，袁超.巖土體滲流自由面問題的重心插值無網(wǎng)格方法［J］.巖土力學(xué)，2013（7）：1867-1873.

［6］鄭珊珊.求解滲流自由面的逐步剖分法［D］.煙臺(tái)：煙臺(tái)大學(xué)，2017.

［7］Kim B J，Ahn C，Choe H J.Direct computation for American put option and free boundary using finite differencemethod［J］.Japan Journal of Industrial＆Applied Mathematics，2013，30（1）：21-37.

［8］Devi P K，Naidu V G.A new finite difference front tracking method for two phase 1-D moving boundary problems［J］.Procedia Engineering，2015，127：1034-1040.

［9］林鴻熙，宋麗平，江良.美式期權(quán)自由邊界的計(jì)算方法［J］.數(shù)學(xué)的實(shí)踐與認(rèn)識(shí)，2014，44（24）：141-145.

［10］Gawlik E S，Lew A J.Unified analysis of finite elementmethods for problems with moving boundaries［J］.Siam Journal on Numerical Analysis，2015，53（6）：2822-2846.

［11］鄧建霞.有自由面滲流分析的加密高斯點(diǎn)單元傳導(dǎo)矩陣調(diào)整法［D］.成都：四川大學(xué)，2004.

［12］鐘艷麗，嚴(yán)月月，鐘思超，等.求解自由邊界問題的自適應(yīng)投影方法［J］.重慶工商大學(xué)學(xué)報(bào)（自然科學(xué)版），2017（5）：7-12.

［13］余春君.Signorini問題的無網(wǎng)格投影迭代算法［D］.重慶：重慶師范大學(xué)，2016.

［14］AliabadiM H，Ortiz E L.Numerical treatment ofmoving and free boundary value problems with the tau method［J］.Computers＆Mathematics with Applications，1998，35（8）：53-61.

［15］李樹忱，王兆清.高精度無網(wǎng)格重心插值配點(diǎn)法——算法、程序及工程應(yīng)用［M］.北京：科學(xué)出版社，2012.