热传导反问题在航空航天, 机械工程, 冶金以及建筑消防等领域具有广泛的应用. 一般来说热传导反问题分为以下几类: 逆时热传导问题, 反边值问题, 反系数问题, 反边界问题及反源问题, 目前对于此类问题的研究已有一些重要进展. 在文献[1]中, 作者结合边界元和Tikhonov正则化方法研究了随时间变化的热传导方程反源问题. 文献[2]采用动力系统方法求解热传导方程的反问题, 并与Tikhonov正则化方法进行了比较. 文献[3]使用Tikhonov正则化和截断奇异值分解正则化方法, 研究了热传导方程的源项识别问题, 利用确定位置的测量数据反演未知热源强度. 在文献[4]中, 作者比较了两种求解热传导方程初始值反问题的数值方法. 文献[5]研究了通过圆盘上终值时刻的温度场观测数据来反演初始温度的逆问题, 得到了初始数据和观测数据之间的积分表示关系. 文献[6]作者利用Landweber迭代正则化方法研究了热传导方程反源问题, 在最优控制框架的基础上, 得到了目标函数最小化的必要条件. 文献[7]中, 作者利用谱配置方法与正则化技术相结合的数值方法求解了四分之一平面内热传导反问题. 文献[8]研究了二维逆时热传导反问题的改进正则化方法, 得到了H$ \ddot{\rm o} $lder型的误差估计. 文献[9]利用第一类Chebyshev展开提出了一种非线性逆时热传导反问题求解的积分方程方法. 文献[10]中, 作者研究了圆柱区域上非定常热传导方程的正反问题, 分析了热传导系数反演问题的解关于数据测量误差和热耦合误差的敏感性. 文献[11]中, 作者根据温度和外部边界曲线上的热通量, 重构了双连通有界域的内部边界曲线, 开发了一种迭代算法求解积分方程, 获得了相应积分算子的Fréchet导数, 并且证明了线性化方程的唯一可解性, 通过三角积分法实现完全离散化, 结合Tikhonov正则化求解了反问题. 文献[12]中, 作者利用广义和修正的广义Tikhonov正则化方法求解了时间分数阶扩散方程的逆源问题. 文献[13]中, 作者提出了一种无需使用初始数据即可识别二维热传导方程中的运动边界的算法. 尽管目前针对不同类型的热传导反问题已有不少结果, 但关于三维热传导反问题的分析和数值研究相对较少, 特别是非规则区域和非均匀介质中的相关的热传导反问题还有待进一步深入研究.

该文主要研究一类三维矩形截面梁结构热传导反问题的数值算法, 此问题源于建筑物梁结构受热问题的研究. 我们分析了规则区域上的三维热传导反问题, 利用分离变量方法建立任意时刻温度场与初始温度场之间的对应关系, 给出明确的反演公式. 为了避免“逆问题犯罪”有效求解反问题, 对于正问题求解本文使用有限元对空间离散, 借助有限差分法对时间离散, 简要给出了热传导问题的数值计算方法, 对不同时刻节点处的温度进行数值求解. 针对数值求解反演方程的不适定性, 我们分别采用了Tikhonov正则化和终值数据扰动正则化方法反演初始温度场, 并从数值上分析了正则化参数以及数据观测时间对反演结果的影响. 结论表明, 随着时间推移温度场逐渐衰减, 此时噪声数据(误差水平)带来的影响逐渐增大, 反演效果逐渐变差, 甚至完全得不到初始温度场, 这种反演结果和实际物理过程一致. 论文结构如下: 第2, 3节给出三维热传导正问题的有限元解法及数值实验; 第4节中研究三维逆时热传导问题, 证明了反问题局部稳定性, 第5节给出了反问题求解的两种数值方法, 并通过数值结果验证了算法的有效性, 最后在第6节中给出了总结.

三维热传导问题的控制方程可以表示为^{[14, 15]}

(2.1) $ \begin{equation} \frac{\partial}{\partial x}\left(k\frac{\partial u}{\partial x}\right)+\frac{\partial}{\partial y}\left(k\frac{\partial u}{\partial y}\right)+\frac{\partial}{\partial z}\left(k\frac{\partial u}{\partial z}\right)+G=\rho c\frac{\partial u}{\partial t}, \quad 0<t<T, \end{equation} $

其中$ k $为热传导系数, $ G $为内热源, $ \rho $为介质密度, $ c $为介质的比热容, $ u(x, y, z, t) $为$ t $时刻介质某点处的温度, $ \Omega $为三维空间上的有界区域, $ T $为正的常数.

初始条件为

