1、收稿日期:2022 10 10;修订日期:2022 12 08作者简介:谢硕平(1999)，男，硕士研究生，研究方向为反问题建模与计算。基金项目:国家自然科学基金项目(11961002，11861007);江西省自然科学基金项目;江西省研究生创新专项资金项目(YC2022 s617)。*通信作者:胡彬(1982)，女，硕士，副教授，研究生导师，研究方向为数学物理反问题研究及数值计算。E mail:bhu ecut edu cn。第 41 卷第 1 期2023 年 2 月江西科学JIANGXISCIENCEVol 41 No 1Feb 2023doi:1013990/j issn1001 367
2、9 202301003基于离散特征系统的抛物方程源项反演方法谢硕平，胡彬*，张文，王梓鉴，黄雯(东华理工大学理学院，330013，南昌)摘要:主要研究离散特征系统下抛物方程源项反演的对数型正则化方法。首先，用有限差分法离散椭圆算子，利用分块矩阵的特点计算出椭圆算子的离散特征值和相应的特征向量;然后，将它们应用到抛物型方程源项反演的对数型正则化方法中。通过数值实验表明，对数型正则化方法可以通过离散特征值及其对应的特征向量成功实现。关键词:抛物方程;源项反演;椭圆算子;有限差分;特征系统中图分类号:O175 26文献标识码:A文章编号:1001 3679(2023)01 011 06Source
3、Term Inversion of Parabolic EquationBased on Discretecharacteristic SystemXIE Shuoping，HU Bin*，ZHANG Wen，WANG Zijian，HUANG Wen(School of Science，East China University of Technology，330013，Nanchang，PC)Abstract:This paper mainly studies the logarithmic regularization method for source term inversionof
4、 parabolic equations in discrete characteristic systems Firstly，the finite difference method is usedto discretize the elliptic operator，and the discrete eigenvalues and corresponding eigenvectors of theelliptic operator are calculated by using the characteristics of the block matrix;Then，they are ap
5、-plied to the logarithmic regularization method for source term inversion of parabolic equations Nu-merical experiments show that the logarithmic regularization method can be successfully implementedthrough discrete eigenvalues and their corresponding eigenvectorsKey words:parabolic equation;source
6、term inversion;elliptic operator;finite difference;charac-teristic system0引言抛物型方程反问题的研究是一种应用性极强的研究。过去的几十年中，由于抛物型偏微分方程在研究热传导、扩散等物理现象的巨大理论价值和在工程数学、物理等方向的巨大实用价值而受到了许多研究者的关注。例如，抛物型偏微分方程可用于探究石油地质过程中沉积岩系地温分布情况和有关参数识别1，多孔介质中多组分反应性溶质运移反问题2，地表水污染的点源识别3。然而，在实际问题中，往往需要通过边界数据或额外数据来解决方程系数、未知非齐次项、初始分布等问题，于是产生了各种抛
7、物型偏微分方程反问题4 15。这里主要考虑以下抛物型方程的反源问题:ut(x，t)=di，j=1xjaij(x)ux()i c(t)u+f(x)，(x，t)(0，T)，u(x，t)=0，(x，t)0，T，u(x，0)=0，x(1)其中 aij(x)C1()，c(t)C(0，T)且 c(t)0，aij(x)=aji(x)，1i，jd，且存在大于零的常数 和 ，使得0 di=1i2di，j=1aij(x)ijdi=1i2，d，x 。记Lu:=di，j=1xjaij(x)ux()i，x (2)则算子 L 是一致椭圆算子。本文考虑的源项反演问题是:在定解问题(1)中，给定物理场 u(x，t)在 T 时
8、刻的测量值g()x L2()，由此重建未知源项 f()x，其中测量数据 g()x 满足g()x g()x L2()(3)这里，g()x=ux，()T 是精确值，为误差水平。众所周知，由于在实际测量数据中一定存在的微小误差会导致在重建未知源时产生不可控的高振荡误差，所以抛物型方程的反源问题是不适定的。对此，可以采用一些正则化方法来解决逆源问题。这一问题受到了各个领域研究者的广泛关注，他们采用了不同的方法。例如，有作者基于与 Lebesgue 测量相关联的傅里叶变换提出采用傅里叶正则化方法对热源进行了识别，并得到了误差估计4;还有作者提出利用简化的 Tikhonov正则化方法5，从数据 g()x
9、中确定热源 f()x，得到了正则解与精确解之间的 Hlder 型稳定性估计;在二维热源的确定中6，作者利用截断积分和傅里叶变换的方法，构造正则化解并得到显式误差估计;抛物方程中与空间相关的未知源的数值重构还有很多7 10。本文的安排如下，第 1 节介绍源项反演的对数型正则化方法，第 2 节给出如何求解算子离散后矩阵的特征值和对应特征向量的过程，第 3 节给出了 2 个数值算例，以说明通过使用离散特征值和相应特征向量恢复数值解的效果，最后得出了一些结论。1源项反演的对数型正则化方法根据文献 16，对于负一致对称椭圆算子 L满足以下引理。引理 1:负一致对称椭圆算子 L 的特征值有:1)L 的每个
10、特征值都是实数，且只有可数个特征值(其中重特征值按重数计)。2)L 的所有特征值 kk=1可排列成(重特征值按重数重复排列)0 123，且当 k 时，k。3)存 在 L2()中 的 一 组 标 准 正 交 基 wkk=1，使得 wk H10()是对应于特征值 k的特征函数，即有Lwk=kwkin wk=0on。定理1:设 f(x)L2()，则定解问题(1)的解可表达成u(x，t)=k=1fkektC(t)t0eks+C(s)dswk(x)，(4)其中 k和 wk(x)是负一致对称椭圆算子 L 特征值与特征函数，C(t)=t0c(s)ds，fk=(f，wk)。注:该定理可直接由文献 16 的 7
11、 1 节的结论得到。对于 精 确 的 终 端 数 据 g()x=ux，()T L2()，令 gk=(g，wk)。于是，由g(x)=u(x，T)=k=1fkekTC(T)T0eks+C(s)dswk(x)，(5)则有gk=fkekTC(T)T0eks+C(s)ds。于是得到fk=gkT0ek(Ts)e(C(T)C(s)ds，从而得到21江西科学2023 年第 41 卷f(x)=k=1gkT0ek(Ts)e(C(T)C(s)dswk(x)(6)对于带噪声的测量数据 g(x)，源项对应的反演结果记为 f(x)，且f(x)=k=1gkT0ek(Ts)e(C(T)C(s)dswk(x)，(7)其中 gk
12、=(g，wk)。由积分中值定理，可得limkT0ek(Ts)e(C(T)C(s)ds=limke(C(T)C()T0ek(Ts)ds=limke(C(T)C()1 ekTk=0。(8)因此，在利用式(7)反演源项时，测量数据 g(x)中的噪声将被趋于无穷的特征值 k放大，即源项反演是不适定的。为了克服源项反演的不适定性，近期针对一类热传导方程的热源反演问题提出了一种对数型正则化方法11。根据该对数型正则化方法的思想，即可得本文源项反演的正则化解为f，(x)=k=1gkT0ek(Ts)e(C(T)C(s)dsT0ek(Ts)e(C(T)C(s)()ds2+ln(1+k)wk(x)，(9)其中 0
13、 为一给定常数，为正则化参数。但是，在利用式(18)反演源项 f(x)得到其近似的正则化解时，必须要事先知道椭圆算子的特征系统，这对于定义在 上的一般椭圆算子是非常困难的。本文给出的正则化解只需形如(2)简单形式的椭圆算子的特征系统，且它利于离散为矩阵的形式并方便计算出离散情形下的特征系统，即连续特征系统的近似;然后，利用离散特征系统构建源项反演的正则化解。2离散椭圆算子特征系统的计算方法为简单起见，以下仅考虑 aij(x)=1，且分别以 =(0，1)和 =(0，1)(0，1)为例予以说明离散特征系统的计算方法。对于 =(0，1)，在等距网格点xixi=ih，h=1n，i=0，1，2，n上离散
14、一维椭圆算子 Lu=2ux2为2ux2u(xi1)2u(xi)+u(xi+1)h2。令 uiu(xi)。结合齐次边界条件和上述离散方法，则负椭圆算子 L 的特征值问题离散为矩阵特征值问题U=U，(10)其中=1h22112112112(n1)(n1)，U=u1u2un2un1。对于 =(0，1)(0，1)，在等距网格点(xi，yj)|xi=ih，yj=jh，h=1n，i，j=0，1，2，n)上用五点有限差分格式离散二维椭圆算子 Lu=2ux2+2uy2为2ux2u(xi1，yj)2u(xi，yj)+u(xi+1，yj)h2u(xi，yi1)2u(xi，yi)+u(xi，yi+1)h2。令 ui
15、j u(xi，yj)，U=(u11，u21，un1，1，u12，u22，u1，n1，u2，n1，un1，n1)T。结合齐次边界条件和上述离散方法，则二维负椭圆算子 L 的特征值问题离散为矩阵特征值问题AU=U，(11)其中 I 是(n 1)阶单位矩阵，A 是(n 1)2阶矩阵，B 是(n 1)阶矩阵，且A=1h2B I IB I IB I IB，31第 1 期谢硕平等:基于离散特征系统的抛物方程源项反演方法B=4 1 14 1 14 114。(12)实际上，矩阵特征问题(10)和(11)容易由下述引理所给方法计算其特征值和相应的特征向量。引理 212:设 C 是 M M 的三对角矩阵，且C=a
16、bcabcabca，则 C 的特征值为 i=a+2bccosiM+()1，对应的特征向量是ui=(cb)1/2sin(iM+1)，(cb)sin(2iM+1)，(cb)M/2sin(MiM+1)T，i=1，2，M。(13)引理 317:设矩阵 A 和 B 如式(12)所定义，矩阵 Q 定义为Q=0 1 10 1 10 110(n1)(n1)。设矩阵 B 的特征值和特征向量分别为 i和ui，i=1，2，n 1，Q 的特征值和特征向量分别为 j和 vj，j=1，2，n 1。显然，它们均可由引理 2 中的方法计算得到。则矩阵 A 的特征值为 k=i+j，对应的特征向量 wk是矩阵 uivTj按列拉直的列向量，k=i+(j 1)(n 1)。根据引理2 可得，一维椭圆算子 Lu 离散后矩阵 的特征值和对应的特征向量分别为i=2h21+cosi()n，i=1，2，n 1，(14)ui=sinin，sin2in，sinn()1i()nT，i=1，2，n 1。(15)再根据引理 3 可求出二维椭圆算子 Lu 离散后矩阵的特征值为k=2h22+cosin+cosj()n，k=i+(j 1)(n 1)，i
