热传导方程作为三大偏微分方程之一, 它广泛的应用于物理化学、航天航空[1]、国防军事[2]及机械制造[3]等各个领域. 通常热传导方程分为直接热传导问题(Direct heat conduction problems(DHCPs), 简称正问题)和反热传导问题(Inverse heat conduction problems(IHCPs), 简称反问题). 热传导方程反问题具体可分为三大类：(1) 初始条件反问题. 即通过已知条件确定初始时刻的温度分布; (2) 侧边反问题[4]. 即通过已知条件确定某一边界上的温度分布; (3) 源项反问题. 即通过已知条件确定源项的温度分布. 自反问题出现以来, 人们对于热传导方程反问题的研究从未间断, Djerrar[5]等人在基于Laplace变换求解正问题的基础上, 得到侧边反问题的模型, 并利用Tikhonov正则化进行反演求解; M.Nili Ahmadabadi[6]等人讨论了热源性仅与空间相关的逆热源问题; Nakamura et al.[7]通过变换技术解决了初始条件逆问题; 肖庭延[8]等人对研究反问题经典的正则化方法进行了改进, 提出了迭代正则化方法; 葛美宝[9]等人讨论了利用中心差分法求解一类含源项的热传导方程反问题; 任常青[10]利用拟逆法和积分方程方法对热传导方程的反问题进行了研究; 刘继军[11]也对拟逆法求解热传导逆时问题做出了详细的阐述.

综观国内外对反问题的研究进展, 可以看出人们对一维热传导方程反问题的研究较为深入, 但是对高维热传导方程的反问题研究较少. 本文将给出一种反演二维热传导方程初始条件反问题的数值求解方法：通过Crank-Nicolson-Galerkin有限元法对二维热传导方程进行离散处理, 并利用GMRES算法的改进算法–RRGMRES进行求解. 最后对本文提出的算法进行了具体的数值模拟.

考虑如下的二维热传导方程：

$ \begin{equation} \left\{ \begin{array}{l} \frac{{\partial u}}{{\partial t}} = \frac{{{\partial ^2}u}}{{\partial {x^2}}} + \frac{{{\partial ^2}u}}{{\partial {y^2}}}\;\;\quad0 < x < p, 0 < y < q, t > 0\\ {\left. {\frac{{\partial u}}{{\partial x}}} \right|_{x = 0}} = {\left. {\frac{{\partial u}}{{\partial x}}} \right|_{x = p}} = 0\;\;\quad0 < y < q, t > 0\\ {\left. {\frac{{\partial u}}{{\partial y}}} \right|_{y = 0}} = {\left. {\frac{{\partial u}}{{\partial y}}} \right|_{y = q}} = 0\;\;\quad0 < x < p, t > 0\\ u(x, y, 0) = f(x, y)\;\;\quad0 < x < p, 0 < y < q \end{array} \right. \end{equation} $

(2.1)

其中$ u(x, y, t) $为温度分布, $ f(x, y) $为初始条件, $ x, y $为空间变量, $ t $为时间变量. 当初始条件$ f(x, y) $已知, 需要求出任意$ t $时刻的温度分布$ u(x, y, t) $, 即为热传导方程的正问题.

若已知

$ \begin{equation} u\left( {x, y, T} \right) = g\left( {x, y} \right), \end{equation} $

(2.2)

其中$ T > 0 $, 反求$ t=0 $时的$ u(x, y, 0) $, 即求初始条件$ f(x, y) $, 这就是二维热传导方程初始条件反问题.

结合方程(2.1), 不难发现：对每一个初始条件$ f(x, y) $均存在方程(2.1)在$ T $时刻的一个解函数$ u(x, y, T) $, 即存在一个从$ f(x, y) $到$ u(x, y, T) $的非线性算子$ {A_1} $(即$ {A_1}\left[ {f(x, y)} \right] = u(x, y, T) $); 由(2.2)式可知, 存在一个从$ u(x, y, T) $到$ g(x, y) $的线性算子$ {A_2} $(即$ {A_2}\left[ {u(x, y, T)} \right] = g(x, y) $)). 因此, 求解该初始条件反问题即是求解下列的非线性算子方程[12]：

$ \begin{equation} {A_2}\left[ {{A_1}\left[ {f(x, y)} \right]} \right] = g(x, y). \end{equation} $

(2.3)