(2.2) $ \begin{equation} u(x, y, z, t)|_{t=0}=g(x, y, z), \end{equation} $

边界条件为

(2.3) $ \begin{equation} \frac{\partial u(x, y, z, t)}{\partial v}=h(u_f-u)|_{\partial \Omega}, \end{equation} $

其中$ v=({{\boldsymbol{{ l }}}}_x, {{\boldsymbol{{ l }}}}_y, {{\boldsymbol{{ l }}}}_z) $是边界上单位外法向量, $ u_f $为外界温度, $ u $为物体温度, $ h $为对流换热系数, 当对流换热系数为零时对应第二类边界条件, 当此系数无穷大时对应第一类边界条件.

该文所研究的正问题是已知初始条件, 边界条件以及相关物理参数, 求不同时刻空间$ \Omega $中温度场的分布. 下面简述正问题的数值求解方法, 首先将空间进行有限元离散, 将温度函数表示为

(2.4) $ \begin{equation} u(x, y, z, t)=\sum\limits_{i=1}^{n}N_i(x, y, z)u_i(t), \end{equation} $

其中$ N_i $为单元形状函数, $ n $为单元中的节点数, $ u_i(t) $为随时间变化的节点温度. 通过Galerkin方法可以将方程(2.1)写为如下形式

(2.5) $ \begin{equation} \int_{\Omega}N_i\left[\frac{\partial}{\partial x}\left(k\frac{\partial u}{\partial x}\right)+\frac{\partial}{\partial y}\left(k\frac{\partial u}{\partial y}\right)+\frac{\partial}{\partial z}\left(k\frac{\partial u}{\partial z}\right)+G-\rho c\frac{\partial u}{\partial t}\right]{\rm d}\Omega=0, \ i=1, 2, \cdots , n. \end{equation} $

利用高斯公式得

(2.6) $ \begin{eqnarray} &&\int_{\Omega}N_i\left(G-\rho c\frac{\partial u}{\partial t}\right){\rm d}{\Omega}-\int_{\Omega}\left(k\frac{\partial N_i}{\partial x}\frac{\partial u}{\partial x}+k\frac{\partial N_i}{\partial y}\frac{\partial u}{\partial y}+k\frac{\partial N_i}{\partial z}\frac{\partial u}{\partial z}\right){\rm d}\Omega{}\\ &&+ \int_{\Gamma}N_ik\left(\frac{\partial u}{\partial x}{{{\boldsymbol{ l }}}}_x+\frac{\partial u}{\partial y}{{{\boldsymbol{ l }}}}_y+\frac{\partial u}{\partial z}{{{\boldsymbol{ l }}}}_z\right){\rm d}\Gamma=0, \end{eqnarray} $

这里$ \Gamma $表示区域$ \Omega $的边界, 注意

(2.7) $ \begin{equation} \int_{\Gamma}N_i\left(\frac{\partial u}{\partial x}{{{\boldsymbol{ l }}}}_x+\frac{\partial u}{\partial y}{{{\boldsymbol{ l }}}}_y+\frac{\partial u}{\partial z}{{{\boldsymbol{ l }}}}_z\right){\rm d}\Gamma=\int_{\Gamma}N_ih(u_f-u){\rm d}\Gamma, \end{equation} $

利用公式(2.4), (2.6)和(2.7), 将(2.5)式转化为

(2.8) $ \begin{eqnarray} &&\sum\limits_{j=1}^{n}\left( \int_{\Omega}k(\nabla N_i\cdot \nabla N_j){\rm d}\Omega+\int_{\Gamma}khN_iN_j{\rm d}\Gamma\right)u_j(t)+\sum\limits_{j=1}^{n}\left(\int_{\Omega}\rho cN_iN_j{\rm d}\Omega \right)\frac{\partial u_j(t)}{\partial t}{}\\ &=&\int_{\Omega}N_iG{\rm d}\Omega-\int_{\Gamma}N_ikhu_f{\rm d}\Gamma, \quad i=1, 2, \cdots , n, \end{eqnarray} $

其中$ i, j $为结点, $ \nabla $是梯度算子, 将(2.4)式写为如下形式

(2.9) $ \begin{equation} u={{{\boldsymbol{ N }}}}{{{\boldsymbol{ u }}}}= \left[\begin{array}{cccc} N_1&N_2&\cdots&N_n \end{array}\right] \left[\begin{array}{cccc} u_1\\u_2\\\vdots\\u_n \end{array}\right], \end{equation} $

由(2.9)式可得, 温度梯度矩阵为

