對(duì)流擴(kuò)散方程最優(yōu)控制問(wèn)題的重心插值配點(diǎn)格式

黃蓉, 姚夢(mèng)麗, 翁智峰

(華僑大學(xué) 數(shù)學(xué)科學(xué)學(xué)院, 福建 泉州 362021)

對(duì)流擴(kuò)散方程描述的最優(yōu)控制問(wèn)題被廣泛應(yīng)用于許多領(lǐng)域,如水污染處理[1]、空氣污染[2]等.對(duì)流擴(kuò)散方程模擬一個(gè)化學(xué)或生物過(guò)程,涉及的物種相互之間會(huì)發(fā)生擴(kuò)散、對(duì)流,因此,尋找穩(wěn)定、高效的數(shù)值求解方法具有十分重要的實(shí)際意義.目前,許多學(xué)者已經(jīng)提出許多數(shù)值求解格式,如間斷伽遼金方法[3]、雜交間斷伽遼金方法[4]、有限元方法[5-8]、雙線(xiàn)性偽譜方法[9]、勒讓德-伽遼金譜方法[10-11]、自適應(yīng)間斷伽遼金方法[12]、譜伽遼金近似方法等[13-15].本文主要研究由對(duì)流擴(kuò)散方程控制的最優(yōu)控制問(wèn)題.

1 預(yù)備知識(shí)

設(shè)區(qū)域Ω∈Rn(n=1,2)是帶有利普希茨邊界?Ω的空間有界域.考慮以下無(wú)約束最優(yōu)控制問(wèn)題,即

(1)

重心插值配點(diǎn)格式廣泛應(yīng)用于數(shù)值求解各類(lèi)微分方程,如平面彈性問(wèn)題[16]、Volterra積分方程[17]、Allen-Cahn方程[18-19]、Burgers方程[20].重心插值配點(diǎn)格式是一種新型的無(wú)網(wǎng)格方法,能以機(jī)器精度任意逼近光滑函數(shù),具有操作簡(jiǎn)單、計(jì)算有效、精度高等優(yōu)勢(shì).然而,對(duì)重心插值配點(diǎn)格式求解微分方程的研究相對(duì)較少.Yi等[21]采用重點(diǎn)插值配點(diǎn)格式求解時(shí)間分?jǐn)?shù)階電報(bào)方程,并給出理論分析.文獻(xiàn)[22-23]采用重心有理插值配點(diǎn)格式分別求解熱傳導(dǎo)方程、電報(bào)方程,并給出格式的誤差分析.Darehmiraki 等[24]基于重心插值配點(diǎn)格式,求解橢圓對(duì)流擴(kuò)散方程的最優(yōu)控制問(wèn)題,并證明了配點(diǎn)格式的收斂性.

2 重心插值配點(diǎn)格式

2.1 重心Lagrange插值配點(diǎn)格式

假設(shè)存在n+1個(gè)互異節(jié)點(diǎn)xj,函數(shù)f(x)在節(jié)點(diǎn)xj處的函數(shù)值為fj(j=0,1,…,n).令插值多項(xiàng)式p(x)在節(jié)點(diǎn)處成立,p(xj)=fj,根據(jù)多項(xiàng)式的唯一性,p(x)可改寫(xiě)為L(zhǎng)agrange插值形式,即

(2)

式(2)中:Lj(x)是Lagrange插值的基函數(shù),滿(mǎn)足基函數(shù)的性質(zhì),有

(3)

(4)

將式(4)代入式(2),可得

(5)

當(dāng)p(x)=1時(shí),有

(6)

結(jié)合式(5),(6),則重心Lagrange插值公式為

(7)

2.2 重心有理插值配點(diǎn)格式

根據(jù)一維重心插值配點(diǎn)格式的推導(dǎo)過(guò)程,對(duì)于定義在區(qū)間[a,b]上的函數(shù)f(x),給定n+1個(gè)插值節(jié)點(diǎn)a=x0

