考虑二维平面圆形区域$\Omega$上的Laplace方程的边值反问题

$\left\{ {\begin{array}{*{20}{l}} { - \Delta u = 0\;{\rm{in}}\;\;\Omega ,}\\ {u = {u_0}\;{\rm{on}}\;\;{\Gamma _0},}\\ {{u_n} = {g_0}\;{\rm{on}}\;\;{\Gamma _0},} \end{array}} \right.$

(1.1)

其中$u_n=\dfrac{\partial u}{\partial n}$是$\Omega$的边界$\Gamma$上的外法向导数, $u_0, g_0$是$\Gamma_0$上的已知函数, $\Gamma=\Gamma_0 \bigcup \Gamma_1$, 并且$\Gamma_0 \bigcap \Gamma_1 = \varnothing$.当$u$和$u_n$分别满足周期性条件$\displaystyle\int_\Gamma u(s) ds=0$和相容性条件$\displaystyle\int_\Gamma u_n(x) dx=0$时, 要确定边界$\Gamma_1$上的未知信息, 即$u$和$u_n$在部分边界$\Gamma_1$上的值.

Laplace方程的边值反问题也称为Cauchy问题, 这类问题在工程技术如:地质勘探、卫星测量、空间遥感、目标识别以及医学成像等领域有着深刻的应用背景.根据已有的关于解的信息反演初值或边界条件通常是很困难的, 目前理论上还没有统一的方法, 因此偏微分方程的反问题研究具有很大的挑战.

问题1.1已经被证明了是严重不适定问题, 即初值的任意很小的改变都会引起解的巨大改变. Belgacem Belgacem^[1]对Cauchy问题的严重不适定性进行了理论分析; Calderón和Engl^{[2, 3]}研究了Cauchy问题弱解的存在性与唯一性; 文献[4]证明Cauchy问题是条件稳定的, 即在一个附加的有界条件下, 问题的解连续依赖初始数据.关于Laplace方程边值反问题的数值方法目前已经有了一些成果^[5-7], 其中Li^[7]利用圆周上Hilbert变换的特殊性质构造了Laplace方程边值反问题的积分方程解法.本文主要研究奇异积分在偏微分方程边值反问题中的应用, 构造求解这类问题的一种新的数值方法, 并用数值实验验证方法的有效性.

设$u_0=u|_{\Gamma_0}$和$g_0=u_n|{\Gamma_0}$是部分边界$\Gamma_0$上的已知函数, $u_1=u|_{\Gamma_1}$和$g_1=u_n|{\Gamma_1}$是部分边界$\Gamma_1$上的待定函数, $u$和$u_n$分别满足周期性条件$\displaystyle\int_\Gamma u(s) ds=0$和相容性条件$\displaystyle\int_\Gamma u_n(x) dx=0$.根据自然边界归化原理^[8]可知, 单位圆内区域$\Omega$上的调和方程满足自然积分方程

${u_n}(x) = - \frac{1}{{4\pi }} = \int_0^{2\pi } {\frac{{u(s)}}{{{{\sin }^2}\frac{{x - s}}{2}}}ds,\;\;\;x \in [0,2\pi ],} $

(2.1)

其反演公式

$u(s) = - \frac{1}{\pi }\int_0^{2\pi } {\ln } \left| {2\sin \frac{{x - s}}{2}} \right|{u_n}(x)dx,\;\;\;s \in [0,2\pi ]$

(2.2)

是一个弱奇异积分方程.由方程(2.1) 和(2.2) 可得

$ - \frac{1}{{4\pi }}\int_{{\Gamma _1}} {\frac{{{u_1}(s)}}{{{{\sin }^2}\frac{{x - s}}{2}}}} ds = {g_0}(x) + \frac{1}{{4\pi }} = \int_{{\Gamma _0}} {\frac{{{u_0}(s)}}{{{{\sin }^2}\frac{{x - s}}{2}}}ds} ,\;\;\;x \in {\Gamma _0}$

(2.3)

和

$-\dfrac{1}{\pi} \int_{\Gamma_1} \ln \left| 2 \sin \dfrac{x-s}{2} \right |g_1(x) dx=u_0(s)+\dfrac{1}{\pi} \int_{\Gamma_0} \ln \left| 2 \sin \dfrac{x-s}{2} \right |g_0(x) dx,\ \ \ s\in \Gamma_0.$

(2.4)

考虑超奇异积分

${\cal K}\varphi (\beta ) = - \frac{1}{{4\pi }} = \int_0^{2\pi } {\frac{{\varphi (\theta )}}{{{{\sin }^2}\frac{{\theta - \beta }}{2}}}d\theta } $

(2.5)

和弱奇异积分

$\mathcal{N} \varphi(\beta)= -\frac{1}{\pi}\int_0^{2\pi} \ln \left[ 2\sin \frac{\theta-\beta}{2}\right] \varphi(\theta) d\theta,$

(2.6)

其中$\varphi(\theta)$是以$2\pi$为周期的函数.对于任意的周期函数$\varphi(\theta)$, 其$n$阶三角多项式插值可定义为

$\varphi_n(\theta)=\frac{1}{2n+1}\sum\limits_{k=0}^{2n} \varphi(\theta_k)[1+2\sum\limits_{m=1}^{n} \cos(\theta-\theta_k)],$

(2.7)

用$\varphi_n(\theta)$分别替换(2.5) 和(2.6) 式中的$\varphi(\theta)$, 得到超奇异和弱奇异积分的计算公式分别为

${{\cal K}_n}\varphi (\beta ) = - \frac{1}{{4\pi }} = \int_0^{2\pi } {\frac{{{\varphi _n}(\theta )}}{{{{\sin }^2}\frac{{\theta - \beta }}{2}}}d\theta } = \sum\limits_{k = 0}^{2n} {{\omega _k}} (\beta )\varphi ({\theta _k})$

(2.8)

和

$\mathcal{N}_n \varphi(\beta)= -\frac{1}{\pi}\int_0^{2\pi} \ln \left[2\sin \frac{\theta-\beta}{2}\right] \varphi_n(\theta) d\theta=\sum\limits_{k=0}^{2n} \omega_k ' (\beta)\varphi(\theta_k), $

(2.9)

其中

$\omega_k (\beta)=\frac{2}{2n+1}\sum\limits_{m=1}^{n} \frac 1m \cos m(\beta-\theta_k),$

(2.10)

$\omega_k '(\beta)= \frac{2}{2n+1}\sum\limits_{m=1}^{n} m\cos m(\beta-\theta_k).$

(2.11)

当$\varphi(\theta)$足够光滑时, 公式(2.8) 和(2.9) 是非常精确的, 有关它们的详细讨论可以参见文献[9].

考虑超奇异积分方程

$ - \frac{1}{{4\pi }}{\rm{ = }}\int_0^{2\pi } {\frac{{\varphi (\theta )}}{{{{\sin }^2}\frac{{\theta - \beta }}{2}}}d\theta } = f(\beta ),\;\;\beta \in [0,2\pi ],$

(2.12)

并假设它满足相容性条件$\displaystyle\int_0^{2\pi} f(\beta) d\beta=0$和周期性条件$\displaystyle\int_0^{2\pi} \varphi(\theta) d\theta=0$.利用三角插值积分公式(2.8) 逼近超奇异积分方程(2.12), 取配置点$\beta_i(0\le i \le 2n)$并排列所得到的方程, 由此可以得到如下线性方程组

$\sum\limits_{j=0}^{2n} \omega_j(\beta_i) \varphi_j = f(\beta_i), \ \ \ 0\le i \le 2n,$

(2.13)

或写成矩阵的形式

$\mathbb{A} \Phi = F,$

(2.14)

其中

$\mathbb{A}=(a_{ij})_{2n\times 2n}, \ \ \ a_{ij}=\omega_j(\beta_i), i, j=0, 1, \cdots, 2n, \\ \Phi=(\varphi_0, \varphi_1, \cdots, \varphi_{2n})^T, \ \ \ F=(f(\beta_0), f(\beta_1), \cdots, f(\beta_{2n}))^T,$

并且$\varphi_i(i=0, 1, \cdots, 2n)$表示函数$\varphi$在配置点$\beta_i$处的近似值.矩阵$\mathbb{A}$一般是奇异的, 为了求解线性方程组(2.13) 或(2.14), 通常需要选取某种正则化方法去求解.

对于弱奇异积分方程

$-\frac{1}{\pi}\int_0^{2\pi} \ln \left[2\sin \frac{\theta-\beta}{2}\right] \varphi(\theta) d\theta = f(\beta), \ \ \beta \in [0, 2\pi],$

(2.15)

利用三角插值积分公式(2.8) 逼近超奇异积分方程(2.15), 取一系列配置点$\beta_i (0\le i \le 2n)$并排列所得到的方程, 由此可以得到如下线性方程组

$\sum\limits_{j=0}^{2n} \omega_j '(\beta_i) \varphi_j = f(\beta_i), \ \ \ 0\le i \le 2n,$

(2.16)

或写成矩阵的形式

$\mathbb{B} \Phi = F,$

(2.17)

其中

$\mathbb{B}=(b_{ij})_{2n\times 2n}, \ \ \ b_{ij}=\omega_j '(\beta_i), i, j=0, 1, \cdots, 2n, \\ \Phi=(\varphi_0, \varphi_1, \cdots, \varphi_{2n})^T, \ \ \ F=(f(\beta_0), f(\beta_1), \cdots, f(\beta_{2n}))^T,$

并且$\varphi_i(i=0, 1, \cdots, 2n)$表示函数$\varphi$在配置点$\beta_i$处的近似值.矩阵$\mathbb{B}$一般也是奇异的, 为了求解它也需要正则化之后才能求解.

利用配置法求解(2.3) 或(2.4) 会得到两个线性方程组, 其矩阵形式为

$\mathbb{A} u = F,$

(2.18)

其中$u$是方程(2.3) 或(2.4) 中未知量近似值的列向量, $F$是由方程(2.3) 或(2.4) 右端表达式计算的边界$\Gamma_0$上的已知值.方程(2.18) 通常是一个严重的病态方程组, 假如边界$\Gamma_0$上没有足够多的信息或方程的右端项含有一定的噪声, 无法用方程(2.18) 求解另一部分边界$\Gamma_1$上所需要的信息.

为了求解方程(2.18), 采用Tiknonov正则化方法.令

$\hat{u}_{\alpha}=\arg \min \{\alpha ||\mathbb{A}u-F||^2+||Lu||^2\},$

(2.19)

其中$\alpha>0$是正则化参数, $L$是离散的一阶线性微分算子.上式的解可以写成矩阵的形式如下

$\hat{u}_{\alpha}=(\alpha L^T L+ \mathbb{A}^T \mathbb{A})^{-1}(\mathbb{A}^T F).$

(2.20)

求解这个方程即可得到边值反问题的近似解.这里关于最优正则化参数$\alpha$的选取可以参考文献[10].

下面给出两个数值算例来验证上节中提出的边值反问题的数值计算方法.

数值算例1 考虑边值反问题(1.1), 其中$u_0=3\cos 2x+4\sin 2x$, $u_n=6\cos 2x+8\sin 2x$, 分别计算了$\Gamma_0=[0, \pi]$, $\Gamma_0=[0, \pi/2]$和$\Gamma_0=[0, \pi/4]$三种情况, 其结果见图 1-3.

从图 1可以看出, 当$\Gamma_0=[0, \pi]$, 即已知边界占整个边界一半时, 由超奇异积分方程和弱奇异积分方程配置法复原得到的未知边界上的信息与其理论值高度一致, 这说明本文所提出的算法是有效的.从图 2和图 3的结果可以看出, 当已知边界小于一半时, 由配置法复原得到的未知边界上的值与理论值有一定的出入, 但近似计算结果仍然复原了原边界曲线上的波动特征.

在实际工程中, 通常测量数据都会有一定的误差, 即$\Gamma_0$上的已知值会含有一定的噪声.由于边值反问题已经被证明是一个严重的病态问题, 任何微小的扰动都会引起解的巨大改变.下面给出一个算例测试噪声对数值解得影响.

数值算例2 考虑例1中的边值反问题.分别计算了$\Gamma_0=[0, \pi/2]$和$\Gamma_0=[0, \pi/4]$两种情况, 同时对$\Gamma_0$上的已知信息$u_0$和$g_0$加入$2\%$的随机误差, 其结果见图 4-5.

从图 4和图 5的结果可以看出, 由积分方程配置法复原的数值解与精确解匹配还是非常好.同时把图 4和图 5的结果与没有添加噪声的复原结果相比, 其变化并不太大.这也从侧面印证了算法的有效性和可行性.