(2.10) $ \begin{equation} {{{\boldsymbol{ g }}}} = \left[\begin{array}{cccc} { } \frac{\partial u}{\partial x}\\ { }\frac{\partial u}{\partial y}\\ { }\frac{\partial u}{\partial z} \end{array}\right] = \left[\begin{array}{cccc} { } \frac{\partial N_1}{\partial x}&{ }\frac{\partial N_2}{\partial x}&\cdots&{ }\frac{\partial N_n}{\partial x}\\ { } \frac{\partial N_1}{\partial y}&{ }\frac{\partial N_2}{\partial y}&\cdots&{ }\frac{\partial N_n}{\partial y}\\ { } \frac{\partial N_1}{\partial z}&{ } \frac{\partial N_2}{\partial z}&\cdots&{ }\frac{\partial N_n}{\partial z} \end{array}\right] \left[\begin{array}{cccc} u_1\\u_2\\\vdots\\u_n \end{array}\right] = [{{{\boldsymbol{ B }}}}]{{{\boldsymbol{ u }}}}, \end{equation} $

考虑到

(2.11) $ \begin{equation} {{{\boldsymbol{ g }}}}^T[{{{\boldsymbol{ D }}}}]{{{\boldsymbol{ g }}}}= \left[\begin{array}{cccc} { }\frac{\partial u}{\partial x}&{ }\frac{\partial u}{\partial y}&{ }\frac{\partial u}{\partial z} \end{array}\right] \left[\begin{array}{cccc} k&0&0\\0&k&0\\0&0&k \end{array}\right] \left[\begin{array}{cccc} { } \frac{\partial u}{\partial x}\\ { }\frac{\partial u}{\partial y}\\ { }\frac{\partial u}{\partial z} \end{array}\right] =k\left[\left(\frac{\partial u}{\partial x}\right)^2+\left(\frac{\partial u}{\partial y}\right)^2+\left(\frac{\partial u}{\partial z}\right)^2\right], \end{equation} $

其中$ {{\boldsymbol{{ D }}}}=k{{\boldsymbol{{ I }}}} $, $ {{\boldsymbol{{ I }}}} $表示单位矩阵, (2.8)式的矩阵形式为

(2.12) $ \begin{equation} [{{{\boldsymbol{ C }}}}]\left(\frac{\partial{{{\boldsymbol{ u }}}}}{\partial t}\right)+[{{{\boldsymbol{ K }}}}]{{{\boldsymbol{ u }}}}=\{{{{\boldsymbol{ f }}}}\}, \end{equation} $

其中

(2.13) $ \begin{equation} [{{{\boldsymbol{ C }}}}]=\int_{\Omega}\rho c[{{{\boldsymbol{ N }}}}]^T[{{{\boldsymbol{ N }}}}]{\rm d}\Omega , \end{equation} $

(2.14) $ \begin{equation} [{{{\boldsymbol{ K }}}}]=\int_{\Omega}[{{{\boldsymbol{ B }}}}]^T[{{{\boldsymbol{ D }}}}][{{{\boldsymbol{ B }}}}]{\rm d}\Omega+\int_{\Gamma} kh[{{{\boldsymbol{ N }}}}]^T[{{{\boldsymbol{ N }}}}]{\rm d}\Gamma \end{equation} $

和

(2.15) $ \begin{equation} \{{{{\boldsymbol{ f }}}}\}=\int_{\Omega}G[{{{\boldsymbol{ N }}}}]^T{\rm d}\Omega+\int_{\Gamma}khu_f[{{{\boldsymbol{ N }}}}]^T{\rm d}\Gamma, \end{equation} $

上述方程仅在空间中离散, 还需要对时间项进行离散, 这里使用有限差分法. 将$ [0, T] $进行等步长剖分, 步长为$ \Delta t $, $ t_n=n\Delta t $, 根据泰勒公式, 可以用$ t_n $时刻的温度$ u^n $将$ t_{n+1} $时刻的温度$ u^{n+1} $近似表示出来

(2.16) $ \begin{equation} \frac{\partial u^n}{\partial t}=\frac{u^{n+1}-u_n}{\Delta t}+o(\Delta t). \end{equation} $

引入参数$ \theta(0<\theta<1) $, 记

(2.17) $ \begin{equation} u^{n+\theta}=\theta u^{n+1}+(1-\theta)u^n, \ \ \ \{{{{\boldsymbol{ f }}}}\}^{n+\theta}=\theta \{{{{\boldsymbol{ f }}}}\}^{n+1}+(1-\theta)\{{{{\boldsymbol{ f }}}}\}^n. \end{equation} $

