重构导热方程源项的反问题

Reversing Inverse Problem of Source Term of Heat Conduction Equation

摘要: 本文将终端观测值应用在依赖热源的导热方程的反问题中。该问题在应用科学领域具有重要的应用。本文基于最优控制框架的基础上，建立了控制函数最小化的存在性和必要条件。应用Landweber迭代算法在该问题中，并得到了数值模拟结果。

Abstract: In the paper, terminal observations are applied to the inverse problem of heat-transfer equations that rely on heat sources. This problem has important applications in the field of applied science. Based on the optimal control framework, the existence and necessary conditions of minimization of control functions are established. The Landweber iterative algorithm is applied in this problem, and the numerical simulation results are obtained.

1. 引言

偏微分方程 [1] 是纯粹数学和应用数学的一个重要分支，它是描述与刻画物理过程、系统状态、社会与生物现象的有力工具。由给定的方程、定解区域以及相应的初、边值条件确定方程的解，这就是所谓的偏微分方程的正问题，而偏微分方程的反问题 [2] 是指已知或部分已知方程的解反求方程中的未知量。与正问题相比，大多数的反问题都是不适定 [3] 的。近年来，反问题已经成为应用数学领域迅速发展的一门理论，在医疗、油气勘探、信号探测等方面都有重要的应用。此类问题的一个主要特点是在部分边界上可能缺失边界条件。与此同时，系数的退化性也会导致方程的解没有足够的正则性。

近年来，抛物型方程系数反问题的研究得到了一系列重要的成果。例如文献 [4] 讨论了利用边界上给出的附加数据确定如下热传导方程的源项系数 $f\left(x\right)$ 的反问题：

$u_t-\Delta u=f\left(x\right)x,\left(x,t\right)\in Q.$

例如文献 [5] 利用Tikhonov正则化方法研究了Black-Scholes方程 [6]

$C_t+\frac{1}{2}{\sigma }^2{S}^2{C}_{SS}+\left(r-q\right)S{C}_S-rC=0,$

波动率的稳定性和收敛性，并对理论分析的结果进行了数值模拟。

文献 [7] 研究了利用附加条件重构热源 $p\left(t\right)$ 的反问题，得到了p的唯一性和稳定性的数值解。抛物型方程源项系数的确定，无论是在科学领域还是工程领域都有着广泛的应用。在弹性力学、流体力学 [8] 、温度场的测量 [9] [10] 、热传导和系统耦合性研究 [11] 等都会出现这些反问题。本文主要考虑如下问题：

问题P本文讨论的是二阶导热方程源项的反演问题