(8)

將zk(x)改寫(xiě)為L(zhǎng)agrange插值公式,即

(9)

(10)

結(jié)合式(8)～(10),重心有理插值公式為

(11)

2.3 重心插值配點(diǎn)格式的微分矩陣

設(shè)p(x)為函數(shù)f(x)的重心Lagrange插值公式,則p(x)關(guān)于x求導(dǎo),可得

(12)

3 對(duì)流擴(kuò)散方程最優(yōu)控制問(wèn)題的離散格式

3.1 最優(yōu)性條件

針對(duì)對(duì)流擴(kuò)散方程最優(yōu)控制問(wèn)題,采用Lagrange乘子求導(dǎo)數(shù)法,推出最優(yōu)控制問(wèn)題的最優(yōu)性條件.對(duì)于最優(yōu)控制問(wèn)題(1),定義p為區(qū)域Ω上的Lagrange乘子,即

(13)

對(duì)式(13)進(jìn)行泛函變分,對(duì)p求Frechet導(dǎo)數(shù),推導(dǎo)得出狀態(tài)方程,即

(14)

對(duì)v求Frechet導(dǎo)數(shù),導(dǎo)出伴隨方程為

(15)

類(lèi)似地,對(duì)u求Frechet導(dǎo)數(shù),則最優(yōu)性方程為

γu-p=0.

(16)

將對(duì)流擴(kuò)散最優(yōu)控制問(wèn)題轉(zhuǎn)化為代數(shù)方程組,即

(17)

3.2 對(duì)流擴(kuò)散方程最優(yōu)控制問(wèn)題的離散格式

采用重心插值配點(diǎn)格式離散方程組,求解對(duì)流擴(kuò)散方程最優(yōu)控制問(wèn)題的最優(yōu)化條件,選取區(qū)域Ω=[0,1]×[0,1],α=(α1,α2)T,則化簡(jiǎn)后的最優(yōu)控制系統(tǒng)為

(18)

令空間x,y方向的節(jié)點(diǎn)分別是M,N,設(shè)v(x,y),p(x,y)的重心插值為

(19)

式(19)中:ψi(x),φj(y)分別是空間x,y方向上的基函數(shù);vi,j=v(xi,yj);pi,j=p(xi,yj).

考慮v(x,y)對(duì)空間方向變量x,y求k+t階偏導(dǎo)數(shù),即

(20)

偏導(dǎo)數(shù)在節(jié)點(diǎn)(xp,yr)處的函數(shù)近似值為

(21)

類(lèi)似地,p(x,y)對(duì)空間方向變量x,y求k+t階偏導(dǎo)數(shù),則偏導(dǎo)數(shù)在節(jié)點(diǎn)(xp,yr)處的近似值為

(22)

將式(21),(22)代入式(18)中,則對(duì)流擴(kuò)散最優(yōu)控制系統(tǒng)的重心插值配點(diǎn)格式離散格式為

(23)

式(23)的微分矩陣形式為

(24)

3.3 相容性分析

設(shè)函數(shù)u(x,y)運(yùn)用重心Lagrange插值法逼近的數(shù)值解為pK,S(x,y),誤差函數(shù)e(x,y)為

e(x,y)=u(x,y)-pK,S(x,y).

(25)

重心Lagrange插值公式的逼近性質(zhì),如引理1所示.

引理1[21]若u(x,y)∈C(n+1)(Ω),Ω是非空、具有Lipschitz連續(xù)邊界的開(kāi)區(qū)域,則

(26)

類(lèi)似地,有

(27)

根據(jù)引理1,可得定理1.

定理1設(shè)對(duì)流擴(kuò)散方程最優(yōu)控制問(wèn)題中的狀態(tài)變量函數(shù)、伴隨變量函數(shù)分別是v(x,y),p(x,y),且v(x,y),p(x,y)∈C(n+1)(Ω),Ω=[a,b]×[c,d], 數(shù)值解分別是v(xk,ys),p(xk,ys),則