将(2.12)式离散后得

(2.18) $ \begin{equation} [{{{\boldsymbol{ C }}}}]\left(\frac{u^{n+1}-u^n}{\Delta t}\right)+[{{{\boldsymbol{ K }}}}](\theta u^{n+1}+(1-\theta)u^n)=\{{{{\boldsymbol{ f }}}}\}^{n+\theta}, \end{equation} $

整理得到全离散格式

(2.19) $ \begin{equation} ([{{{\boldsymbol{ C }}}}]+\theta\Delta t[{{{\boldsymbol{ K }}}}])[{{{\boldsymbol{ u }}}}]^{n+1}=([{{{\boldsymbol{ C }}}}]-(1-\theta)\Delta t[{{{\boldsymbol{ K }}}}])[{{{\boldsymbol{ u }}}}]^n+\Delta t\left(\theta \{{{{\boldsymbol{ f }}}}\}^{n+1}+(1-\theta)\{{{{\boldsymbol{ f }}}}\}^n\right). \end{equation} $

给定三维空间中的一个正方体$ \Omega\in[0, \pi]^3 $, 其密度$ \rho=1 $, 比热容$ c=1 $, 热传导系数$ k=0.04 $, 无内热源, 初始温度为$ g(x, y, z)=\sin x\sin y\sin z $, 考虑第一类边界条件的情形, 即$ u|_{\partial\Omega=0} $. 首先对空间$ \Omega $进行剖分, 沿$ x, y, z $轴方向的步长分别为$ h_1=h_2=h_3=\pi/J $, $ J\in Z^+ $为网格剖分数, $ x_i=ih_1, y_j=jh_2, z_k=kh_3 $, 时间$ t\in[0, T] $, 时间步长$ \Delta t=0.1 $, 定义相应的离散范数

(3.1) $ \begin{equation} \|v\|_{L^\infty}=\max\limits_{1\leq i, j, k\leq J-1}|v(x_i, y_j, z_k)|, \quad \|v\|_{L^2}=\left[\sum\limits_{i, j, k=1}^{J-1}|v(x_i, y_j, z_k)|^2\right]^{\frac{1}{2}}. \end{equation} $

利用公式(2.19)进行求解, 可求出内部所有节点的温度$ u_e(x_i, y_j, z_k, t_n) $, $ u(x_i, y_j, z_k, t_n) $为内部所有节点的精确温度, 所得数值解与精确解的误差比较见表 1和图 1. 其中

$ e_1(x_i, y_j, z_k, T):=u_e(x_i, y_j, z_k, T)-u(x_i, y_j, z_k, T), $

$ \|e_2(x_i, y_j, z_k, T)\|:=\frac{\|e_1(x_i, y_j, z_k, T)\|_{L^2}}{\|u(x_i, y_j, z_k, T)\|_{L^2}}, $

$ error1(x_i, y_j, z_k, T):=|e_1(x_i, y_j, z_k, T)|, $

$ error2(x_i, y_j, z_k, T):=\frac{|e_1(x_i, y_j, z_k, T)|}{|u(x_i, y_j, z_k, T)|}. $

对于$ J=10, 20, 30, 40 $这四种情形, 表 1中分别给出了不同时刻的$ L^2 $ -范数和$ L^\infty $ -范数意义下的误差, 图 1表示当$ h_1=h_2=h_3=\pi/20 $, $ \Delta t=0.1 $, $ T=6 $时, 精确解, 数值解和误差图像. 数值实验表明正问题数值解法是正确和有效的.

逆时热传导问题的是通过观测数据$ u(x, y, z, T) $求初始温度$ u(x, y, z, 0) $的分布^{[16, 18]}. 该问题通常是不适定的, 即问题的解不连续依赖于输入数据, 故需要使用正则化方法. 本文考虑由

