Helmholtz方程Cauchy问题的积分方程方法研究

Study on Integral Equation Method for the Cauchy Problem of the Helmholtz Equation

摘要: 本文研究了在环状区域上利用积分方程方法求解Helmholtz方程Cauchy问题。首先利用Green公式将方程解用积分方程表示，然后利用跳跃关系将区域内点趋于边界，得到两条边界上的积分方程组，采用单双层位势算子表示积分方程组，采用核裂解的方法离散积分方程组中的奇异积分，并使用Tikhonov正则化结合Morozov偏差原理求解病态方程。两种添加了噪声的数值算例验证了算法的有效性。

Abstract: In this paper, we study the problem of solving the Cauchy problem of Helmholtz’s equation on a toroidal region using the integral equation method. Firstly, the solution of the equation is expressed in terms of integral equations using Green’s formula, then the jump relation is used to converge the points in the region to the boundary to get the set of integral equations on the two boundaries, the single and double-layer potential operator is used to express the set of integral equations, the kernel cleavage method is used to discretize the singular integrals in the set of integral equations, and the pathological equation is solved by using the Tikhonov regularization combined with the Morozov’s deviation principle. Two numerical examples with added noise verify the effectiveness of the algorithm.

1. 引言

考虑二维环形区域Ω上Helmholtz方程的边值问题