(28)

證明:定義線(xiàn)性微分算子D1,D2為

(29)

式(29)對(duì)應(yīng)的離散格式為

(30)

為簡(jiǎn)化分析過(guò)程,僅先分析與狀態(tài)量函數(shù)v(x,y)相關(guān)的項(xiàng),將式(29),(30)相減,可得

D1v(x,y)-D1(xk,ys)=-ε[v(x,y)-vx,x(xk,ys)]-ε[v(x,y)-vy,y(xk,ys)]+

α1[v(x,y)-vx(xk,ys)]+α2[v(x,y)-vy(xk,ys)]-

[v(x,y)-v(xk,ys)]=R1+R2+R3+R4+R5.

(31)

式(31)中:R1=-ε[v(x,y)-vx,x(xk,ys)];R2=-ε[v(x,y)-vyy(xk,ys)];R3=α1[v(x,y)-vx(xk,ys)];R4=α2[v(x,y)-vy(xk,ys)];R5=-[v(x,y)-v(xk,ys)].

根據(jù)引理1,有

(32)

類(lèi)似地,可得

(33)

(34)

同理,可推得

(35)

文獻(xiàn)[22]的定理的推導(dǎo)過(guò)程類(lèi)似定理1,采用重心有理配點(diǎn)格式求解二維最優(yōu)控制問(wèn)題的相容性誤差,即定理2.

定理2設(shè)v(x,y),p(x,y)分別是對(duì)流擴(kuò)散控制問(wèn)題中的狀態(tài)量函數(shù)、伴隨量函數(shù),v(x,y),p(x,y)∈C(n+1)(Ω),Ω=[a,b]×[c,d],v(xk,ys),p(xk,ys)分別是v(x,y),p(x,y)運(yùn)用重心有理配點(diǎn)格式求解的數(shù)值解,則

(36)

式(36)中:C=c·max{ε,α1,α2,(1/γ)},c為常數(shù);h1,h2分別是空間x,y方向的步長(zhǎng).

4 數(shù)值算例

4.1 算例1

為便于分析,定義最大相對(duì)誤差為

(37)

一維、二維的最優(yōu)控制問(wèn)題分別采用切比雪夫重心插值配點(diǎn)格式、有限差分法兩種離散方法,比較算例的數(shù)值結(jié)果,驗(yàn)證配點(diǎn)格式的有效性及高精度.對(duì)于一維最優(yōu)控制問(wèn)題,選取的真解為

v(x)=x(x-1)ex,p(x)=2x2(x-1)ex.

(38)

選取區(qū)域Ω=[0,1],ε=10-6,α=10-3,γ=10-4,令節(jié)點(diǎn)數(shù)為M=12,重心Lagrange插值配點(diǎn)格式、重心有理插值配點(diǎn)格式的數(shù)值解圖,分別如圖1,2所示.

(a) v(x) (b) p(x) (c) u(x)

(a) v(x) (b) p(x) (c) u(x)

由圖1,2可知:對(duì)于狀態(tài)量、伴隨量、控制量,采用兩種重心插值配點(diǎn)格式求解的數(shù)值解圖均逼近解析解圖,表明該數(shù)值算法是穩(wěn)定的.

分別選取ε,α,γ的不同剖分,重心Lagrange插值配點(diǎn)格式、重心有理配點(diǎn)格式、差分法的最大相對(duì)誤差(e),如表1,2,3所示.狀態(tài)量的收斂階(R)對(duì)比,如圖3所示.

表1 重心Lagrange插值配點(diǎn)格式的最大相對(duì)誤差(M)

表2 重心有理配點(diǎn)格式的最大相對(duì)誤差(M)

表3 差分法求解的最大相對(duì)誤差(M)

圖3 狀態(tài)量的收斂階對(duì)比(m)