(4.1) $ \begin{equation} \left\{\begin{array}{ll} u_t=a^2(u_{xx}+u_{yy}+u_{zz}), & (x, y, z, t)\in\Omega\times(0, T), \\ u|_{\partial\Omega} =0, & t\in(0, T), \\ u(x, y, z, 0)=g(x, y, z), & (x, y, z) \in \Omega \end{array}\right. \end{equation} $

产生的逆时问题. 利用$ u(x, y, z, T)=f(x, y, z) $的扰动值反演初始温度分布$ g(x, y, z) $的近似值, 或者反演$ u(x, y, z, t) $, 其中$ 0<t<T $. 记精确数据$ f $的扰动数据为$ f^{\delta} $, 且满足

(4.2) $ \begin{equation} \|f^{\delta}-f\|_{L^2(\Omega)}\leq\delta. \end{equation} $

利用分离变量法, 可得

(4.3) $ \begin{equation} {{\bf K}_{\bf 1}^{{\bf T-t}}}\cdot u(\cdot, t)(x, y, z)=\int_{\Omega}K_1(x, y, z, \xi, \eta, \zeta, T-t)u(\xi, \eta, \zeta, t){\rm d}\xi {\rm d}\eta {\rm d}\zeta =u(x, y, z, T), \end{equation} $

其中核函数为

(4.4) $ \begin{eqnarray} &&K_1(x, y, z, \xi, \eta, \zeta, t){}\\ &=&\frac{8}{\pi^3}\sum\limits_{l=1}^{\infty}e^{-a^2l^2t}\sin lx\sin l\xi \sum\limits_{m=1}^{\infty}e^{-a^2m^2t}\sin my\sin m\eta \sum\limits_{n=1}^{\infty}e^{-a^2n^2t}\sin nz\sin n\zeta. \end{eqnarray} $

为证明反问题的条件稳定性, 我们令$ {{\boldsymbol{{ X }}}}=(x, y, z) $, $ {{\boldsymbol{{ Y }}}}=(\xi, \eta, \zeta) $, 根据本征函数解法, 可以将正问题的解简写为

(4.5) $ \begin{equation} u({{\boldsymbol{{ X }}}}, t)=\sum\limits_{r=1}^\infty G_ru_r({{\boldsymbol{{ X }}}})e^{-\lambda_r a^2t}, \end{equation} $

这里常数$ G_r $为

(4.6) $ \begin{equation} G_r=\int_{\Omega}g({{\boldsymbol{{ X }}}})u_r({{\boldsymbol{{ X }}}}){\rm d}{{\boldsymbol{{ X }}}}, \end{equation} $

这里的$ u_r({{\boldsymbol{{ X }}}}) $和$ \lambda_r $分别是Laplace算子的第$ r $个正交本征函数和特征值.

令$ t'=t_0-t $, $ w({{\boldsymbol{{ X }}}}, t')=u({{\boldsymbol{{ X }}}}, t_0-t) $, 逆时问题等价于

(4.7) $ \begin{equation} \left\{\begin{array}{ll} w_{t'}=-a^2\Delta w, & ({{\boldsymbol{{ X }}}}, t')\in\Omega\times(0, T), \\ w|_{\partial\Omega} =0, & t'\in(0, T), \\ w|_{t'=0}=u({{\boldsymbol{{ X }}}}, t_0), & {{\boldsymbol{{ X }}}} \in \Omega. \end{array}\right. \end{equation} $

该问题的解可用$ \{u_r({{\boldsymbol{{ X }}}}), \lambda_r\} $表示为

(4.8) $ \begin{equation} ({\bf K}_1^{\bf T})^{-1}u({{{\bf{ X }}}}, t_0)=g({{\boldsymbol{{ X }}}})=w({{\boldsymbol{{ X }}}}, t_0)=\sum\limits_{r=1}^{\infty}w_r u_r({{\boldsymbol{{ X }}}})e^{\lambda_ra^2t_0}, \end{equation} $

这里常数$ w_r $为

(4.9) $ \begin{equation} w_r=\int_{\Omega}u({{\boldsymbol{{ X }}}}, t_0)u_r({{\boldsymbol{{ X }}}}){\rm d}{{\boldsymbol{{ X }}}}. \end{equation} $

定义

(4.10) $ \begin{equation} Q(M)=\left\{g({{\boldsymbol{{ X }}}}):g\in L^2(\Omega), \sum\limits_{r=1}^{\infty}G_r^2e^{2\lambda_ra^2T}\leq M^2\right\}, \end{equation} $

$ M $为已知正的常数, 要建立反问题的条件稳定性, 需对未知函数集做一个先验假设.

定理4.1  对于$ g({{\boldsymbol{{ X }}}})\in Q(M) $, (4.1)式的解满足

(4.11) $ \begin{equation} \|g(\cdot)\|_{L^2(\Omega)}\leq M^{t_0/T}\|u(\cdot, t_0)\|_{L^2(\Omega)}^{1-t_0/T}, \quad \forall t_0\in(0, T). \end{equation} $

证  类似(4.8)式可得

(4.12) $ \begin{equation} w({{\boldsymbol{{ X }}}}, T)=\sum\limits_{r=1}^{\infty}w_ru_r({{\boldsymbol{{ X }}}})e^{\lambda_ra^2T}, \end{equation} $

这里常数$ w_r $为