$\{\begin{array}{ll}{u}_t-a\left(x\right)u_x{}_x+b\left(x\right)u=f\left(x\right),\hfill & \left(x,t\right)\in Q=\left(0,l\right)\times \left(0,T\right],\hfill \\ {u|}_{x=0}={u|}_{x=l}=0,\hfill & t\in \left(0,T\right],\hfill \\ {u|}_{t=0}=\phi \left(x\right),\hfill & x\in \left(0,l\right)，\hfill \end{array}$ (1)

这里 $a\left(x\right),b\left(x\right)$ 和 $\phi$ 在 $\left(0,l\right)$ 上是给定光滑函数， $f\left(x\right)$ 在方程(1)中是一个未知的右端项。

给定如下附加条件

${u|}_{t=T}=g\left(x\right),x\in \left[0,l\right],$

其中 $g\left(x\right)$ 是一个已知的函数，它满足齐次Dirichlet边界条件。我们需要由此同时确定未知函数对 $\left(u,f\right)$ 。

2. 控制问题

本文中，我们假设 $a\left(x\right),b\left(x\right)\in C^{\alpha }\left(0,l\right)$ ， $a\left(x\right)\ge a_0>0$ ， $b\left(x\right)>0$ ， $\phi \left(x\right)\in C^{2,\alpha }\left(0,l\right)$ ， $0<\alpha <1$ ，且 $a_0$ 是一个常数， $\phi \left(x\right)$ 满足齐次Dirichlet边界条件。

由Schauder抛物方程理论我们知道，对于方程(1)任意给定的函数 $f\left(x\right)\in C^{\alpha }\left(0,l\right)$ ，存在一个方程的唯一解 $u\left(x,t\right)\in C^{2+\alpha ,1+\left(\alpha /2\right)}\left(0,l\right)$ ，我们转而考虑以下最优控制问题P1：

寻找 $\tilde{f}\left(x\right)\in {\rm A}$ ，使得：

$J\left(\tilde{f}\right)=\underset{a\in {\rm A}}{\min }J\left(f\right),$ (2)

这里

$J\left(f\right)=\frac{1}{2}{\displaystyle {\int }_0^l{\left|u\left(x,T;f\right)-g\left(x\right)\right|}^2\text{d}x}+\frac{N}{2}{\displaystyle {\int }_0^l{\left|\nabla f\right|}^2\text{d}x},$ (3)

${\rm A}=\left\{f\left(x\right)|\left|f\right|\le M,f\in H^1\left(0,l\right)\right\},$ (4)

$u\left(x,t;f\right)$ 是方程(1)对应于给定系数 $f\left(x\right)\in {\rm A}$ 的解，N是正则化参数，M是一个给定的常数。

假设终端观测数据 $g\left(x\right)$ 满足以下条件

$g\left(x\right)\in L^2\left(0,l\right).$ (5)

因此，由Schauder抛物方程理论，我们对问题(1)有以下的估计。

引理2.1令 $u\left(x,t\right)$ 为(1)的解， $f\left(x\right)\in {\rm A}$ 是给定函数，对于有以下估计

$\underset{\bar{Q}}{\max }\left|u\right|{\displaystyle {\int }_0^l{u}^2\text{d}x\text{d}t}+{\displaystyle {\int }_0^T{\displaystyle {\int }_0^l\left({\left|u_t\right|}^2+{\left|u_{xx}\right|}^2\right)\text{d}x\text{d}t}}\le C$ ，其中C为一个常数。

引理2.1证明是明显的。

从引理2.1和(5)可以看出，控制泛函(3)对 $\forall f\left(x\right)\in {\rm A}$ 是有意义的。

定理2.1存在一个的极小元 $\tilde{f}\in {\rm A}$ ，即

证明可以看出 $J\left(f\right)$ 是非负的。因此 $J\left(f\right)$ 有最大下界 $\underset{f\in {\rm A}}{\inf }J\left(f\right)$ ，令 $\left(u_n,f_n\right)$ 是一个极小化序列

$\underset{f\in {\rm A}}{\inf }J\left(f\right)\le J\left(f_n\right)\le \underset{f\in {\rm A}}{\inf }J\left(f\right)+\frac{1}{n},$

由于 $J\left(f_n\right)\le C$ 和J的特殊结构，我们推导 ${\Vert \nabla f_n\Vert }_{L^{\frac{1}{2}}\left(0,l\right)}\le C$ (C与n无关)，利用Sobolev嵌入定理，我们得到 ${\Vert f_n\Vert }_{C^{\frac{1}{2}}\left(0,l\right)}\le C$ 。

因此，我们可以选择一个和 $u_n$ 的子序列，同样用 $f_n$ 和 $u_n$ 表示，这样

$f_n\left(x\right)\to \tilde{f}\left(x\right)\in C^{\frac{1}{2}}\left(0,l\right)$ ，得到 $C^{\alpha }\left(0,l\right)\left(0\le \alpha \le \frac{1}{2}\right)$ ，

$u_n\left(x,t\right)\to \tilde{u}\left(x,t\right)$ ，得到 $C^{\alpha ,\frac{\alpha }{2}}\left(\bar{Q}\right)\cap C_{loc}^{2+2,1+\frac{\alpha }{2}}\left(Q\right)$ ，

很容易验证 $\left(\tilde{f}\left(x\right),\tilde{u}\left(x,t\right)\right)$ 满足(1.1)，利用Lebesgue控制收敛定理和 $L^2$ 上的弱半连续定理，我们可以获得

$J\left(\tilde{f}\right)\le \underset{n\to \infty }{\lim }\inf J\left(f_n\right)=\underset{f\in {\rm A}}{\min }J\left(f\right),$

因此 $J\left(\tilde{f}\right)=\underset{f\in {\rm A}}{\min }J\left(f\right)$ 。

定理2.1证完。

3. 必要条件

定理3.1若f是最优控制问题(2)的解，则存在一个三元函数组 $\left(u,v;f\right)$ 满足下列系统

(6)

$\{\begin{array}{ll}{v}_t-a\left(x\right)v_x{}_x+bv=h-f,\hfill & \left(x,t\right)\in Q,\hfill \\ v\left(0,t\right)=v\left(l,t\right)=0,\hfill & t\in \left(0,T\right],\hfill \\ v\left(x,0\right)=0,\hfill & x\in \left(0,l\right),\hfill \end{array}$ (7)

且

$N{\displaystyle {\int }_0^l\nabla f\cdot \nabla \left(h-f\right)\text{d}x}+{\displaystyle {\int }_0^{l}v\left(h-f\right)\text{d}x}\ge 0.$ (8)

证明对于 $\forall h\in {\rm A}$ ， $\text{0}\le \delta \le \text{1}$ ，令 $f_{\delta }\equiv \left(1-\delta \right)f+\delta h\in {\rm A}$ ，则

$J_{\delta }\equiv J\left(f_{\delta }\right)=\frac{1}{2}{\displaystyle {\int }_0^l{\left|u\left(x,T;f_{\delta }\right)-g\left(x\right)\right|}^2\text{d}x}+\frac{N}{2}{\displaystyle {\int }_0^l{\left|\nabla f_{\delta }\right|}^2\text{d}x}.$ (9)

令 $u_{\delta }$ 是方程(1)的解，其中 $f=f_{\delta }$ ，我们有

${\frac{\text{d}{J}_{\delta }}{\text{d}\delta }|}_{\delta =0}={{\displaystyle {\int }_0^l\left[u\left(x,T;f\right)-g\left(x\right)\right]\frac{\partial u_{\delta }}{\partial \delta }|}}_{\delta =0}\text{d}x+N{\displaystyle {\int }_0^l\nabla f\cdot \nabla \left(h-f\right)\text{d}x}\ge 0.$ (10)

令 ${\tilde{u}}_{\delta }\equiv \frac{\partial u_{\delta }}{\partial \delta }$ ，直接计算得到以下方程

$\{\begin{array}{l}{\tilde{u}}_{\delta ,t}-a{\tilde{u}}_{\delta ,xx}+b{\tilde{u}}_{\delta }=h-f,\\ {{\tilde{u}}_{\delta }|}_{x=0}=0,{{\tilde{u}}_{\delta }|}_{x=l}=0,\\ {\tilde{u}}_{\delta }\left(x,0\right)=0.\end{array}$ (11)

令 $\xi ={{\tilde{u}}_{\delta }|}_{\delta =0}$ ，则 $\xi$ 满足

$\{\begin{array}{l}{\xi }_t-a{\xi }_x{}_x+b\xi =h-f,\hfill \\ \xi \left(0,t\right)=\xi \left(l,t\right)=0,\hfill \\ \xi \left(x,0\right)=0.\hfill \end{array}$ (12)

从(10)中，我们得到

${\displaystyle {\int }_0^l\left[u\left(x,T;f\right)-g\left(x\right)\right]\xi \left(x,T\right)\text{d}x}+N{\displaystyle {\int }_0^l\nabla f\cdot \nabla \left(h-f\right)\text{d}x}\ge 0.$ (13)

记 $L\xi ={\xi }_t-a{\xi }_x{}_x-b\xi$ ，令v是下列问题的解

$\{\begin{array}{l}{L}^{\ast }v=-v_t-a{v}_x{}_x+bv=0,\hfill \\ v\left(x,T\right)=u\left(x,T\right)-g\left(x\right),\hfill \\ v\left(0,t\right)=v\left(l,t\right)=0.\hfill \end{array}$ (14)

其中算子 $L^{\ast }$ 是L的共轭算子，于是利用Green公式和分部积分容易得到：

$\begin{array}{c}0={\displaystyle {\int }_0^T{\displaystyle {\int }_0^l\xi L^{\ast }v\text{d}x\text{d}t}}=-{\displaystyle {\int }_0^l\xi \left(x,T\right)\left[u\left(x,T\right)-g\left(x\right)\right]\text{d}x}+{\displaystyle {\int }_0^T{\displaystyle {\int }_0^{l}vL\xi \text{d}x\text{d}t}}\\ =-{\displaystyle {\int }_0^l\xi \left(x,T\right)\left[u\left(x,T\right)-g\left(x\right)\right]\text{d}x}+{\displaystyle {\int }_0^T{\displaystyle {\int }_0^{l}v\left(h-f\right)\text{d}x\text{d}t}}\end{array}$ (15)

由(13)和(15)可得到

$N{\displaystyle {\int }_0^l\nabla f\cdot \nabla \left(h-f\right)\text{d}x}+{\displaystyle {\int }_0^{l}v\left(h-f\right)\text{d}x}\ge 0.$

定理3.1证完。

4. 迭代格式步骤

第一步：选择一个初始的迭代函数初始函数 $f=f_0\left(x\right)$ ，可以任意的选取，本篇论文中我们选择

$f_0\left(x\right)=x\left(l-x\right),x\in \left(0,l\right)$ ;

第二步：通过解初边值问题(1)得到 $u_0\left(x,t\right)$ ，其中 $f=f_0\left(x\right)$ ；

第三步：解共轭方程(14)得到 $v_0\left(x,t\right)$ ，其中 $u\left(x,T\right)=u_0\left(x,T\right)$ ；

第四步：令 $f_1=f_0\left(x\right)-\alpha v_0\left(x,T\right)$ ，其中 $\alpha \ge 0$ ，并且令 $u_1\left(x,T\right)$ 是(1)的解，此时 $f=f_1\left(x\right)$ ；

第五步：选择一个任意小的正常数 $\epsilon$ 作为误差限，计算 $\Vert u_1\left(x,T\right)-g\left(x\right)\Vert$ 并且和 $\epsilon$ 比较大小。如果：

$\Vert u_1\left(x,T\right)-g\left(x\right)\Vert \le \epsilon$ ，

终止迭代，这时取 $f=f_1\left(x\right)$ ；如果：

$\Vert u_1\left(x,T\right)-g\left(x\right)\Vert \ge \epsilon$ ，

那么继续执行第三步，并且令 $f_1\left(x\right)$ 为新的迭代初值继续执行归纳的准则直到迭代满足终止条件。

一般来说，如果输入数据是精确的，Landweber迭代法的迭代次数越多，那么输出的数据精度也越高，对于有扰动的情况，在初始的迭代过程中就存在误差，这个计算误差一开始会随着迭代次数的增加而减少。因此，必须在适当的时刻终止迭代。所以，要选择合适的参数，使得迭代格式既是精度高又有好的稳定性的方法。

5. 数值实验

在这个数值实验中，取 $a\left(x\right)=1,b=0$ ，较正参数 $\alpha =1$ ，源项 $f\left(x\right)=\left(2{\text{π}}^2-{\text{π}}^4\right)\cdot \sin \text{π}x$ ，初值条件

$\phi \left(x\right)=\left(1-{\text{π}}^2\right)\sin \text{π}x,x\in \left[0,1\right].$

正问题的解析解为 $u=\left(2-{\text{π}}^2-{\text{e}}^{-{\text{π}}^{2}t}\right)\sin \text{π}x,\left(x,t\right)\in \left[0,1\right]\times \left[0,1\right]$ 。因此，终端观测值为

$g\left(x\right)=\left(2-{\text{π}}^2-{\text{e}}^{-{\text{π}}^2}\right)\sin \text{π}x,x\in \left[0,1\right].$

重构结果如图1所示，我们可以看到，当迭代次数为500次时，源项能够被很好的反演

Figure 1. Source item numerical inversion

图1. 源项数值反演

下面我们来考虑加噪的情况，给出如下的加噪格式

$g^{\delta }=u^{\delta }\left(x,1\right)=u\left(x,1\right)\left[1+\delta \times random\left(x\right)\right]$ ，

在图2中我们可以看到，当噪声水平在 $\delta =5\%$ 和 $\delta =7\%$ 时，未知源项也可以被很好的重构。

Figure 2. Numerical value inversion of noisy source term

图2. 加噪的源项数值反演

备注应当指出，在不适定问题的数值模拟中，正则化参数N扮演者一个非常重要的角色。