注意到对于上述(2.2)所示的附加条件, $ T $时刻的温度分布$ u(x, y, t) $通常是通过测量得到的, 必然会存在一定的误差. 因此, 反问题则可简述为通过测量数据估计未知初始条件$ f(x, y) $.

由于反问题的求解是建立在正问题高精度求解的基础上, 下面利用Crank-Nicolson-Galerkin有限元法对二维热传导方程的正问题进行求解.

假设函数$ u(x, y, t) $充分光滑, $ \Omega = [0, p] \times [0, q] $. 对方程(2.1)第一式的两边分别乘以$ H_0^1(\Omega ) $中的函数$ v(x, y) $, 再在$ \Omega $上积分：

$ \begin{equation} \iint\limits_{\Omega}u_{t}(x, y, t)v(x, y)dxdy=\iint\limits_{\Omega}\Delta u(x, y, t)v(x, y)dxdy. \end{equation} $

(3.1)

对上式的右端项利用格林公式, 并结合边界条件, 同时注意到$ v \in H_0^1(\Omega ) $, 则(3.1) 式可以转化为

$ \begin{equation} \frac{d}{{dt}}{\iint\limits_{\Omega}u(x, y, t)v(x, y)}dxdy=-{\iint\limits_{\Omega}u_{x}(x, y, t)v_{x}(x, y)+u_{y}(x, y, t)v_{y}(x, y)}dxdy. \end{equation} $

(3.2)

为便于书写, 引入如下的内积定义：

对于函数$ \varphi , \psi \in {L^2}(\Omega ) $, 定义$ {L^2}(\Omega ) $中的内积

$ \begin{equation} \left( {\varphi , \psi } \right)=\iint\limits_{\Omega}{\varphi (x, y)\psi (x, y)}dxdy. \end{equation} $

(3.3)

对$ \varphi , \psi \in H_0^1(\Omega ) $, 定义$ H_0^1(\Omega ) \times H_0^1(\Omega ) $上的双线性泛函：

$ \begin{equation} {A^ * }(\varphi , \psi )={\iint\limits_{\Omega}\varphi_{x}(x, y, t)\psi_{x}(x, y)+\varphi_{y}(x, y, t)\psi_{y}(x, y)}dxdy. \end{equation} $

(3.4)

则求解方程(2.1)的正问题可以转化为如下的Galerkin弱形式：

求$ u:t \in [0, T] \to u(t) \in H_0^1(\Omega ) $, 使得