(4.13) $ \begin{equation} w_r=\int_{\Omega}u({{\boldsymbol{{ Y }}}}, t_0)u_r({{\boldsymbol{{ Y }}}}){\rm d}{{\boldsymbol{{ Y }}}}= \sum\limits_{s=1}^{\infty}G_se^{-\lambda_sa^2t_0}\int_{\Omega}u_s({{\boldsymbol{{ Y }}}})u_r({{\boldsymbol{{ Y }}}}){\rm d}{{\boldsymbol{{ Y }}}}. \end{equation} $

通过Parserval不等式和$ \{u_r\} $的正交性可得

(4.14) $ \begin{eqnarray} \|w(T)\|^2&=&\int_{\Omega}|w({{\boldsymbol{{ X }}}}, T)|^2{\rm d}{{\boldsymbol{{ X }}}}=\sum\limits_{r=1}^{\infty}e^{2\lambda_ra^2T}w_r^2{}\\ &=&\sum\limits_{r=1}^{\infty}e^{2\lambda_ra^2T}\left(\sum\limits_{s=1}^{\infty}G_se^{-\lambda_sa^2t_0}\int_{\Omega}u_s({{\boldsymbol{{ Y }}}}) u_r({{\boldsymbol{{ Y }}}}){\rm d}{{\boldsymbol{{ Y }}}}\right)^2{}\\ &=&\sum\limits_{r=1}^{\infty}G_r^2e^{2\lambda_ra^2(T-t_0)}\leq \sum\limits_{r=1}^{\infty}G_r^2e^{2\lambda_ra^2T}\leq M^2. \end{eqnarray} $

对于$ g({{\boldsymbol{{ X }}}})\in Q(M) $, 由于$ \log G(t'):=\log \|w(t')\|^2 $是一个凸函数(证明参见文献[18]), 故

(4.15) $ \begin{equation} \beta\log G(t_1')+(1-\beta)\log G(t_2')\geq \log G(\beta t_1'+(1-\beta)t_2'), \end{equation} $

$ \forall \beta \in(0, 1) $, $ t_1', t_2'\in[0, T] $, 取$ \beta=1-t_0/T $, $ t_1'=0 $, 和$ t_2'=T $, 带入(4.15)式, 可得

(4.16) $ \begin{equation} G(t_0)\leq G(T)^{t_0/T}G(0)^{1-t_0/T}. \end{equation} $

由于$ g({{\boldsymbol{{ X }}}})=w({{\boldsymbol{{ X }}}}, t_0) $和$ w({{\boldsymbol{{ X }}}}, 0)=u({{\boldsymbol{{ X }}}}, t_0) $, 定理4.1得证.

由定理4.1, 我们得到了反问题的条件稳定性.

定理4.2  对于$ g_i({{\boldsymbol{{ X }}}})\in Q(M) $, 我们假设$ u_i({{\boldsymbol{{ X }}}}, t_0)={{\bf K}}_1^{{\bf T}}g_i({{\boldsymbol{{ X }}}}) $, $ i=1, 2 $有

(4.17) $ \begin{equation} \|g_1-g_2\|_{L^2(\Omega)}\leq (2M)^{t_0/T}\|{{\bf K}}_1^{{\bf T}}g_1-{{\bf K}}_1^{{\bf T}}g_2\|_{L^2(\Omega)}^{1-t_0/T}, \quad \forall t_0\in(0, T) \end{equation} $

成立.

证  容易验证, 当$ g_i({{\boldsymbol{{ X }}}})\in Q(M) $, $ i=1, 2 $时, 有$ (g_1({{\boldsymbol{{ X }}}})-g_2({{\boldsymbol{{ X }}}}))\in Q(2M) $, 根据定理4.1和算子$ {\bf K}_1^{\bf T} $线性性可推出(4.17)式.

本节主要采用两种正则化方法对第4节的逆时热传导问题进行数值求解.

方法一  Tikhonov正则化方法^[19].

公式(4.3)建立了不同时刻温度场和$ T $时刻温度场之间的关系. 通过求解方程(4.3)计算$ t=0 $时的温度场, 即为初始温度场逆时问题; 计算$ t\in(0, T) $的温度场, 则对应其它时刻的逆时问题. 对于给定的扰动数据$ f^{\delta} $, 无论扰动$ \delta $有多小, 逆时问题的解也可能不存在, 因此$ u(x, y, z, t) $的近似值$ u^{\delta}(x, y, z, t) $不能直接由(4.3)式右端换为$ f^{\delta} $的方程求解. 故使用正则技术求解如下正则化方程

