基于谱方法和单纯形算法的一类偏微分方程参数反演研究.pdf

第 卷第期V o l ,N o 滨州学院学报J o u r n a l o fB i n z h o uU n i v e r s i t y 年月A p r,【微分方程与动力系统研究】基于谱方法和单纯形算法的一类偏微分方程参数反演研究收稿日期:基金项目:河北省高等学校科学研究计划项目(Z C )第一作者简介:王福昌(),男,山东定陶人,教授,硕士,主要从事应用数学研究.E m a i l:f z m a t h c o m王福昌,贺财宝(防灾科技学院 基础部,河北 三河 )摘要:根据观测数据反演偏微分方程参数具有重要的应用价值.通过基于快速傅立叶变换的谱方法实现对偏微分方程快速高精度求解,与观测数据结合建立待优化的目标函数,再用带边界约束的N e l d e r M e a d单纯形优化方法进行参数反演.通过算例证实了算法的有效性.关键词:快速傅立叶变换;谱方法;N e l d e r M e a d单纯形法;偏微分方程;参数反演中图分类号:O ;O 文献标识码:A D O I:/j c n k i 偏微分方程和偏微分方程组研究中一个具有实际意义的问题是,在确定了偏微分方程形式和获得了实际观测数据后,如何对偏微分方程中的参数进行有效的估计或解释 .一般地,偏微分方程组没有解析解或者解析解表达式非常复杂,因此主要考虑数值解.基于快速傅立叶变换的谱方法是一种求解某些偏微分方程的有效高精度算法,计算速度快、精度高 ,反演参数的迭代算法需要求大量的偏微分方程数值解,采用这一方法可以在一般电脑上方便快速地运行出结果.为了加快优化算法的计算速度,考虑观察数据噪声的影响,建立优化目标函数时,使用N e l d e r M e a d(NM)单纯形优化方法 .基本假设和模型实际的偏微分方程类型多而复杂,这里选取偏微分方程nutnL uN(u)()开展研究.其中,uu(x,t,p)为x,t的函数,pD为方程参数,D为参数范围,L代表线性算子,N(u)为非线性项.已知初始条件为u(x,t,p),在周期性边界条件下求某一时刻的u(x,t,p).偏微分方程一般情况下没有解析解,常用有限差分法求数值解.然而由于在有限差分法中,每个位置的导数由邻近点计算求得,故精度不够高.基于快速傅立叶变换的谱方法使用所有已知点计算某一处的导数,可以提高精度.另外,由于C F L条件限制,有限差分方法中时间步长不能太大,傅立叶谱方法是在先由快速傅立叶变换将偏微分方程变为常微分方程后,再利用变步长的常微分方程数值解法求解,因而能在较高精度下尽快得到数值结果,故本文采用该方法计算偏微分方程数值解,步骤如下.步骤对于形如式()的偏微分方程,可利用变量替换将nutn化为一阶方程组,若n,则忽略此步.步骤给定方程参数p后,在x域上对方程()做傅立叶变换得到ut(k)uFN(u),其中滨州学院学报第 卷uu(k,t)代表u(x,t)在x域上的傅立叶变换.为说明对线性算子L和非线性项N(u)做傅立叶变换的细节,可令线性算子Laxbxc(a,b,c为常数),N(u)uuuxf(x)ux.对线性部分有FL ua(ik)b(ik)cu,其中a(k)a(ik)b(ik)c.对非线性部分有FN(u)FFuFuF(i k)uf(x)Fiku .于是关于x的偏导数都通过傅立叶变换及逆变换简化为u(k,t)和k的代数运算,再将u和k离散化,可化偏微分方程为常微分方程.步骤利用时间步进法(如R u n g e K u t t a法等)数值计算离散化的关于u(k,t)的微分方程组,默认条件为周期性边界条件,初值条件u(k,t)Fu(x,t).步骤将上一步得到结果从频域变换为空间域,即u(x,t)Fu(k,t).注意,若空间再增加一个维度,基本步骤不变,只是用二维傅立叶变换替换掉一维傅立叶变换.另外,用傅立叶谱方法求解偏微分方程,隐含周期性边界条件.周期性边界条件可将具有时空周期性物理问题化为单元进行处理,对一些特定非周期性问题,需要修改其他条件来确保边界处的函数值恒为某一特定常数,以排除相邻周期间的干扰.参数反演方法 带边界约束的NM单纯形算法单纯形算法是利用单纯形各顶点的目标函数值,在连续改变几何图形的过程中,逐步以目标函数值较小的顶点取代目标函数值最大的顶点,从而求出最优点(最小值)的方法,属于直接法.假设求最小值的最优化问题模型为 m i nxDf(x),()其中,可行域DRn为下界是xLxL,xL,xn L、上界是xUxU,xU,xn U 的超矩形区域.下面给出该问题的改进NM单纯形算法.步骤给定初始点xx,x,xn,做中心化、反正弦和平移变换.为了利用边界xL,xU信息,先做中心化变换xcxxLxUxL,则每个分量变为到之间的实数,再做反正弦和平移变换,得到xu a r c s i n(m a x(,m i nxc),设置单纯形算法参数:反射系数、紧缩系数(,),扩张系数、收缩系数(,),设定精度.步骤用变换后的初始点xu在新区域G 内生成有n个顶点单纯形.构造单纯形的方法不唯一,这里让xu作为单纯形的个顶点,其余n个顶点xu,xu,xn u中的第i(i,n)个顶点xi u的第i(i,n)个非零分量是xu第i(i,n)个非零分量乘以 ,其余分量与xu相同.若xu第i(i,n)个分量为,则生成顶点xi u时,可以让加 ,即令xi u的第i个分量为 .步骤把n个单纯形顶点做反变换xi Ds i n(xi u)(xUxL)xL,i,n,并计算对应的函数值f(xi D),i,n.对上面得到的n个单纯形顶点进行反变换到原来可行域G中,计算对应目标函数值,并按目标函数值的大小重新编号,满足f(xs)f(xs)f(xn s)f(x(n)s)停止迭代,并输出xs作为最优点,否则执行下一步.步骤计算xnx(n)s(x(n)sxn s),若f(xn)f(xs)则执行步骤.否则,当f(xn)f(xn)时执行步骤,当f(xn)f(xn)时执行步骤.步骤计算xnx(n)s(xnx(n)s),若f(xn)f(xs),令xxn,执行步骤.步骤令xnxn,执行步骤.步骤令xnxi|f(xi)m i nf(xn),f(xn),计算xnxn(xnxn).若f(xn)第期王福昌,贺财宝基于谱方法和单纯形算法的一类偏微分方程参数反演研究f(xn),则令xnxn,执行步骤;否则,执行步骤.步骤令xixs(xi sxs),执行步骤.算法设计通过计算机模拟给出反演偏微分方程参数的方法检验方法的有效性.建立待优化的目标函数.先预设给出参数pe、时间t和空间x的分割方案,计算出u(t,x,pe)在这些节点处的值,然后加上随机扰动,模拟观测数据ur(t,x,pe);接着给定一组参数猜测值p,用傅立叶谱方法求出在这些节点处的值uc(t,x,p),与上面得到的观测数据做差,按照最小二乘准则建立目标函数ur(t,x,pe)uc(t,x,p),参数即为待优化的变量p,可以预先设定取值范围D,降低搜索难度.m i npDur(t,x,pe)uc(t,x,p).()对于待优化的多元函数(),先设定参数边界,再使用带边界约束的NM单纯形算法,即可求出参数p的最优值使得目标函数取得最小值.按照这种方法,选取了两个典型的偏微分方程组案例进行计算机模拟测试.计算机模拟测试通过二维黏性B u r g e r s方程和二维S c h n a k e n b e r g模型的参数反演来检验算法的有效性.二维黏性B u r g e r s模型参数反演B u r g e r s方程是流体力学中一个非常重要和基本的偏微分方程,应用广泛.二维黏性B u r g e r s方程形式为 ut(xy)uu(xy)u.()其中,u为速度,x和y为空间坐标,t为时间,和为模型参数.对式()做二维傅立叶变换,得ut(kxky)uFFu Fi(kxky)u .取参数,初始条件取为u(x,y,)s e c h(xy).图模拟了在t,时刻的u的等图形和加上随机扰动后的点.为了反演模型参数,对其加上随机扰动,图中“”即表示在不同时刻对u扰动后的观测值,目标是利用扰动后的随机点反演二维B u r g e r s模型中的个参数.为便于寻找,设定初始值猜测为 ,参数寻优的范围分别是 ,.图二维B u r g e r s方程的数值解与随机扰动的点根据扰动点的坐标数据,建立优化目标函数.目标函数中含有两个变量,利用边界约束的NM单纯形方法进行优化,得到的参数估计值为 ,与模型的参数值相比,非常接近.二维S c h n a k e n b e r g模型参数反演斑图是一类普遍存在于自然界、在时间或空间上具有某种规律的非均匀宏观结构.反应扩散系统是斑图理论中研究的最为广泛的系统,广泛应用于化学、生物学、物理学、医学和金融学等领域.S c h n a k e n b e r g模型是反应扩散系统中的一个有趣的模型,数学形式为滨州学院学报第 卷 ut(xy)u(auuv),vtd(xy)v(buv).()其中,u和v可以看作两种化学反应物质的浓度,x和y为空间坐标,t为时间,a,b,d,为常数.对式()做二维傅立叶变换,可将其转化为偏微分方程组ut(kxky)uFaFuFuFv ,vtd(kxky)vFbFuFv .参数取值为:a,b,d,.为了得到靶形波,设置初始条件:u在x y平面原点处为,其他位置处为;v在x y平面上均为.为了反演模型参数,在上述参数下计算t,时刻的u,并对其加上随机扰动,图中“”表示在不同时刻对u扰动后的观测值,目标是利用扰动后的随机点反演S c h n a k e n b e r g模型中的个参数.为便于寻找,设定初始值猜测为a ,b,d,参数寻优的范围分别是a ,b,d,.图二维S c h n a k e n b e r g方程的数值解与随机扰动的点经过编程计算,得到的参数估计为a ,b ,d ,与原来模型的参数值非常接近.结论本文针对两类偏微分方程(组)的参数估计问题进行了讨论,给出了一种比较有效的参数估计方法,这种方法也可以推广到更一般的类似情形.当然,这种方法是基于快速傅立叶变换的谱方法,对于复杂边界和谱方法不能使用的情形,需要将问题改造成能使用谱方法的问题或者使用其他方法快速求出偏微分方程(组)的数值解.一般情况下偏微分数值解计算较慢,为了加快优化算法计算速度和精度,可以在NM单纯形优化算法中借助专业知识,设置尽可能精确的参数范围和初值猜测值.为了克服初值的敏感性,也可尝试使用其他智能优化算法,只是要增加较多的寻优过程的计算时间.在建立优化目标函数时要注意两点:一是这里假设时间t和位置(x,y)是精确的,观测误差集中于u(t,x,y)和v(t,x,y),如果观测点的位置不在数值解的分割节点处,可以使用插值的方法来估计;二是要注意如果观测值被污染导致有离群值存在,会导致参数估计误差增大,这时考虑使用稳健的准则,如最小一乘或者删除离群值,这也是使用NM单纯形算法的原因,这种优化方法不要求优化的目标函数有光滑性.参考文献:I S AKOVV I n v e r s ep r o b l e m s f o rp a r t i a l d i f f e r e n t i a l e q u a t i o n sM N e wY o r k:S p r i n g e r,MA L E NG I E RB,VANK E E RR P a r a m e t e re s t i m a t i o ni nc o n v e c t i o nd o m i n a t e dn o n l i n e a rc o n v e c t i o n d i f f u s i o np r o b l e m sb yt h er e l a x a t i o nm e t h o da n dt h ea d j o i n t e q u a t i o nJ J o u r n a lo f c o m p u t a 第期王福昌,贺财宝基于谱方法和单纯形算法的一类偏微分方程参数反演研究t i o n a l a n da p p l i e dm a t h e m a t i c s,():D E NGXY,L I AOQX P a r a m e t e r e s t i m a t i o n f o rp a r t i a l d i f f e r e n t i a l e q u a t i o n sb yc o l l a g e b a s e dn u m e r i c a l a p p r o x i m a t i o nJ M a t h e m a t i c a l p r o b l e m s i ne n g i n e e r i n g,():MA R I NO,C ON S T ANT I N E S C U E,S M I TH B As c a l a b l em a t r i x f r e es p e c t r a le l e m e n ta p p r o a c hf o ru n s t e a d yP D Ec o n s t r a i n e do p t i m i z a t i o nu s i n gP E T S c/T AOJ J o u r n a lo fc o m p u t a t i o n a ls c i e n c e,:P E D R AMP,M I R Z A E IM,G OU S HEHSS U s i n gs p e c t r a lm e t h o da s a na p p r o x i m a t i o n f o r s o l v i n gh y p e r b o l i cP D E sJ C o m p u t e rp h y s i c sc o mm u n i c a t i o n s,(/):C ANUT O C,HU S S A I N IM Y,QUA R T E R ON IA,e ta l F u n d a m e n t a l so fs p e c t r a lm e t h o d sf o rP D E sM/S p e c t r a lm e t h o d s i nf l u i dd y n a m i c s B e r l i n:S p r i n g e r,N E L D E RJA,ME A DR AS i m p l e xm e t h o df o r f u n c t i o nm i n i m i z a t i o nJ T h ec o m p u t e r j o u r n a l,():L AGA R I A SCJ,R E E D SAJ,WR I GHT H M,e ta l C o n v e r g e n c ep r o p e r t i e so ft h eN e l d e r M e a ds i m p l e xm e t h o d i nl o wd i m e n s i o n sJ S I AMj o u r n a l o fo p t i m i z a t i o n,():P a r a m e t e r s I n v e r s i o nS t u d yo faK i n do fP a r t i a lD i f f e r e n t i a lE q u a t i o nB a s e do nS p e c t r a lM e t h o da n dN e l d e r M e a dS i m p l e xA l g o r i t h mWANGF u c h a n g,HEC a i b a o(D e p a r t m e n t o fB a s i cC o u r s e s,I n s t i t u t e o fD i s a s t e rP r e v e n t i o n,S a n h e ,C h i n a)A b s t r a c t:P a r a m e t e r s i n v e r s i o no fp a r t i a ld i f f e r e n t i a le q u a t i o nf r o mo b s e r v e dd a t ah a sg r e a tv a l u e s B a s e do nF F Ts p e c t r a lm e t h o da n dN e l d e r M e a ds i m p l e xa l g o r i t h m,ap r a c t i c a la l g o r i t h mi sp r e s e n t e d G i v e nt h eg u e s sv a l u e so fp a r a m e t e r s,t h ep r e d i c t i o nv a l u e so fp a r t i a l d i f f e r e n t i a l e q u a t i o n ss o l u t i o n sc a nb ec o m p u t e db y t h eF F Ts p e c t r a l a l g o r i t h m,a n d t h eo p t i m i z a t i o no b j e c t i v e f u n c t i o nc a nb e c o n s t r u c t e du s i n gd i f f e r e n c eb e t w e e no b s e r v e dd a t aa n dp r e d i c t i o nv a l u e s S o,t h eo p t i m a lp a r a m e t e r sa n do r d e ra r es e a r c h e db yN e l d e r M e a ds i m p l e xo p t i m i z a t i o nm e t h o d F i n a l l y,b yc o m p u t e rs i m u l a t i o n,as e to f f e a s i b l ei n v e r s i o nm e t h o d i sp r e s e n t e df o rs t u d ys c i e n t i f i ca n de n g i n e e r i n gp r o b l e m s K e y w o r d s:F F T;s p e c t r a lm e t h o d;N e l d e r M e a ds i m p l e xm e t h o d;p a r t i a ld i f f e r e n t i a lq u a t i o n;p a r a m e t e r s i n v e r s i o n(责任编辑:贾晶晶)