$\{\begin{array}{l}\Delta u\left(x\right)+k^{2}u\left(x\right)=0,x\in \Omega \hfill \\ u\left(x\right)=f\left(x\right),x\in L_1\hfill \\ \frac{\partial u}{\partial \overrightarrow{n}}\left(x\right)=g\left(x\right),x\in L_1\hfill \end{array}$ (1.1)

其中，$L_0$ 和$L_1$ 是二维空间上的两条闭曲线，且满足$L_0\cap L_1=\varnothing$ ，$\partial \text{Ω }=L_0\cup L_1$ ，记$L_0$ 和$L_1$ 围成的区域为Ω，$\overrightarrow{n}$ 表示指向$\partial \text{Ω}$ 外的单位法向量。求解的问题为当已知(1.1)时，求${u|}_{L_0}$ 和${\frac{\partial u}{\partial \overrightarrow{n}}|}_{L_0}$ 。函数定义域Ω如图1所示。

Figure 1. Circular area Ω

图1. 环状区域Ω

对于问题(1.1)中，已知部分边界上的边界条件去求未知边界上的边界条件也称为反问题或Cauchy问题，Helmholtz方程Cauchy问题有几个重要的应用，例如，在检测材料内部缺陷、探测地下资源、非侵入式诊断(如肿瘤检测、组织特性分析)、桥梁、飞机、建筑等大型结构的损伤检测等。

1923年，著名数学家Hadamard提出了问题不适定性的概念，即方程的解需同时满足存在性、唯一性和稳定性。但由于现实技术、仪器噪声、数据精度等问题会导致初始条件存在误差，而问题(1.1)的不适定性主要体现在解的不稳定，即初始数据的任意小的波动都会导致最终解的巨大波动或解的震荡行为，针对于Helmholtz方程Cauchy问题，2003年Marin等将边界元与共轭梯度法相结合求解问题[1]，2008年Jin等利用平面波法求解该问题[2]，2011年Cheng等利用最优滤子法求解问题，并得到了Hölder型误差估计[3]，2011年Lin等利用基本解法和TR、TSVD和DSVD三种正则化方法求解方程[4]，2014年Berntsson等改进Kozlov的交替迭代法求解问题[5]，2018年孙瑶等利用间接积分方程方法求解Cauchy问题，并根据Morozov偏差原理选取正则化参数[6]，本文主要研究利用单双层位势函数得到方程解的表达式，并通过核裂解的方式避免核函数的奇异性，再利用正则化方法求解出方程的解，最后用数值实验验证方法的有效性。

2. 积分方程的导出

2.1. 理论准备

对于问题1.1，若Ω是$C^2$ 类有界域，$\overrightarrow{n}$ 表示指向$\partial \text{Ω}$ 外的单位法向量，$u\in C^2\left(\text{Ω}\right)\cap C^1\left(\overline{\text{Ω}}\right)$，可以利用Green公式去推导出解的积分方程形式：

$u\left(x\right)={\displaystyle {\int }_{\partial D}\left\{ \frac{\partial u}{\partial \overrightarrow{n}}\left(y\right)\Phi \left(x,y\right)-\frac{\partial \Phi \left(x,y\right)}{\partial \overrightarrow{n}\left(y\right)}u\left(y\right)\right\}\text{d}s\left(y\right)},x\in \text{Ω},$ (2.1)

这里$\frac{\partial u}{\partial \overrightarrow{n}}$ 定义为：

$\frac{\partial u}{\partial \overrightarrow{n}}\left(x\right)=\underset{h\to 0}{\lim }\overrightarrow{n}\left(x\right)\cdot \nabla u\left(x-h\overrightarrow{n}\left(x\right)\right),x\in \partial \text{Ω},$ (2.2)

其中，$\Phi \left(x,y\right)$ 为Helmholtz方程的基本解，表达式为：

$\Phi \left(x,y\right)=\frac{i}{4}{H}_0^{\left(1\right)}\left(k\left|x-y\right|\right),x,y\in R^2,x\ne y.$ (2.3)

定义1 设D是$C^2$ 类有界域，$\overrightarrow{n}$ 表示指向$\partial D$ 外的单位法向量，可积函数$\phi \in C\left(\partial D\right)$ ，我们称如下两个形式的积分

$\begin{array}{l}u\left(x\right):={\displaystyle {\int }_{\partial D}\phi \left(y\right)\Phi \left(x,y\right)\text{d}s\left(y\right)},x\in R^2\backslash \partial D,\\ v\left(x\right):={\displaystyle {\int }_{\partial D}\phi \left(y\right)\frac{\partial \Phi \left(x,y\right)}{\partial \overrightarrow{n}\left(y\right)}\text{d}s\left(y\right)},x\in R^2\backslash \partial D,\end{array}$ (2.4)

分别为密度函数$\phi$ 的单层位势和双层位势。

定理1 [7] (跳跃关系)设$\partial D$ 是$C^2$ 类有界域，$\overrightarrow{n}$ 表示指向$\partial D$ 外的单位法向量，可积函数$\phi \in C\left(\partial D\right)$ ，以$\phi$ 为密度函数的单层位势u在全空间$R^2$ 上连续，存在一个连续依赖于$\partial D$ 的常数C满足

${\Vert u\Vert }_{\infty ,R^2}\le C{\Vert \phi \Vert }_{\infty ,R^2}$ (2.5)

那么在边界上就有

$u\left(x\right)=\underset{a\to x}{\lim }u\left(a\right)={\displaystyle {\int }_{\partial D}\phi \left(y\right)\Phi \left(x,y\right)\text{d}s\left(y\right)},x\in \partial D,$ (2.6)

以$\phi$ 为密度函数的双层位势v在极限意义下，可以从D连续延拓到$\overline{D}$ 和从$R^2\backslash D$ 连续延拓到$R^2\backslash \overline{D}$ ，有

$v_{\pm }\left(x\right)={\displaystyle {\int }_{\partial D}\phi \left(y\right)\frac{\partial \Phi \left(x,y\right)}{\partial \overrightarrow{n}\left(y\right)}\text{d}s\left(y\right)}\pm \frac{1}{2}\phi \left(x\right),x\in \partial D,$ (2.7)

其中

$v_{\pm }\left(x\right):=\underset{h\to 0^{+}}{\lim }v\left(x\pm h\overrightarrow{n}\left(x\right)\right),x\in \partial D,$ (2.8)

以反常积分的形式存在。

2.2. 边界方程组的导出

回到问题1.1，Helmholtz方程解可表示为：

$\begin{array}{l}u\left(x\right)={\displaystyle {\int }_{L_1}\left\{ \frac{\partial u}{\partial \overrightarrow{n}}\left(y\right)\Phi \left(x,y\right)-\frac{\partial \Phi \left(x,y\right)}{\partial \overrightarrow{n}\left(y\right)}u\left(y\right)\right\}\text{d}s\left(y\right)}\\ -{\displaystyle {\int }_{L_0}\left\{\frac{\partial u}{\partial \overrightarrow{n}}\left(y\right)\Phi \left(x,y\right)-\frac{\partial \Phi \left(x,y\right)}{\partial \overrightarrow{n}\left(y\right)}u\left(y\right)\right\}\text{d}s\left(y\right)},x\in \Omega .\end{array}$ (2.9)

这里，$\overrightarrow{n}$ 表示指向$\partial D$ 外的单位量。

在上式中，任取$z\in \Omega$ ，令$x\to z\in L_0$ 时，右端在$L_1$ 上的积分式没有奇异性，但在$L_0$ 上的积分存在奇异性，所以利用跳跃关系[7]和边界条件有：

$\begin{array}{l}u\left(x\right)+2{\displaystyle {\int }_{L_0}\left\{\frac{\partial u}{\partial \overrightarrow{n}}\left(y\right)\Phi \left(x,y\right)-\frac{\partial \Phi \left(x,y\right)}{\partial \overrightarrow{n}\left(y\right)}u\left(y\right)\right\}\text{d}s\left(y\right)}\\ =2{\displaystyle {\int }_{L_1}\left\{g\left(y\right)\Phi \left(x,y\right)-\frac{\partial \Phi \left(x,y\right)}{\partial \overrightarrow{n}\left(y\right)}f\left(y\right)\right\}\text{d}s\left(y\right)},x\in L_0\end{array}$ (2.10)

同理利用跳跃关系[7]和边界条件有：

$\begin{array}{l}2{\displaystyle {\int }_{L_1}\left\{\frac{\partial u}{\partial \overrightarrow{n}}\left(y\right)\Phi \left(x,y\right)-\frac{\partial \Phi \left(x,y\right)}{\partial \overrightarrow{n}\left(y\right)}u\left(y\right)\right\}\text{d}s\left(y\right)}\\ = 2{\displaystyle {\int }_{L_0}\left\{g\left(y\right)\Phi \left(x,y\right)-\frac{\partial \Phi \left(x,y\right)}{\partial \overrightarrow{n}\left(y\right)}f\left(y\right)\right\}\text{d}s\left(y\right)}-f\left(x\right),x\in L_1.\end{array}$ (2.11)

设可积函数$\phi \in C\left(\partial \text{Ω}\right)$ ，我们定义下列积分算子：

$\begin{array}{l}{S}^{i,j}\left[\phi \right]\left(x\right):=2{\displaystyle {\int }_{L_i}\Phi \left(x,y\right)\phi \left(y\right)\text{d}s\left(y\right)},x\in L_j,\\ K^{i,j}\left[\phi \right]\left(x\right):=2{\displaystyle {\int }_{L_i}\frac{\partial \Phi \left(x,y\right)}{\partial \overrightarrow{n}\left(y\right)}\phi \left(y\right)\text{d}s\left(y\right)},x\in L_j.\end{array}$ (2.12)

其中，$i,j=0,1$ 。设I是单位算子，我们可以将积分方程组用算子表示为：

$\{\begin{array}{c}\left(I-K^{0,0}\right)u+S^{0,0}\frac{\partial u}{\partial \overrightarrow{n}}=S^{1,0}g-K^{1,0}f,\\ S^{0,1}\frac{\partial u}{\partial \overrightarrow{n}}-K^{0,1}u=S^{1,1}g-\left(K^{1,1}+I\right)f.\end{array}$ (2.13)

注意到$K^{0,0}$ 的核函数具有弱奇异性，根据文献[8]，可知算子$I-K^{0,0}$ 是可逆算子，由表达式(2.10)可以得到：

$u={\left(I-K^{0,0}\right)}^{-1}\left(S^{1,0}g-K^{1,0}f-S^{0,0}\frac{\partial u}{\partial \overrightarrow{n}}\right).$ (2.14)

再将表达式(2.14)代入表达式(2.11)，可以得到：

$\begin{array}{l}{S}^{0,1}\frac{\partial u}{\partial \overrightarrow{n}}-K^{0,1}{\left(I-K^{0,0}\right)}^{-1}\left(S^{1,0}g-K^{1,0}f-S^{0,0}\frac{\partial u}{\partial \overrightarrow{n}}\right)\\ =S^{1,1}g-\left(K^{1,1}+I\right)f,\end{array}$ (2.15)

整理表达式(2.15)可以得到：

$\begin{array}{l}\left[S^{0,1}+K^{0,1}{\left(I-K^{0,0}\right)}^{-1}{S}^{0,0}\right]\frac{\partial u}{\partial \overrightarrow{n}}\\ =S^{1,1}g-\left(K^{1,1}+I\right)f+K^{0,1}{\left(I-K^{0,0}\right)}^{-1}\left(S^{1,0}g-K^{1,0}f\right).\end{array}$ (2.16)

从积分方程组可以看出，方程(2.11)是关于$u$ 的第二类Fredholm积分方程，而方程(2.10)是关于$\frac{\partial u}{\partial \overrightarrow{n}}$ 的第一类积分方程，我们可以利用正则化理论求解方程(2.16)。再利用求解出$\frac{\partial u}{\partial \overrightarrow{n}}$ 代入式(2.16)可以计算出$u$ 。

3. 数值方法求解

3.1. 积分算子离散化

首先，我们假定数值算例的曲线充分光滑，定义二维空间上的曲线方程

$\begin{array}{l}\Gamma :=\left\{x:x\left(t\right)=\left(x_1\left(t\right),x_2\left(t\right)\right),t\in \left[0,2\pi \right]\right\},\\ \Gamma :=\left\{y:y\left(\tau \right)=\left(y_1\left(\tau \right),y_2\left(\tau \right)\right),\tau \in \left[0,2\pi \right]\right\}.\end{array}$ (3.1)

将积分算子参数化得到：

$\begin{array}{l}{S}^{i,j}\left[\phi \right]\left(x\right)=2{\displaystyle {\int }_0^{2\pi }\Phi \left(x\left(t\right),y\left(\tau \right)\right)\left|y^{\prime }\left(\tau \right)\right|\phi \left(y\left(\tau \right)\right)\text{d}\tau },\\ K^{i,j}\left[\phi \right]\left(x\right)=2{\displaystyle {\int }_0^{2\pi }\frac{\partial \Phi \left(x\left(t\right),y\left(\tau \right)\right)}{\partial \overrightarrow{n}\left(y\left(\tau \right)\right)}\left|y^{\prime }\left(\tau \right)\right|\phi \left(y\left(\tau \right)\right)\text{d}\tau }.\end{array}$ (3.2)

记

$\begin{array}{c}M\left(t,\tau \right)=2\cdot \Phi \left(x\left(t\right),y\left(\tau \right)\right)\left|y^{\prime }\left(\tau \right)\right|=\frac{i}{2}{H}_0^{\left(1\right)}\left(k\left|x\left(t\right)-y\left(\tau \right)\right|\right)\left|y^{\prime }\left(\tau \right)\right|,\\ K\left(t,\tau \right)=2\cdot \frac{\partial \Phi \left(x\left(t\right),y\left(\tau \right)\right)}{\partial \overrightarrow{n}\left(y\left(\tau \right)\right)}\left|y\text{'}\left(\tau \right)\right|\\ =\frac{ki}{2}{H}_1^{\left(1\right)}\left(k\left|x\left(t\right)-y\left(\tau \right)\right|\right) \left|y\text{'}\left(\tau \right)\right|\overrightarrow{n}\left(y\left(\tau \right)\right)\cdot \frac{\left(x\left(t\right)-y\left(\tau \right)\right)}{\left|x\left(t\right)-y\left(\tau \right)\right|}.\end{array}$ (3.3)

其中，$\overrightarrow{n}$ 表示单位外法向量，$H_0^{\left(1\right)}$ 和$H_1^{\left(1\right)}$ 分别为第一类零阶和一阶Hankel函数，且有${\left(H_0^{\left(1\right)}\right)}^{\prime }=-H_1^{\left(1\right)}$ 。

当$x,y\in L_j$ 时，$M\left(t,\tau \right)$ 和$K\left(t,\tau \right)$ 在$t=\tau$ 处有奇异性，针对核函数$M\left(t,\tau \right)$ 和$K\left(t,\tau \right)$ 的奇异性，根据文献[8]，我们对核函数进行裂解：

$\begin{array}{l}M\left(t,\tau \right)=M_1\left(t,\tau \right)\ln \left(4{\sin }^2\frac{t-\tau }{2}\right)+M_2\left(t,\tau \right),\\ K\left(t,\tau \right)=K_1\left(t,\tau \right)\ln \left(4{\sin }^2\frac{t-\tau }{2}\right)+K_2\left(t,\tau \right),\end{array}$ (3.4)

其中

$\begin{array}{l}{M}_1\left(t,\tau \right)=\{\begin{array}{l}-\frac{1}{2\pi }{J}_0\left(k\left|x\left(t\right)-y\left(\tau \right)\right|\right)\left|y^{\prime }\left(\tau \right)\right|,t\ne \tau \hfill \\ 0,t=\tau \hfill \end{array},\\ K_1\left(t,\tau \right)=\{\begin{array}{l}-\frac{k}{2\pi }{J}_1\left(k\left|x\left(t\right)-y\left(\tau \right)\right|\right) \left|y^{\prime }\left(\tau \right)\right|\overrightarrow{n}\left(y\left(\tau \right)\right)\cdot \frac{\left(x\left(t\right)-y\left(\tau \right)\right)}{\left|x\left(t\right)-y\left(\tau \right)\right|} ,t\ne \tau \hfill \\ 0,t=\tau \hfill \end{array}\end{array}$ (3.5)

其中，$J_0$ 和$J_1$ 分别为零阶和一阶Bessel函数。

$M_2\left(t,\tau \right)$ 和$K_2\left(t,\tau \right)$ 表达式为：

$M_2\left(t,\tau \right)=M\left(t,\tau \right)-M_1\left(t,\tau \right)\ln \left(4{\sin }^2\frac{t-\tau }{2}\right),$

$K_2\left(t,\tau \right)=K\left(t,\tau \right)-K_1\left(t,\tau \right)\ln \left(4{\sin }^2\frac{t-\tau }{2}\right).$

特别地，当$t=\tau$ 时，有

$M_2\left(t,t\right)=\left\{\frac{i}{2}-\frac{\gamma }{\pi }-\frac{1}{\pi }\ln \frac{k\left|x^{\prime }\left(t\right)\right|}{2}\right\}\left|x^{\prime }\left(t\right)\right|,$

$K_2\left(t,t\right)=\frac{1}{2\pi }\frac{\left( {x^{″}}_1\left(t\right){x^{\prime }}_2\left(t\right)-{x^{″}}_2\left(t\right){x^{\prime }}_1\left(t\right)\right)}{{\left|x^{\prime }\left(t\right)\right|}^2}.$

其中，$\gamma$ 为Euler常数。

对连续积分利用复合矩形公式离散，奇异积分利用三角插值多项式[9]，

$\begin{array}{l}{\displaystyle {\int }_0^{2\pi }\phi \left(\tau \right)\text{d}\tau }\approx \frac{\pi }{N}{\displaystyle \sum }_{j=0}^{2N-1}\phi \left({\tau }_j\right),\\ {\displaystyle {\int }_0^{2\pi }\phi \left(\tau \right)\ln \left(4{\sin }^2\frac{t-\tau }{2}\right)\text{d}\tau }\approx {\displaystyle \sum }_{j=0}^{2N-1}{R}_j^N\left(t\right)\phi \left({\tau }_j\right).\end{array}$ (3.6)

其中

$R_j^N\left(t\right)=-\frac{\pi }{N}\left[{\displaystyle \sum }_{m=1}^{N-1}\frac{1}{m}\cos \left(m\left(t-{\tau }_j\right)\right)+\frac{1}{2N}\cos N\left(t-{\tau }_j\right)\right].$ (3.7)

记$u_i=u\left(x\left(t_i\right)\right)$ ，由于$x,y\in L_m$ ，$M\left(t,\tau \right)$ 存在奇异性，所以将$M\left(t,\tau \right)$ 核裂解可以得到：

$\begin{array}{c}{S}^{m,m}{u}_i={\displaystyle {\int }_0^{2\pi }\left\{M_1\left(t_i,\tau \right)\ln \left(4{\sin }^2\frac{t_i-\tau }{2}\right)u\left(y\left(\tau \right)\right)+M_2\left(t_i,\tau \right)u\left(y\left(\tau \right)\right)\right\}\text{d}\tau }\\ ={\displaystyle \sum }_{j=0}^{2N-1}\left\{R_j^N\left(t_i\right)M_1\left(t_i,{\tau }_j\right)+\frac{\pi }{N}{M}_2\left(t_i,{\tau }_j\right)\right\}{u}_j.\end{array}$ (3.8)

记$U={\left(u_0,u_1,\cdots ,u_{2N-1}\right)}^{\text{T}},P^{m,m}={\left(p_{ij}^{m,m}\right)}_{2N\times 2N}$ 。遍历$i=0,1,\cdots ,2N-1$ ，可以得到算子参数化的矩阵形式：

$S^{m,m}u\approx P^{m,m}U,x,y\in L_m$ (3.9)

其中，$p_{ij}^{m,m}=R_j^N\left(t_i\right)M_1\left(t_i,{\tau }_j\right)+\frac{\pi }{N}{M}_2\left(t_i,{\tau }_j\right)$ 。

同理可以得到$K^{m,m}\left(m=0,1\right)$ 的参数化的矩阵形式为：

$K^{m,m}u\approx Q^{m,m}U ,x,y\in L_m$ (3.10)

其中，$Q^{m,m}={\left(q_{ij}^{m,m}\right)}_{2N\times 2N}$ ，$q_{ij}^{m,m}=R_j^N\left(t_i\right)K_1\left(t_i,{\tau }_j\right)+\frac{\pi }{N}{K}_2\left(t_i,{\tau }_j\right)$ 。

由于$x\in L_m,y\in L_n,M\left(t,\tau \right)$ 不存在奇异性，对$S^{m,n}{u}_i$ 采用复合矩形离散方法：

$S^{m,n}{u}_i={\displaystyle {\int }_0^{2\pi }M\left(t_i,\tau \right)u\left(y\left(\tau \right)\right)\text{d}\tau }={\displaystyle \sum }_{j=0}^{2N-1}\frac{\pi }{N}M\left(t_i,{\tau }_j\right)u_j,j=0,1,\cdots ,2N-1.$ (3.11)

算子参数化的矩阵形式：

$S^{m,n}u\approx P^{m,n}U,x\in L_m,y\in L_n.$ (3.12)

此时$p_{ij}^{m,n}=\frac{\pi }{N}M\left(t_i,{\tau }_j\right)$ 。同理$q_{ij}^{m,n}=\frac{\pi }{N}K\left(t_i,{\tau }_j\right)$ 。

3.2. 积分方程离散化

对于积分方程组

$\{\begin{array}{c}\left(I-K^{0,0}\right)u+S^{0,0}\frac{\partial u}{\partial \overrightarrow{n}}=S^{1,0}g-K^{1,0}f,\\ S^{0,1}\frac{\partial u}{\partial \overrightarrow{n}}-K^{0,1}u=S^{1,1}g-\left(K^{1,1}+I\right)f.\end{array}$ (3.13)

将其离散化可以得到：

$\{\begin{array}{c}{A}_{1}U+B_{1}V=F_1\\ A_{2}U+B_{2}V=F_2\end{array},$ (3.14)

这里，$U,V$ 是待求结果，$F_1,F_2$ 为常数向量。记$f_i=f\left(y\left({\tau }_i\right)\right),g_i=g\left(y\left({\tau }_i\right)\right),F={\left\{f_i\right\}}_{2N\times 1},G={\left\{g_i\right\}}_{2N\times 1}$ ，对照离散化矩阵方程组和积分方程组组系数可得：

$\{\begin{array}{c}{A}_1=\left(I-Q^{0,0}\right) ,B_1=P^{0,0},F_1=P^{1,0}G-Q^{1,0}F,\\ A_2=-Q^{0,1} ,B_2=P^{0,1},F_2=P^{1,1}g-\left(Q^{1,1}+I\right)F.\end{array}$ (3.15)

由$\left(I-K^{0,0}\right)$ 是可逆算子，推广到矩阵可知$A_1$ 是可逆矩阵，由此可以得到：

$U=A_1^{-1}\left(F_1-B_{1}V\right),$ (3.16)

代入$A_{2}U+B_{2}V=F_2$ ，可以得到：

$\left(B_2-A_2{A}_1^{-1}{B}_1\right)V=F_2-A_2{A}_1^{-1}{F}_1.$ (3.17)

记$C=B_2-A_2{A}_1^{-1}{B}_1,d=F_2-A_2{A}_1^{-1}{F}_1$ ，方程组(3.16)、(3.17)可写作：

$\{\begin{array}{l}CV=d,\hfill \\ U=A_1^{-1}\left(F_1-B_{1}V\right).\hfill \end{array}$ (3.18)

待求解的方程组(3.14)等价于求解方程组(3.18)。

3.3. 正则化

将积分方程组离散得到的方程组(3.18)中$CV=d$ 是严重病态方程，当数据包含噪声时，右端项的波动会造成解的剧烈震荡，为了求解方程(3.18)，本文采用Tikhonov正则化，令

$V_{\alpha }=\arg \min \left\{{\Vert CV-d\Vert }^2+\alpha {\Vert Iu\Vert }^2\right\}.$ (3.19)

上式的解$V_{\alpha }$ 为：

$V_{\alpha }={\left(C^{T}C+\alpha I\right)}^{-1}\cdot \left(C^{T}d\right).$ (3.20)

求解出的$V_{\alpha }$ 就是原边值问题的近似解，这里最优正则化参数的选取参照Morozov偏差原理[10]，由此，我们可以得到U的表达式：

$U=A_1^{-1}\left(F_1-B_1{V}_{\alpha }\right).$ (3.21)

3.4. 误差分析

对于本文数值求解的问题，误差来源于离散积分时导致的结果误差：(3.6)中的连续型积分离散和奇异积分对核函数三角插值引入的误差，对于这个误差，文献[8]中指出了解的收敛性(即数值解收敛于解析解)，并给出了误差估计：

假设积分方程本身有唯一解，且核函数以及边界条件连续，那么有

1) 对于所有足够大的n，近似线性系统(3.13)，即近似方程(3.14)有唯一解；

2) 当$n\to \infty$ 时，近似解$u^n\left(t\right)$ 一致收敛于积分方程的解$u\left(t\right)$ ；

3) (3.6)的求积分误差的收敛阶传递到$u^n\left(t\right)$ ~$u\left(t\right)$ 上；

特别地，第三条意味着解析核函数以及解析边界条件连续的情况下，近似误差与三角插值误差一致为$O\left(e^{-n\sigma }\right)$ ，这里$\sigma$ 表示表示实解析函数可以全纯延拓到复平面上的平行带宽度的一半。即存在正常数C和$\sigma$ 使得

$\left|u^n\left(t\right)-u\left(t\right)\right|\le C{e}^{-n\sigma },0\le t\le 2\pi ,$

对于所有n成立。原则上，常数C是可计算的，但通常很难求值。在大多数实际情况中，通过将节点数2n加倍，然后借助指数收敛阶比较粗网格和细网格的结果，即根据节点数2n加倍会使近似解中正确数字的位数加倍这一事实，来判断计算解的精度就足够了。

而算法的稳定性依赖于Tikhonov正则化，也就是说数值方法的稳定性只需分析正则化算法即可，具体分析可见文献[11]。

4. 数值算例

4.1. 数值样例

下面给出两个数值算例验证Helmholtz方程Cauchy问题(1.1)，由于现实问题中的数据往往是含噪声的，我们对第一、第二类边界条件添加噪声。

例1 考虑求解边界$L_1:=\left\{\left(x_1\left(t\right),x_2\left(t\right)\right)=\left(5\cos \left(t\right),5\sin \left(t\right)\right),t\in \left[0,2\pi \right]\right\},$ $L_0:=\left\{\left(x_1\left(t\right),x_2\left(t\right)\right)=\left(3\cos \left(t\right),2\sin \left(t\right)\right),t\in \left[0,2\pi \right]\right\}$ 。

例2 考虑求解边界$L_1:=\left\{\left(x_1\left(t\right),x_2\left(t\right)\right)=(5\cos \left(t\right)-\sin \left(t\right),5\sin \left(t\right)+\cos \left(t\right),t\in \left[0,2\pi \right]\right\},$ $L_0:=\left\{\left(x_1\left(t\right),x_2\left(t\right)\right)=\left(3\cos \left(t\right)+1.5\sin \left(t\right),3\sin \left(t\right)\right),t\in \left[0,2\pi \right]\right\}$ 。

本文分别计算了波数k = 3和k = 13时的边界复原结果，并计算了添加噪声$\delta =1\%,3\%,5\%$ 时的结果，由于不添加噪声时数值解和解析解的相对误差只有0.01%，这里不再给出。

从图2~5可以看出，在噪声水平比较低的时候数值解与理论值一致，噪声水平增大，数值解在尖点处数据出现微小波动。从数值结果上可以印证算法的有效性、可行性和鲁棒性。

(a) $u$ (b) $\frac{\partial u}{\partial \overrightarrow{n}}$

Figure 2. Example 1: Numerical and analytical solution images for the upper added noise $\delta =1\%,3\%,5\%$ on L₀ at k = 3

图2. 例1：k = 3时L₀上添加噪声$\delta =1\%,3\%,5\%$ 的数值解与解析解图像

(a) $u$ (b) $\frac{\partial u}{\partial \overrightarrow{n}}$

Figure 3. Example 1: Numerical and analytical solution images for the upper added noise $\delta =1\%,3\%,5\%$ on L₀ at k = 10

图3. 例1：k = 10时L₀上添加噪声$\delta =1\%,3\%,5\%$ 的数值解与解析解图像

(a) $u$ (b) $\frac{\partial u}{\partial \overrightarrow{n}}$

Figure 4. Example 2: Numerical and analytical solution images for the upper added noise $\delta =1\%,3\%,5\%$ on L₀ at k = 3

图4. 例2：k = 3时L₀上添加噪声$\delta =1\%,3\%,5\%$ 的数值解与解析解图像

(a) $u$ (b) $\frac{\partial u}{\partial \overrightarrow{n}}$

Figure 5. Example 2: Numerical and analytical solution images for the upper added noise $\delta =1\%,3\%,5\%$ on L₀ at k = 10

图5. 例2：k = 10时L₀上添加噪声$\delta =1\%,3\%,5\%$ 的数值解与解析解图像

(a) $u$ (b) $\frac{\partial u}{\partial \overrightarrow{n}}$

Figure 6. Example 3: Numerical and analytical solution images for the upper added noise $\delta =1\%,3\%,5\%$ on L₀ at k = 5

图6. 例3：k = 5时L₀上添加噪声$\delta =1\%,3\%,5\%$ 的数值解与解析解图像

(a) $u$ (b) $\frac{\partial u}{\partial \overrightarrow{n}}$

Figure 7. Example 3: Numerical and analytical solution images for the upper added noise $\delta =1\%,3\%,5\%$ on L₀ at k = 6

图7. 例3：k = 6时L₀上添加噪声$\delta =1\%,3\%,5\%$ 的数值解与解析解图像

但算法在波数较大或复杂区域的拟合程度会有一定的下降，如例3考虑求解边界$L_1:=\left\{\left(x_1\left(t\right),x_2\left(t\right)\right)=(6\cos \left(t\right)+3.9\cos \left(2t\right)-3.9,9\sin \left(t\right),t\in \left[0,2\pi \right]\right\},$ $L_0:=\left\{\left(x_1\left(t\right),x_2\left(t\right)\right)=\left(3\cos \left(t\right)+1.5\sin \left(t\right),3\sin \left(t\right)\right),t\in \left[0,2\pi \right]\right\}$ 。当k ≤ 5.5时，数值解与解析解误差较小。但是当k ≥ 6时，数值解和解析解误差增大。

从图6、图7可以看出，此数值解法对于不同边界上的同波数误差情况也不相同，在波数增大后数据误差变大。

4.2. 结果分析

设$u^{an}$ 和$u^{num}$ 分别表示解析解和数值解，那么我们可以定义第一、二类边界条件上的标准化相对RMS误差为：

$E=\frac{\sqrt{{{\displaystyle \sum }}_{i=0}^{2N-1}{\left|u^{an}\left(t_i\right)-u^{num}\left(t_i\right)\right|}^2}}{\max \left|u^{an}\left(t_i\right)\right|}.$

对于正则化参数选取，本文采用Morozov偏差原理给出，为了验证Morozov偏差原理的有效性，计算区域选用例1，本文给出在初始数据噪声$\delta =1\%,3\%,5\%$ 时，正则化参数$\alpha ={10}^{-2},{10}^{-3},{10}^{-4},{10}^{-5}$ 和Morozov偏差原理计算出的正则化参数下的相对RMS误差。表中最后一列表示Morozov偏差原理计算出的正则化参数。

Table 1. Relative error for noise δ = 1%

表1. 噪声δ = 1%时相对误差

正则化参数	10^−2	10^−3	10^−4	10^−5	1.81 × 10^−3
k = 3	0.0339	0.0054	0.0101	0.0185	0.0049
k = 10	0.2049	0.0295	0.0080	0.0208	0.0102

Table 2. 噪声δ = 3%时相对误差

表2. Relative error for noise δ = 3%

正则化参数	10^−2	10^−3	10^−4	10^−5	5.5×10^−3
k = 3	0.0336	0.0095	0.0253	0.0449	0.0213
k = 10	0.2057	0.0341	0.0480	0.0703	0.0296

Table 3. Relative error for noise δ = 5%

表3. 噪声δ = 5%时相对误差

正则化参数	10^−2	10^−3	10^−4	10^−5	9.0×10^−3
k = 3	0.0370	0.0203	0.0391	0.1349	0.0166
k = 10	0.2059	0.0374	0.0544	0.0891	0.0289

从表1~3中可以看出，不同噪声的Morozov偏差原理计算出的正则化参数不相同，但其相对误差均最小，同时正则化参数呈现先下降后上升的趋势，在最优参数处相对误差最小。