$ \begin{equation} \left\{ \begin{array}{l} \frac{d}{{dt}}\left( {u(t), v} \right) + {A^ * }(u(t), v) = 0\;\; \quad \forall v \in H_0^1(\Omega )\\ u(0) = {u_0} = f(x, y) \end{array} \right. \end{equation} $

(3.5)

下面对问题(3.5)进行求解：

首先考虑如下的半离散问题：给定$ {u_{0, h}} \in {V_h} $, 求函数$ {u_h}:t[0, T] \to {u_h}(t) \in {V_h} $, 使得

$ \begin{equation} \left\{ \begin{array}{l} \frac{d}{{dt}}\left( {{u_h}(t), v} \right) + {A^ * }({u_h}(t), v) = 0\;\;\quad \forall v \in V\\ {u_h}(0) = {u_{0, h}} = f(x, y) \end{array} \right. \end{equation} $

(3.6)

假设其中$ {V_h} $是$ V $的有限维子空间, 其维数为$ d = d(h) $. 为便于求解, 现对方程(3.6)稍作变形. 假设空间$ {V_h} $的一组基为$ \left\{ {{\varphi _1}, {\varphi _2}, \cdots , {\varphi _d}} \right\} $, 则$ {u_h}(t) \in {V_h} $可表示为

$ \begin{equation} {u_h}(t) = \sum\limits_{j = 1}^d {{\xi _j}(t){\varphi _j}} , t \in [0, T] \end{equation} $

(3.7)

初值可表示为

$ \begin{equation} {u_{0, h}}(t) = \sum\limits_{j = 1}^d {{\xi _{0, j}}(t){\varphi _j}}. \end{equation} $

(3.8)

将(3.7)、(3.8)式代入到半离散问题(3.6)式中, 可将问题(3.6)表示为

$ \begin{equation} \left\{ \begin{array}{l} \sum\limits_{j = 1}^d {\left( {{\varphi _j}, {\varphi _i}} \right)\frac{{d{\xi _j}(t)}}{{dt}}} + \sum\limits_{j = 1}^d {{A^ * }\left( {{\varphi _j}, {\varphi _i}} \right){\xi _j}(t)} = 0\qquad \;\;1 \le i \le d\\ {\xi _{_i}}(0) = {\xi _{0, i}} \end{array} \right. \end{equation} $

(3.9)

将(3.9)式写成矩阵向量形式为

$ \begin{equation} \left\{ \begin{array}{l} M\frac{{d\xi (t)}}{{dt}} + K\xi (t) = 0\\ \xi (0) = {\xi _0} \end{array} \right. \end{equation} $

(3.10)

其中$ M = \left[ {({\varphi _j}, {\varphi _i})} \right] $, $ K = \left[ {{A^ * }({\varphi _j}, {\varphi _i})} \right] $, $ \xi (t) = {\left[ {{\xi _1}(t), {\xi _2}(t), \cdots , {\xi _d}(t)} \right]^T} $, $ {\xi _0} = {\left[ {{\xi _{0, 1}}, {\xi _{0, 2}}, \cdots , {\xi _{0, d}}} \right]^T} $. 此时, 方程(3.10)为含有$ d $个未知数的常微分方程, 现对常微分方程(3.10)利用$ \theta $方法对时间离散, 可得到如下的全离散方程：

$ \begin{equation} \left\{ \begin{array}{l} \frac{1}{{\Delta t}}M\left( {{\xi ^{n + 1}} - {\xi ^n}} \right) + K\left[ {\theta {\xi ^{n + 1}} + (1 - \theta ){\xi ^n}} \right] = 0\\ {\xi ^0} = {\xi _0} \end{array} \right. \end{equation} $

(3.11)

即

$ \begin{equation} (M + \theta \Delta tK){\xi ^{n + 1}} = (M - \theta \Delta tK){\xi ^n}\qquad\;\;n = 0, 1, 2, \cdots N - 1. \end{equation} $

(3.12)

当取$ \theta = \frac{1}{2} $时, (3.12)式称为Crank-Nicolson格式. 根据(3.12)式, 依次迭代得出$ \xi (t) $, 并将其代入到(3.7)式中, 进而得出热传导方程(2.1)正问题的数值解.

对初始条件进行反演时, 利用Crank-Nicolson-Galerkin有限元方法对二维热传导方程进行离散, 最终得到一个线性方程组$ Ax = b $, 其中$ A $是由有限元得到的系数矩阵, $ b = {\left[ {g({x_1}, {y_1}), g({x_2}, {y_1}), \cdots , g({x_n}, {y_1}), g({x_1}, {y_2}), \cdots , g({x_n}, {y_n})} \right]^T} $, $ x $为剖分网格节点上$ f(x, y) $的值, 即$ x = {\left[ {f({x_1}, {y_1}), f({x_2}, {y_1}), \cdots , f({x_n}, {y_1}), \cdots , f({x_n}, {y_n})} \right]^T} $.

限制值域的GMRES算法是对GMRES算法[13]的一种改进, 其基本思想与GMRES算法相同, 均是将离散后得到的线性方程组$ Ax = b $投影到某个Krylov子空间[14]$ {K_m} $上进行求解, 即求$ {x_m} \in {x_0} + {K_m} $, 使得

$ \begin{equation} \left\| {{r_m}} \right\| = \mathop {\min }\limits_{x \in {x_0} + {K_m}} \left\| {b - Ax} \right\|. \end{equation} $

(4.1)

限制值域的GMRES算法则是用$ {K_m} = {K_m}(A, Ab) = span\{ Ab, {A^2}b, \cdots , {A^m}b\} $替代GMRES算法[15]中的Krylov子空间. 其具体算法步骤为

Step1：选初始值$ {x_0} $, 精度$ \varepsilon $, 计算残余向量$ {r_0} = b - A{x_0} $, 以及初始正交向量$ {v_1} = \frac{{A{r_0}}}{{\left\| {A{r_0}} \right\|}} $;

Step2：Arnoldi过程：求出$ \left\{ {{v_i}} \right\}_{i = 1}^m $和$ {\bar H_m} $

2.1 正交化：$ {h_{ij}} = v_i^TA{v_j}, i = 1, 2, \cdots , j $, $ {\hat v_{j + 1}} = A{v_j} - \sum\limits_{i = 1}^j {{h_{ij}}{v_i}} , j = 1, 2, \cdots , m $, $ {h_{j + 1, j}} = {\left\| {{{\hat v}_{j + 1}}} \right\|_2} $;

2.2 标准化：若$ {h_{j + 1, j}} < \varepsilon $, 则转Step4, 否则, 记$ {v_{j + 1}} = {\hat v_{j + 1}}/{h_{j + 1, j}} $, $ j = j + 1 $;

2.3 求出新的$ {V_{j + 1}} $和$ {\bar H_j} $：$ {V_{j + 1}} = \left( {{V_j}, {v_{j + 1}}} \right) $, $ {\bar H_j} = \left[ {\begin{array}{*{20}{c}} {{{\bar H}_{j - 1}}}&{{{\bar h}_j}}\\ 0&{{h_{j + 1, j}}} \end{array}} \right] $,

其中$ {\bar h_j} = {\left[ {\begin{array}{*{20}{c}} {{h_{1j}}}&{{h_{2j}}}& \cdots &{{h_{jj}}} \end{array}} \right]^T}; $

Step3：解最小二乘问题:

$ \begin{array}{l} \left\| {{r_m}} \right\| = \mathop {\min }\limits_{{x_m} \in {x_0} + {K_m}(A, Ab)} \left\| {b - Ax} \right\| = \mathop {\min }\limits_{{z_m} \in {K_m}(A, Ab)} \left\| {{r_0} - Az} \right\| = \mathop {\min }\limits_{{y_m} \in {R^m}} \left\| {{r_0} - A{V_m}y} \right\|\\ \; \; \; \; \; \; \; \; \; \; = \mathop {\min }\limits_{{y_m} \in {R^m}} \left\| {{r_0} - {V_{m + 1}}{{\bar H}_m}y} \right\| = \mathop {\min }\limits_{{y_m} \in {R^m}} \left\| {{V_{m + 1}}(\left\| {{r_0}} \right\|{e_1} - {{\bar H}_m}y)} \right\|\\ \; \; \; \; \; \; \; \; \; \; = \mathop {\min }\limits_{{y_m} \in {R^m}} \left\| {\left\| {{r_0}} \right\|{e_1} - {{\bar H}_m}y} \right\|, \end{array} $

得到$ {y_m} $, 即有$ {x_m} = {x_0} + {z_m} = {x_0} + {V_m}{y_m} $;

Step4：计算$ \left\| {{r_m}} \right\| = \left\| {b - A{x_m}} \right\| $, 若$ \left\| {{r_m}} \right\| < \varepsilon $, 则停止迭代, 输出$ {x_m} $, 否则, 记$ {r_0} = {r_m}, {x_0} = {x_m}, {v_1} = \frac{{{r_m}}}{{\left\| {{r_m}} \right\|}} $, 转Step2.

为了评价本文方法的有效性, 数值模拟时, 定义相对误差

$ \begin{equation} {\mathop{\rm Re}\nolimits} = \frac{{\sqrt {\sum\limits_{i = 1}^n {\sum\limits_{j = 1}^m {{{\left[ {{f_T}({x_i}, {y_j}) - f({x_i}, {y_j})} \right]}^2}} } } }}{{\sqrt {\sum\limits_{i = 1}^n {\sum\limits_{j = 1}^m {{{\left[ {{f_T}({x_i}, {y_j})} \right]}^2}} } } }}, \end{equation} $

(5.1)

其中$ {f_T}({x_i}, {y_j}) $为点$ ({x_i}, {y_j}) $处的精确解, $ f({x_i}, {y_j}) $为点$ ({x_i}, {y_j}) $处的数值解.

数值算例  考虑如下的二维齐次热传导方程

$ \begin{equation} \left\{ \begin{array}{l} \frac{{\partial u}}{{\partial t}} = \frac{{{\partial ^2}u}}{{\partial {x^2}}} + \frac{{{\partial ^2}u}}{{\partial {y^2}}}\;\;\quad0 < x < 1, 0 < y < 1, t > 0\\ {\left. {\frac{{\partial u}}{{\partial x}}} \right|_{x = 0}} = {\left. {\frac{{\partial u}}{{\partial x}}} \right|_{x = 1}} = 0\;\;\quad0 < y < 1, t > 0\\ {\left. {\frac{{\partial u}}{{\partial y}}} \right|_{y = 0}} = {\left. {\frac{{\partial u}}{{\partial y}}} \right|_{y = 1}} = 0\;\;\quad0 < x < 1, t > 0\\ u(x, y, 0) = f(x, y)\;\;\quad0 < x < 1, 0 < y < 1 \end{array} \right. \end{equation} $

(5.2)

模拟时, 不妨先给出初始条件$ f(x, y) $的真解, 利用Crank-Nicolson-Galerkin有限元方法计算出$ u(x, y, T) $的值. 再利用$ u(x, y, T) $反求初始条件$ f(x, y) $, 以验证所提出方法的有效性.

本例中不妨设真解$ u(x, y, t) = {e^{ - 2{\pi ^2}t}}\cos (\pi x)\cos (\pi y) $, 则$ u(x, y, 0) = f(x, y) $, 即$ f(x, y) = \cos (\pi x)\cos (\pi y) $. 此外, 给出附加条件

$ \begin{equation} u(x, y, T) = {e^{ - 2{\pi ^2}T}}\cos (\pi x)\cos (\pi y). \end{equation} $

(5.3)

取$ n = 61 $, 即将区域$ [0, 1] \times [0, 1] $剖分成7200个小三角单元. 表 5-1和表 5-2分别给出了附加条件中$ T = 0.1 $和$ T = 0.5 $时, $ f({x_i}, {y_i}) $在$ i = 1, 2, \cdots , 10 $处的数值解与精确解的结果. 图 5-1和图 5-2分别给出$ T = 0.1 $和$ T = 0.5 $时数值解与精确解的对比. 其中$ T = 0.1 $时的相对误差为1.281711596334807e-04.

下面考虑获取数据受到噪声干扰影响的情况. 本文选取高斯噪声, 即

$ \begin{equation} {b^\delta } = b + \delta \, \text{randn(size}(b)), \end{equation} $

(5.4)

分别取噪声水平$ \delta = 0;0.01;0.001;0.005 $. 图 5-3给出噪声水平为0.001时的精确解与数值解的对比; 表 5-3给出了$ T=0.5 $时, 不同噪声水平下的相对误差.

同样, 对于方程(6), 取$ n=21 $, 即将区域$ [0, 1] \times [0, 1] $剖分成800个小的三角单元. 表 5-4给出了附加条件中$ T=0.5 $时, $ f({x_i}, {y_j}) $在$ {x_i} = 0.1i, y = 0.1i(i = 1, 2, \cdots , 10) $处的数值解与精确解的结果. 其相对误差为0.002806218328866. 图 5-4给出$ T=0.5 $时数值解与精确解的对比.

下面将讨论利用本文的算法在迭代不同次数时数值结果的情况. 表 5-5给出$ T=0.5 $时利用RRGMRES算法对初始条件进行反演时, 选取不同迭代次数的相对误差结果; 图 5-5给出了不同迭代次数$ k $与相对误差间的关系.

由表 5-1和表 5-2可看出本文的离散和求解方法对初始条件的反演较为准确, 其相对误差可达$ {10^{ - 4}} $; 同时由图 5-1、图 5-2可看出数值解与精确解图像的吻合度较高, 即使数据受到一定的噪声影响, 图 5-3也体现出图像较高的吻合度. 此外, 由表 5-1和表 5-4可明显看出在对方程用Crank-Nicolson-Galerkin有限元离散时, 将区域$ [0, 1] \times [0, 1] $剖分得越细, 也即$ n $的取值越大, 其数值解的精度越高. 值得注意的是$ n $的增大将会降低算法的运行速度. 由表 5-5也可以看出, 采用本文算法对初始条件进行反演时, 迭代次数越多, 解的相对误差就越小. 此外, 由图 5-5可看出当迭代次数大致超过120次时, 其相对误差基本趋于稳定, 且精度保持在$ {10^{ - 4}} $.

本文通过Crank-Nicolson-Galerkin有限元方法对含有Neumann边界条件的二维齐次热传导方程进行离散处理, 并给出用该方法求解其正问题的理论分析. 在此基础上, 利用限制值域的GMRES算法对其初始条件进行反演. 在数值实验中, 分别考虑了获取数据未受到噪声干扰和受到噪声干扰以及剖分不同数量三角单元的情况, 模拟结果表明该方法具有一定的可行性和有效性.