由表1～3可知:采用重心插值配點(diǎn)格式求解變量比差分法求解時(shí)的誤差更小,剖分少量節(jié)點(diǎn),可達(dá)到格式高精度.當(dāng)選取的節(jié)點(diǎn)數(shù)相同時(shí),采用重心Lagrange插值比重心有理插值求解方程的誤差更小.由圖3可知:差分法格式的收斂階是二階,配點(diǎn)格式滿(mǎn)足指數(shù)收斂性質(zhì).

4.2 算例2

v(x,y)=sin(πx)sin(πy),p(x,y)=π2sin(πx)sin(πy).

(39)

(40)

選取區(qū)域Ω=[0,1]×[0,1],α=(1,1)T,γ=0.5,ε=0.1,令節(jié)點(diǎn)數(shù)為M=16,N=16,則重心Lagrange插值配點(diǎn)格式、重心整理插值配點(diǎn)格式的精確解、數(shù)值解與誤差,分別如圖4,5所示.由圖4,5可知:采用兩種重心插值配點(diǎn)格式求解狀態(tài)量、伴隨量的精確解圖像與數(shù)值解圖像逼近,且最大相對(duì)誤差精度高,可達(dá)到10-13量級(jí),表明兩種重心插值配點(diǎn)格式均穩(wěn)定.

(a) 狀態(tài)量的精確解 (b) 狀態(tài)量的數(shù)值解 (c) 狀態(tài)量的誤差

分別選取不同的剖分節(jié)點(diǎn)數(shù),重心Lagrange插值配點(diǎn)格式、重心有理配點(diǎn)格式、差分法的最大相對(duì)誤差,如表4,5,6所示.由表4,5可知:最大相對(duì)誤差可達(dá)到10-11量級(jí),表明兩種重心插值配點(diǎn)格式均具有高精度,且前者的求解效果略?xún)?yōu)于后者.隨著剖分變細(xì),3種變量的最大相對(duì)誤差在逐漸減少.由表6可知:選取M=64,N=64,狀態(tài)量、伴隨量、控制量的最大相對(duì)誤差分別達(dá)到10-3,10-4,10-4量級(jí);與經(jīng)典的差分法比較,重心插值配點(diǎn)格式選取更少的節(jié)點(diǎn),即可達(dá)到更高的精度.

表4 重心Lagrange插值配點(diǎn)格式最大相對(duì)誤差(M,N)

表5 重心有理配點(diǎn)格式的最大相對(duì)誤差(M,N)

表6 差分法求解的最大相對(duì)誤差(M,N)

(a) 狀態(tài)量的精確解 (b) 狀態(tài)量的數(shù)值解 (c) 狀態(tài)量的誤差

狀態(tài)量收斂階對(duì)比,如圖6所示.由圖6可知:采用差分法求解方程的收斂階為2階;重心Lagrange插值配點(diǎn)格式、重心有理配點(diǎn)格式的收斂階都呈現(xiàn)指數(shù)遞減的效果,前者的收斂效果優(yōu)于后者.

圖6 狀態(tài)量收斂階對(duì)比(m,n)

5 結(jié)束語(yǔ)

基于Lagrange乘子法,將對(duì)流擴(kuò)散最優(yōu)控制問(wèn)題轉(zhuǎn)化為由狀態(tài)方程、伴隨方程、最優(yōu)性方程三者聯(lián)立形成的代數(shù)方程組,再分別采用重心Lagrange插值配點(diǎn)格式、重心有理插值配點(diǎn)格式離散求解方程組中的狀態(tài)量v、伴隨量p,并對(duì)提出的配點(diǎn)格式進(jìn)行相容性誤差分析.數(shù)值實(shí)驗(yàn)結(jié)果表明,兩種重心插值配點(diǎn)格式均具有高精度的特性,選取切比雪夫節(jié)點(diǎn)時(shí)具有指數(shù)收斂的效果.此外,與經(jīng)典的有限差分格式相比,該配點(diǎn)格式在剖分較少的節(jié)點(diǎn)數(shù)時(shí),可達(dá)到很高的精度.