黃蓉, 姚夢(mèng)麗, 翁智峰
(華僑大學(xué) 數(shù)學(xué)科學(xué)學(xué)院, 福建 泉州 362021)
對(duì)流擴(kuò)散方程描述的最優(yōu)控制問題被廣泛應(yīng)用于許多領(lǐng)域,如水污染處理[1]、空氣污染[2]等.對(duì)流擴(kuò)散方程模擬一個(gè)化學(xué)或生物過程,涉及的物種相互之間會(huì)發(fā)生擴(kuò)散、對(duì)流,因此,尋找穩(wěn)定、高效的數(shù)值求解方法具有十分重要的實(shí)際意義.目前,許多學(xué)者已經(jīng)提出許多數(shù)值求解格式,如間斷伽遼金方法[3]、雜交間斷伽遼金方法[4]、有限元方法[5-8]、雙線性偽譜方法[9]、勒讓德-伽遼金譜方法[10-11]、自適應(yīng)間斷伽遼金方法[12]、譜伽遼金近似方法等[13-15].本文主要研究由對(duì)流擴(kuò)散方程控制的最優(yōu)控制問題.
設(shè)區(qū)域Ω∈Rn(n=1,2)是帶有利普希茨邊界?Ω的空間有界域.考慮以下無約束最優(yōu)控制問題,即
(1)
重心插值配點(diǎn)格式廣泛應(yīng)用于數(shù)值求解各類微分方程,如平面彈性問題[16]、Volterra積分方程[17]、Allen-Cahn方程[18-19]、Burgers方程[20].重心插值配點(diǎn)格式是一種新型的無網(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)控制問題,并證明了配點(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)可改寫為L(zhǎng)agrange插值形式,即
(2)
式(2)中:Lj(x)是Lagrange插值的基函數(shù),滿足基函數(shù)的性質(zhì),有
(3)
(4)
將式(4)代入式(2),可得
(5)
當(dāng)p(x)=1時(shí),有
(6)
結(jié)合式(5),(6),則重心Lagrange插值公式為
(7)
根據(jù)一維重心插值配點(diǎn)格式的推導(dǎo)過程,對(duì)于定義在區(qū)間[a,b]上的函數(shù)f(x),給定n+1個(gè)插值節(jié)點(diǎn)a=x0 (8) 將zk(x)改寫為L(zhǎng)agrange插值公式,即 (9) (10) 結(jié)合式(8)~(10),重心有理插值公式為 (11) 設(shè)p(x)為函數(shù)f(x)的重心Lagrange插值公式,則p(x)關(guān)于x求導(dǎo),可得 (12) 針對(duì)對(duì)流擴(kuò)散方程最優(yōu)控制問題,采用Lagrange乘子求導(dǎo)數(shù)法,推出最優(yōu)控制問題的最優(yōu)性條件.對(duì)于最優(yōu)控制問題(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) 類似地,對(duì)u求Frechet導(dǎo)數(shù),則最優(yōu)性方程為 γu-p=0. (16) 將對(duì)流擴(kuò)散最優(yōu)控制問題轉(zhuǎn)化為代數(shù)方程組,即 (17) 采用重心插值配點(diǎn)格式離散方程組,求解對(duì)流擴(kuò)散方程最優(yōu)控制問題的最優(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) 類似地,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) 設(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ù)邊界的開區(qū)域,則 (26) 類似地,有 (27) 根據(jù)引理1,可得定理1. 定理1設(shè)對(duì)流擴(kuò)散方程最優(yōu)控制問題中的狀態(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) 證明:定義線性微分算子D1,D2為 (29) 式(29)對(duì)應(yīng)的離散格式為 (30) 為簡(jiǎn)化分析過程,僅先分析與狀態(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) 類似地,可得 (33) (34) 同理,可推得 (35) 文獻(xiàn)[22]的定理的推導(dǎo)過程類似定理1,采用重心有理配點(diǎn)格式求解二維最優(yōu)控制問題的相容性誤差,即定理2. 定理2設(shè)v(x,y),p(x,y)分別是對(duì)流擴(kuò)散控制問題中的狀態(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). 為便于分析,定義最大相對(duì)誤差為 (37) 一維、二維的最優(yōu)控制問題分別采用切比雪夫重心插值配點(diǎn)格式、有限差分法兩種離散方法,比較算例的數(shù)值結(jié)果,驗(yàn)證配點(diǎn)格式的有效性及高精度.對(duì)于一維最優(yōu)控制問題,選取的真解為 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)格式滿足指數(shù)收斂性質(zhì). 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)格式均具有高精度,且前者的求解效果略優(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) 基于Lagrange乘子法,將對(duì)流擴(kuò)散最優(yōu)控制問題轉(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á)到很高的精度.2.3 重心插值配點(diǎn)格式的微分矩陣
3 對(duì)流擴(kuò)散方程最優(yōu)控制問題的離散格式
3.1 最優(yōu)性條件
3.2 對(duì)流擴(kuò)散方程最優(yōu)控制問題的離散格式
3.3 相容性分析
4 數(shù)值算例
4.1 算例1
4.2 算例2
5 結(jié)束語(yǔ)