(5.1) $ \begin{equation} \left(\alpha I+({\bf K}_1^{{\bf T-t}})^*{\bf K}_1^{{\bf T-t}}\right)u^{\alpha, \delta}(x, y, z, t)=({\bf K}_1^{{\bf T-t}})^*f^{\delta}, \end{equation} $

其中$ {\bf K}^* $为$ {\bf K} $的伴随算子, 对于给定的$ \alpha >0 $正则化方程有唯一解$ u^{\alpha, \delta}(x, y, z) $, $ t\in[0, T) $. 当$ \alpha, \delta\rightarrow 0 $时, $ u^{\alpha, \delta}(x, y, z, t)\rightarrow u(x, y, z, t) $, 正则化参数$ \alpha=\alpha(\delta) $的选取方法见文献[17].

方法二  终值数据扰动的正则化方法.

对于给定的$ \alpha>0 $, 利用(4.1)构造对应的扰动问题如下

(5.2) $ \begin{equation} \left\{\begin{array}{ll} u_t^\alpha=a^2 (u_{xx}^\alpha+u_{yy}^\alpha+u_{zz}^\alpha), &(x, y, z, t)\in\Omega\times(0, T), \\ u^\alpha|_{\partial\Omega} =0, & 0<t<T, \\ \alpha u^\alpha(x, y, z, 0)+u^\alpha(x, y, z, T)=f(x, y, z), & (x, y, z) \in \Omega. \end{array}\right. \end{equation} $

对精确的输入数据来讨论近似解$ u^{\alpha}(x, y, z, t) $和精确解$ u(x, y, z, t) $的误差, 对给定的$ u^{\alpha}(x, y, z, 0) $满足上述模型前两式$ u^{\alpha}(x, y, z) $的表达式为

(5.3) $ \begin{eqnarray} u^\alpha(x, y, z, t)&=&\sum\limits_{l=1}^{\infty}\sum\limits_{m=1}^{\infty}\sum\limits_{n=1}^{\infty}\langle u^{\alpha}(x, y, z, 0), \psi_{lmn}\rangle e^{-a^2(l^2+m^2+n^2)t}\cdot \psi_{lmn}(x, y, z) {}\\ &:=&{\bf K}_1^{\bf t}u^\alpha(x, y, z, 0), \end{eqnarray} $

其中$ \psi_{lmn}(x, y, z)=\sqrt{\frac{8}{\pi^3}}\sin lx\sin my \sin nz $, $ t\in[0, T] $, $ l, m, n=1, 2\cdots $, $ \langle \cdot, \cdot \rangle $表示$ L^2(\Omega) $上的内积, 有界线性算子$ {\bf K}_1^{\bf t} $由(4.3)–(4.4)式定义, 将(5.3)式带入(5.2)式中的变形初始条件得

(5.4) $ \begin{equation} u^{\alpha}(x, y, z, 0)=(\alpha I+{\bf K}_1^{\bf T})^{-1}f. \end{equation} $

由(5.3)式得

$ u^{\alpha}(x, y, z, t)={\bf K}_1^{\bf t}(\alpha I+{\bf K}_1^{\bf T})^{-1}f, $

注意到

$ {\bf K}_1^{\bf t}(\alpha I+{\bf K}_1^{\bf T})^{-1}=(\alpha I+{\bf K}_1^{\bf T})^{-1}{\bf K}_1^{\bf t}, $

得$ u^{\alpha}(x, y, z, t) $满足

(5.5) $ \begin{equation} \alpha u^{\alpha}(x, y, z, t)+\left( {\bf K}_1^{\bf T}u^{\alpha}(\cdot, t)\right)(x, y, z) =({\bf K}_1^{\bf t}f)(x, y, z), \quad t\in[0, T], \end{equation} $

这是一个关于$ u^{\alpha}(x, y, z) $的第二类积分方程, 由算子$ {\bf K}_1^{\bf t} $的表达式可得

(5.6) $ \begin{equation} \lim\limits_{\alpha\to 0}\|u^{\alpha}(\cdot, t)-u(\cdot, t) \|=0, \quad t\in[0, T], \end{equation} $

对于右端项加上一定的扰动, 可进一步估计收敛速度^[17].

先考虑Tikhonov正则化方法. 给定初始温度$ g(x, y, z)=\sin x\sin y\sin z $, 以及参数$ a^2=0.04 $, 用分离变量法求解正问题可得$ u(x, y, z, t)=\sin x\sin y\sin z e^{-3a^2 t}. $利用$ T=6 $时刻的温度场$ u(x, y, z, T) $的扰动形式为

$ u^{\delta}(x, y, z, T)=u(x, y, z, T)\cdot (1+\delta rand) $

作为反演的输入数据, 反演出初始温度场$ u^{\alpha}(x, y, z, 0) $. 令$ u(x_i, y_j, z_k, 0) $表示初始温度的精确解, 分别取$ \delta=0.001, 0.005, 0.01, 0.05 $, 空间步长为$ h_1=h_2=h_3=\pi/20 $, 时间步长$ \Delta t=0.1 $, 所得反演的数值解与精确解的误差见表 2和图 2. 其中

$ e_3(x_i, y_j, z_k, 0)=u^{\alpha, \delta}(x_i, y_j, z_k, 0)-u(x_i, y_j, z_k, 0), $

$ \|e_4(x_i, y_j, z_k, 0)\|:=\frac{\|e_3(x_i, y_j, z_k, 0)\|_{L^2}}{\|u(x_i, y_j, z_k, 0)\|_{L^2}}, \\ $

$ error3(x_i, y_j, z_k, 0):=|e_3(x_i, y_j, z_k, 0)|, $

$ error4(x_i, y_j, z_k, 0):=\frac{|e_3(x_i, y_j, z_k, 0)|}{|u(x_i, y_j, z_k, 0)|}. $

通过表 2和图 2, 我们可以看出, 当正则化参数$ \alpha $取适当的值时, 数值解与精确解的误差较小, 精度达到$ 10^{-3} $, 实现了逆时反演, 同时也显示了反演算法的有效性.

下面我们考虑终值数据扰动的正则化方法. 给定初始温度$ g(x, y, z)=\sin x\sin y\sin z $, 以及参数$ a^2=0.04 $, 用分离变量法求解正问题可得$ u(x, y, z, t)=\sin x\sin y\sin z e^{-3a^2 t} $. 取不同时刻$ (T>1) $的数值解按照一定噪声水平生成反演的输入数据, 扰动形式为

$ u^{\delta}(x, y, z, T)=u(x, y, z, T)\cdot (1+\delta rand), $

反演$ t=1 $时刻的温度$ u^{\alpha}(x, y, z, t) $. 这里$ u(x_i, y_j, z_k, 1) $表示$ t=1 $时的精确解, 空间网格划分步长为$ h_1=h_2=h_3=\pi/20 $, 正则化参数分别选为$ \alpha=0.0001, 0.0005, 0.001, 0.005 $, 取$ \delta=0 $和$ \delta=0.001 $, 时间步长$ \Delta t=0.1 $, 所得反演的数值解与精确解的误差比较见表 3, 其中

$ e_5(x_i, y_j, z_k, 1)=u^{\alpha}(x_i, y_j, z_k, 1)-u(x_i, y_j, z_k, 1), $

$ \|e_6(x_i, y_j, z_k, 1)\|:=\frac{\|e_5(x_i, y_j, z_k, 1)\|_{L^2}}{\|u(x_i, y_j, z_k, 1)\|_{L^2}}, $

$ error5(x_i, y_j, z_k, 1):=|e_5(x_i, y_j, z_k, 1)|, $

$ error6(x_i, y_j, z_k, 1):=\frac{|e_5(x_i, y_j, z_k, 1)|}{|u(x_i, y_j, z_k, 1)|}. $

表 3和图 3给出噪声水平$ \delta=0 $时的反演结果, 表 4和图 4给出噪声水平$ \delta=0.001 $时的反演结果. 通过反演结果, 我们可以看出在正则化等其它参数相同的情况下, $ T $相对越小(即离$ t $相对越近)反演$ t $时刻的温度场分布误差越小. 实际上, 对于这类热传导问题, 当时间足够长时, 热传导正问题的解衰减到很小, 此时利用观测数据很难实现初始温度场的反演.

本文研究了在数学物理和工程领域中有重要应用的一类三维逆时热传导问题. 基于变分法和有限元对正问题进行了数值求解. 进一步在规则区域内使用分离变量法构造了两个反演算法, 并在一定假设条件下证明了反问题的局部稳定性. 通过数值实验验证了算法的有效性, 并分析了不同参数对反演结果的影响. 本文得到的理论和算法丰富了目前关于三维逆时热传导问题的研究结果. 文中的反演算法针对第一类和第二类边界条件都是可以适用的, 但对于一般的第三类边界条件, 由于不能直接建立分离变量形式的反演公式, 需要尝试使用其它数学理论和方法加以解决. 针对第三类边界条件的三维逆时问题, 我们将在后续的研究中继续考虑其稳定高效的反演方法.