基于牛顿型算法的椭圆方程多参数识别

Multi-Parameter Identification of Elliptic Equations Based on Newton-Type Algorithms

摘要: 针对椭圆型偏微分方程(PDE) Neumann边值问题中的多参数联合识别问题展开研究，重点解决在状态测量数据精度不足的条件下，同时辨识扩散矩阵、源项、边界条件及物理状态的挑战。建立了一个结合能量泛函与Tikhonov正则化的变分框架，推导出Karush-Kuhn-Tucker (KKT)系统，并采用有限元方法对PDE进行空间离散，发展牛顿型算法以实现数值求解。我们的算法在不同测量数据条件下能够有效地辨识多个参数，最终的数值算例验证了该方法的有效性。

Abstract: The study focuses on the multi-parameter joint identification problem in the Neumann boundary value problem for elliptic partial differential equations (PDEs). It specifically addresses the challenge of simultaneously identifying the diffusion matrix, source term, boundary conditions, and physical states under conditions of insufficient accuracy in state measurement data. We establish a variational framework that combines an energy functional with Tikhonov regularization and derive the Karush-Kuhn-Tucker (KKT) system. We also employ the finite element method for spatial discretization of the PDE and develop a Newton-type algorithm for numerical solution. Our algorithm effectively identifies multiple parameters under different measurement data conditions, and the final numerical examples validate the effectiveness of the method.

1. 引言

偏微分方程(PDE)参数识别问题是一类核心的数学物理反问题，其目标是通过系统输出的观测数据来反演模型中的未知参数、源项或初始边界条件[1]。这类问题广泛存在于含水层分析、地球物理勘探、医学成像及环境污染物传输等关键科学领域，其研究在过去三十余年里持续受到广泛关注[2]-[14]。传统的参数识别研究多集中于标量扩散系数，其理论与数值方法已相对成熟。然而，实际系统往往更为复杂，涉及多参数(如对流速度与源强度)或各向异性扩散矩阵的同时识别，这为反问题研究带来了巨大的理论与计算挑战。系统的源项(source term)难以直接测量或已知甚少，与大量假设源项已知的研究不同，源项本身作为未知量与扩散系数、边界条件一同被识别，构成了一个更具普适性但也不适定性更强的反问题。这类问题解的存在性、唯一性及稳定性分析都更为困难，需要发展专门的正则化理论与稳定的数值算法[15]。在数值求解这类优化问题时，基于一阶最优性条件(Karush-Kuhn-Tucker条件)构建的算法是主流途径。其中，牛顿型算法因其在最优解附近具有二阶收敛速度，相比一阶梯度法能显著提升计算效率，尤其适用于大规模或非线性程度高的参数识别问题[16]。

本文借鉴文献[17]在多参数识别研究上的延续与发展，针对扩散矩阵、未知源项及边界条件需同时识别的这一更具一般性的问题，采用一种基于代价泛函凸性分析的新技术路径。该分析为构建稳定的迭代格式提供了理论基础。构建一种高效的牛顿型算法来数值求解由此导出的KKT系统，旨在实现比传统一阶方法更快的收敛速度。后续内容将围绕问题的数学表述、理论分析、算法设计及数值实验展开。

2. 问题介绍

设区域$\Omega$ 为$R^d$ 的开且有界的连通域，其中$1\le d\le 3$ ，并且其边界$\partial \Omega$ 为多边形边界。通过解$u$ 的测量数据$z_{\delta }\in L^2\left(\Omega \right)$ ，本研究关注在椭圆型偏微分方程的Neumann边值问题(1.1)中，扩散矩阵$A$ 、源项$f$ 、边界条件$g$ 以及状态$u$ 的同时识别问题：

$\{\begin{array}{l}-\nabla \cdot \left(A\nabla u\right)=f\text{in}\Omega ,\\ A\nabla u\cdot \overrightarrow{n}=g\text{on}\partial \Omega ,\end{array}$ (1.1)

其中$\overrightarrow{n}$ 是在$\partial \Omega$ 上的单位外法向量，$\left(A,f,g\right)\in {\mathscr{H}}_{ad}:={\mathcal{A}}_{ad}\times {ℱ}_{ad}\times {\mathcal{G}}_{ad}$ 。这里${\mathcal{A}}_{ad},{ℱ}_{ad}$ 和${\mathcal{G}}_{ad}$ 分别定义如下：

${\mathcal{A}}_{ad}:=\left\{A\in L_{sym}^{\infty }\left(\Omega \right)|\alpha {\left|\xi \right|}^2\le A\left(x\right)\xi \cdot \xi \le \beta {\left|\xi \right|}^2,\forall \xi \in R^d\right\}$ ,(1.2)

${ℱ}_{ad}:=L^2\left(\Omega \right)$ ，${\mathcal{G}}_{ad}:=L^2\left(\partial \Omega \right)$ ，

其中$\alpha ,\beta$ 为满足$\beta \ge \alpha >0$ 的常数。

对于$1\le p\le \infty$ ，设

$L_{sym}^p\left(\Omega \right):=\left\{H:={\left(h_{ij}\right)}_{i,j=\overline{1,\cdots ,d}}\in L^p{\left(\Omega \right)}^{d\times d}|H\left(x\right):={\left(h_{ij}\left(x\right)\right)}_{i,j=\overline{1,\cdots ,d}}\in S_d在\Omega 中几乎��成立\right\}$

式中，$S_d$ 定义为表示具有内积$M\cdot N:=\text{trace}\left(MN\right)$ 运算和相应范数${\Vert M\Vert }_{S_d}={\left(M\cdot M\right)}^{1/2}={\left({\displaystyle {\sum }_{i,j=1}^d{m}_{ij}^2}\right)}^{1/2}$ 的所有$d\times d$ 实对称矩阵的集合，其中$M:={\left(m_{ij}\right)}_{i,j=\overline{1,\cdots ,d}}$ 。

在$L_{sym}^p\left(\Omega \right)$ 中，我们使用标量积和范数

${\left(H^1,H^2\right)}_{L_{sym}^p\left(\Omega \right)}={\displaystyle {\sum }_{i,j=1}^d{\left(h_{ij}^1,h_{ij}^2\right)}_{L^2\left(\Omega \right)}}$

${\Vert H\Vert }_{L_{sym}^2\left(\Omega \right)}:={\left({\displaystyle {\sum }_{i,j=1}^d{\Vert h_{ij}\Vert }_{L^2\left(\Omega \right)}^2}\right)}^{\frac{1}{2}}={\left({\displaystyle {\int }_{\Omega }{\Vert H\left(x\right)\Vert }_{S_d}^2\text{d}x}\right)}^{\frac{1}{2}}$ ，

而空间$L_{sym}^{\infty }\left(\Omega \right)$ 配备的范数为${\Vert H\Vert }_{L_{sym}^{\infty }\left(\Omega \right)}:={\max }_{i,j=\overline{1,\cdots ,d}}{\Vert h_{ij}\Vert }_{L^{\infty }\left(\Omega \right)}$ 。

设$\gamma :H^1\left(\Omega \right)\to H^{1/2}\left(\partial \Omega \right)$ 为连续Dirichlet迹算子，$H_{\diamond }^1\left(\Omega \right)$ 为$H^1\left(\Omega \right)$ 的一个闭子空间，由边界上所有均值为零的函数组成，即：

$H_{\diamond }^1\left(\Omega \right):=\left\{u\in H^1\left(\Omega \right){\displaystyle {\int }_{\partial \Omega }\gamma u\text{d}x}=0\right\}$

以此构建完备的Hilbert空间。

引理1.1 对于每个$\left(A,f,g\right)\in {\mathscr{H}}_{ad}$ ，存在唯一的弱解$u$ 满足(1.1)，且存在以下先验估计：

${\Vert u\Vert }_{H^1\left(\Omega \right)}\le C\max \left({\Vert \gamma \Vert }_{ℒ\left(H^1\left(\Omega \right),H^{1/2}\left(\partial \Omega \right)\right)}{\Vert g\Vert }_{L^2\left(\partial \Omega \right)}+{\Vert f\Vert }_{L^2\left(\Omega \right)}\right)\le C_{\begin{array}{c}N\end{array}}\left({\Vert g\Vert }_{L^2\left(\partial \Omega \right)}+{\Vert f\Vert }_{L^2\left(\Omega \right)}\right)$ (1.3)

其中，$C_{\mathcal{N}}:=C\max \left(1,{\Vert \gamma \Vert }_{ℒ\left(H^1\left(\Omega \right),H^{1/2}\left(\partial \Omega \right)\right)}\right)$ 。

证明：由Lax-Milgram引理可得，对于每个$\left(A,f,g\right)\in {\mathscr{H}}_{ad}$ ，$v\in H_{\diamond }^1\left(\Omega \right)$ 满足等式：

${\displaystyle {\int }_{\Omega }A\nabla u\cdot \nabla v\text{d}x}=\left(f,v\right)+\langle g,\gamma v\rangle$ (1.4)

对所有的$v\in H_{\diamond }^1\left(\Omega \right)$ ，这里的运算符$\left(\cdot ,\cdot \right)$ 和$\langle \cdot ,\cdot \rangle$ 分别表示空间$L^2\left(\Omega \right)$ 和$L^2\left(\partial \Omega \right)$ 上的内积。令$v=u$ ，结合Poincaré-Friedrichs不等式$C{\displaystyle {\int }_{\Omega }{u}^2\text{d}x}\le {\displaystyle {\int }_{\Omega }{\left|\nabla u\right|}^2\text{d}x}$ 和强制性条件${\Vert u\Vert }_{H^1\left(\Omega \right)}^2\le C{\displaystyle {\int }_{\Omega }{\left|\nabla u\right|}^2}\le C{\displaystyle {\int }_{\Omega }A\nabla u\cdot \nabla u\text{d}x}$ ，对所有$u\in H_{\diamond }^1\left(\Omega \right),A\in {\mathcal{A}}_{ad}$ 成立，即可得出结论(1.3)。

接下来定义非线性的系数到解算子$\mathcal{U}:{\mathcal{A}}_{ad}\to H_{\diamond }^1\left(\Omega \right)$ ，其将每个$\left(A,f,g\right)\in {\mathscr{H}}_{ad}$ 映射到问题(1.1)的唯一弱解$u_{A,f,g}:=\Phi$ 。为了便于计算纯Neumann问题的数值解，我们对解进行了标准化，使其在边界上的均值为零。

本研究的识别问题现在表述如下：

给定${\Phi }^{†}:={\mathcal{U}}_{A,f,g}\in H_{\diamond }^1\left(\Omega \right)$ ，寻找一个$\left(A,f,g\right)\in {\mathscr{H}}_{ad}$ ，使得${\Phi }^{†}$ 和$A,f,g$ 满足(1.4)式。

这个反问题可能有多个解。根据不适定问题的一般收敛理论[17]，我们将使用唯一最小范数解的概念，其定义为

$\left(A^{†},f^{†},g^{†}\right):=\arg \underset{\left(A,f,g\right)\in I\left({\Phi }^{†}\right)}{\min }ℛ\left(A,f,g\right)$ (1.5)

其中$I\left({\Phi }^{†}\right):=\left\{\left(A,f,g\right)\in {\mathscr{H}}_{ad}|{\mathcal{U}}_{A,f,g}={\Phi }^{†}\right\}$ ，并且

$ℛ\left(A,f,g\right):={\Vert A\Vert }_{L_{sym}^2\left(\Omega \right)}^2+{\Vert f\Vert }_{L^2\left(\Omega \right)}^2+{\Vert g\Vert }_{L^2\left(\partial \Omega \right)}^2$ .

在下面的问题(1.6)中，函数$ℛ$ 作为Tikhonov正则化方法的惩罚项出现，$\left(A^{†},f^{†},g^{†}\right)$ 是识别问题中具有最小惩罚值的解。注意到问题(1.5)的可允许集$I\left({\Phi }^{†}\right)$ 非空，并且在$L_{sym}^2\left(\Omega \right)\times L^2\left(\Omega \right)\times L^2\left(\partial \Omega \right)$ 这个函数上具有凸性和弱闭性，因此，最小值$\left(A^{†},f^{†},g^{†}\right)$ 被唯一地确定。本文使用如下的凸能量泛函

$\left(A,f,g\right)\in {\mathscr{H}}_{ad}\mapsto {\mathcal{J}}_{\delta }^h\left(A,f,g\right):={\displaystyle {\sum }_{l=1}^L{\displaystyle {\int }_{\Omega }A\nabla \left(u_{A,f,g}^{l,h}-z_{\delta }^l\right)\cdot \nabla \left(u_{A,f,g}^{l,h}-z_{\delta }^l\right)\text{d}x}}$ ,

考虑严格凸问题

$\underset{\left(A,f,g\right)\in {\mathscr{H}}_{ad}}{\min }{\mathcal{J}}_{\delta }^h\left(A,f,g\right)+\eta ℛ\left(A,f,g\right)$ (1.6)

的唯一极小化解$\left(A^h,f^h,g^h\right)$ 作为识别问题的离散正则化解。为了便于说明，我们将考虑一组测量数据${\left(z_{\delta }\right)}_{\delta >0}$ 的情形。在本文中，为了方便相关符号表示，我们用${\displaystyle {\int }_{\Omega }\cdots }$ 代替${\displaystyle {\int }_{\Omega }\cdots \text{d}x}$ 。我们使用来自例如[18]的Sobolev空间$H^1\left(\Omega \right),H^2\left(\Omega \right),W^{k,p}\left(\Omega \right)$ 等的标准概念。

注：文献[19]中证明了测量数据个数$L=4$ 可以确保扩散矩阵$A$ 的唯一反演。且在文献[20]中，数值实验表明，$L$ 越大，$A$ 的恢复越准确。所以本文在数值模拟中，选择了5个测量数据，以确保解的唯一性，平衡所提算法的效率和准确性。

3. 有限元离散

3.1. 预备知识

在乘积空间$L_{sym}^2\left(\Omega \right)\times L^2\left(\Omega \right)\times L^2\left(\partial \Omega \right)$ 和$L_{sym}^{\infty }\left(\Omega \right)\times L^2\left(\Omega \right)\times L^2 \left(\partial \Omega \right)$ 中分别使用范数

${\Vert \left(H,l,s\right)\Vert }_{L_{sym}^2\left(\Omega \right)\times L^2\left(\Omega \right)\times L^2\left(\partial \Omega \right)}={\left({\Vert H\Vert }_{L_{sym}^2\left(\Omega \right)}^2+{\Vert l\Vert }_{L_{ }^2\left(\Omega \right)}^2+{\Vert s\Vert }_{L_{ }^2\left(\partial \Omega \right)}^2\right)}^{1/2}$

和

${\Vert \left(H,l,s\right)\Vert }_{L_{sym}^{\infty }\left(\Omega \right)\times L^2\left(\Omega \right)\times L^2\left(\partial \Omega \right)}={\Vert H\Vert }_{L_{sym}^{\infty }\left(\Omega \right)}^{ }+{\Vert l\Vert }_{L_{ }^2\left(\Omega \right)}^{ }+{\Vert s\Vert }_{L_{ }^2\left(\partial \Omega \right)}^{ }.$

注意到，系数到解的算子

$\mathcal{U}:{\mathscr{H}}_{ad}\subset L_{sym}^{\infty }\left(\Omega \right)\times L^2\left(\Omega \right)\times L^2\left(\partial \Omega \right)\to H_{\diamond }^1{}_{ }\left(\Omega \right)$

和解

$\Gamma :=\left(A,f,g\right)\in {\mathscr{H}}_{ad}\to \mathcal{U}\left(\Gamma \right):={\mathcal{U}}_{\Gamma }$

在${\mathscr{H}}_{ad}$ 上是Fréchet可微的，对于每个$\Gamma =\left(A,f,g\right)\in {\mathscr{H}}_{ad}$ 和所有的$\phi \in H^1{\left(\Omega \right)}_{ }$ ，其在方向$\lambda :=\left(H,l,s\right)\in L_{sym}^{\infty }\left(\Omega \right)\times L^2\left(\Omega \right)\times L^2\left(\partial \Omega \right)$ 上的Fréchet导数${\xi }_{\lambda }:={{\mathcal{U}}^{\prime }}_{\Gamma }\left(\lambda \right):={\mathcal{U}}^{\prime }\left(\Gamma \right)\left(\lambda \right)$ 是方程

${\displaystyle {\int }_{\Omega }A\nabla {\xi }_{\lambda }\cdot \nabla \phi }={\displaystyle {\int }_{\Omega }H\nabla {\mathcal{U}}_{\text{Γ}}\cdot \nabla \phi }+\left(l,\phi \right)+\langle s,\gamma \phi \rangle$ (2.1)

在$H_{\diamond }^1\left(\Omega \right)$ 中的唯一弱解。

在$S_d$ 中，我们引入凸子集$\mathcal{K}:=\left\{M\in S_d|\underset{\_}{q}\le M\xi \cdot \xi \le \overline{q}\forall \xi \in R^d\right\}$以及正交投影：$P_{\mathcal{K}}:S_d\to \mathcal{K}$ ，则$\left(A-P_{\mathcal{K}}\left(A\right)\right)\cdot \left(B-P_{\mathcal{K}}\left(A\right)\right)\le 0$ 对所有的$A\in S_d$ 和$B\in \mathcal{K}$ 成立。此外，设$\xi :=\left({\xi }_1,\cdots ,{\xi }_d\right)$ 和$\zeta :=\left({\zeta }_1,\cdots ,{\zeta }_d\right)$ 是$R^d$ 中的任意两个向量，使用符号${\left(\xi \otimes \zeta \right)}_{1\le i,j\le d}\in S_d$ 表示

${\left(\xi \otimes \zeta \right)}_{ij}:=\frac{1}{2}\left({\xi }_i{\zeta }_j+{\xi }_j{\zeta }_i\right) \forall i,j=1,\cdots ,d.$

3.2. 离散化

设${\left({\mathcal{T}}_h\right)}_{0<h<1}$ 为网格大小为$h$ 的区域$\overline{\Omega }$ 的拟正则均匀三角剖分，使得多边形边界$\partial \Omega$ 的每个顶点都是${\mathcal{T}}_h$ 的节点。状态函数的离散空间定义为

${\mathcal{V}}_1^h:=\left\{{\phi }^h\in C\left(\overline{\Omega }\right)\cap H_{\diamond }^1\left(\Omega \right)|{\phi }^h{}_{|T}\in P_1\left(T\right)\forall T\in {\mathcal{T}}_h\right\}$, (2.2)

其中$P_r$ 表示所有次数不超过$r$ 的多项式函数。与连续情况类似，我们有如下结果。

引理2.1 设$\left(A,f,g\right)\in {\mathscr{H}}_{ad}$ ，则变分方程

${\displaystyle {\int }_{\Omega }A\nabla {\Phi }^h\cdot \nabla {\phi }^h}=\left(f,{\phi }^h\right)+\langle g,\gamma {\phi }^h\rangle$ (2.3)

对所有${\phi }^h\in {\mathcal{V}}_1^h$ 存在唯一解${\Phi }^h\in {\mathcal{V}}_1^h$ 且满足先验估计：

${\Vert {\Phi }^h\Vert }_{H^1\left(\Omega \right)}\le C_{\mathcal{N}}\left({\Vert f\Vert }_{L^2\left(\Omega \right)}+{\Vert g\Vert }_{L^2\left(\partial \Omega \right)}\right).$ (2.4)

映射${\mathcal{U}}^h:{\mathscr{H}}_{ad}\subset L_{sym}^{\infty }\left(\Omega \right)\times L^2\left(\Omega \right)\times L^2\left(\partial \Omega \right)\to {\mathcal{V}}_1^h$ ，从每个$\Gamma :=\left(A,f,g\right)\in {\mathscr{H}}_{ad}$ 到(2.3)的唯一解${\mathcal{U}}_{\Gamma }^h:={\Phi }^h$ 称为系数到离散解的算子。这个算子在集合${\mathscr{H}}_{ad}$ 上也是Fréchet可微的。对于每个$\Gamma =\left(A,f,g\right)\in {\mathscr{H}}_{ad}$ 和$\lambda :=\left(H,l,s\right)\in L_{sym}^{\infty }\left(\Omega \right)\times L^2\left(\Omega \right)\times L^2\left(\partial \Omega \right)$ ，Fréchet导数${\xi }_{\lambda }^h:={\mathcal{U}}_{\text{Γ}}^{h^{\prime }}\left(\lambda \right)$ 属于${\mathcal{V}}_1^h$ ，且对${\mathcal{V}}_1^h$ 中的所有${\phi }^h$ 满足等式

${\displaystyle {\int }_{\Omega }A\nabla {\xi }_{\lambda }^h\cdot \nabla {\phi }^h}=-{\displaystyle {\int }_{\Omega }H\nabla {\mathcal{U}}_{\text{Γ}}^h\cdot \nabla {\phi }^h}+\left(l,{\phi }^h\right)+\left(s,\gamma {\phi }^h\right).$ (2.5)

由于椭圆型问题的有限元方法的标准理论[21]，对于任意固定的$\Gamma =\left(A,f,g\right)\in {\mathscr{H}}_{ad}$ ，都有极限

$\underset{h\to 0}{\lim }{\Vert {\mathcal{U}}_{\Gamma }-{\mathcal{U}}_{\Gamma }^h\Vert }_{H^1\left(\Omega \right)}=0$ .(2.6)

设是Clément平滑插值算子，则满足以下性质

$\underset{h\to 0}{\lim }{\Vert \varphi -{\Pi }^h\varphi \Vert }_{H^k\left(\Omega \right)}=0$ 对所有$k\in \left\{0,1\right\}$ (2.7)

${\Vert \varphi -{\Pi }^h\varphi \Vert }_{H^k\left(\Omega \right)}\le C{h}^{l-k}{\Vert \varphi \Vert }_{H^1\left(\Omega \right)}$ 对于$0\le k\le l\le 2$ ， (2.8)

其中$C$ 与$h$ 和$\varphi$ 无关[22]。因此，利用离散算子${\mathcal{U}}^h$ 和插值算子${\Pi }^h$ ，可以引入离散代价函数

${\mathcal{J}}_{\delta }^h\left(A,f,g\right):={\displaystyle {\int }_{\Omega }A\nabla \left({\mathcal{U}}_{A,f,g}^h-{\Pi }^h{z}_{\delta }\right)\cdot \nabla \left({\mathcal{U}}_{A,f,g}^h-{\Pi }^h{z}_{\delta }\right)}$ (2.9)

其中$\left(A,f,g\right)\in {\mathscr{H}}_{ad}$ 。

引理2.2 假设序列${\left({\Gamma }_n\right)}_n:={\left(A_n,f_n,g_n\right)}_n\subset {\mathscr{H}}_{ad}$ 弱收敛于$L_{sym}^{\infty }\left(\Omega \right)\times L^2\left(\Omega \right)\times L^2\left(\partial \Omega \right)$ 中的$\Gamma :=\left(A,f,g\right)$ 。则对于任何固定的$h>0$ ，序列${\left({\mathcal{U}}_{{\text{Γ}}_n}^h\right)}_n\subset {\mathcal{V}}_1^h$ 在$H^1\left(\Omega \right)$ 中收敛到${\mathcal{U}}_{\text{Γ}}^h$ 。

证明：设${\left(A_n\right)}_n$ 有一个由相同符号表示的子序列，该子序列在$L_{sym}^{\infty }\left(\Omega \right)$ 中弱收敛到$A$ 。此外，通过不等式(2.4)，对应的状态序列${\left({\mathcal{U}}_{{\text{Γ}}_n}^h\right)}_n$ 在有限维空间${\mathcal{V}}_1^h$ 中是有界的。则存在一个不重新标记的子序列${\left({\mathcal{U}}_{{\text{Γ}}_n}^h\right)}_n$ 及${\Theta }^h\in {\mathcal{V}}_1^h$ ，使得${\left({\mathcal{U}}_{{\text{Γ}}_n}^h\right)}_n$ 在$H^1\left(\Omega \right)$ 中收敛到${\Theta }^h$ 。由式(2.3)可知

${\displaystyle {\int }_{\Omega }{A}_n\nabla \left({\mathcal{U}}_{{\text{Γ}}_n}^h-{\mathcal{U}}_{\text{Γ}}^h\right)\cdot \nabla {\phi }^h}={\displaystyle {\int }_{\Omega }\left(A-A_n\right)\nabla {\mathcal{U}}_{\text{Γ}}^h\cdot \nabla {\phi }^h}+\left(f_n-f,{\phi }^h\right)+\langle g_n-g,\gamma {\phi }^h\rangle$ (2.10)

对于所有${\phi }^h\in {\mathcal{V}}_1^h$ 。取${\phi }^h={\mathcal{U}}_{{\Gamma }_n}^h-{\mathcal{U}}_{\Gamma }^h$ ，由(1.5)得到

$\begin{array}{c}C{\Vert {\mathcal{U}}_{{\Gamma }_n}^h-{\mathcal{U}}_{\Gamma }^h\Vert }_{H^1\left(\Omega \right)}^2\le {\displaystyle {\int }_{\Omega }\left(A-A_n\right)\nabla {\mathcal{U}}_{\Gamma }^h\cdot \nabla \left({\mathcal{U}}_{{\Gamma }_n}^h-{\Theta }^h+{\Theta }^h-{\mathcal{U}}_{\Gamma }^h\right)} \\ +\left(f_n-f,{\mathcal{U}}_{{\Gamma }_n}^h-{\Theta }^h+{\Theta }^h-{\mathcal{U}}_{\Gamma }^h\right)\\ +\langle g_n-g,\gamma \left({\mathcal{U}}_{{\Gamma }_n}^h-{\Theta }^h+{\Theta }^h-{\mathcal{U}}_{\Gamma }^h\right)\rangle \end{array}$

$\begin{array}{c}\le C{\Vert {\mathcal{U}}_{{\Gamma }_n}^h-{\Theta }^h\Vert }_{H^1\left(\Omega \right)}+{\displaystyle {\int }_{\Omega }\left(A-A_n\right)\nabla {\mathcal{U}}_{\Gamma }^h\cdot \nabla \left({\Theta }^h-{\mathcal{U}}_{\Gamma }^h\right)}\\ +\left(f_n-f,{\Theta }^h-{\mathcal{U}}_{\Gamma }^h\right)+\langle g_n-g,\gamma \left({\Theta }^h-{\mathcal{U}}_{\Gamma }^h\right)\rangle \end{array}$ (2.11)

由于$L_{sym}^{\infty }\left(\Omega \right)$ 中的$A_n\rightharpoonup A$ ，我们得到$\underset{n\to \infty }{\lim }{\displaystyle {\int }_{\Omega }\left(A-A_n\right)\nabla {\mathcal{U}}_{\Gamma }^h\cdot \nabla \left({\Theta }^h-{\mathcal{U}}_{\Gamma }^h\right)}=0$ 。因此让$n$ 趋于$\infty$ ，从(2.11)中得到，$\underset{n\to \infty }{\lim }{\Vert {\mathcal{U}}_{{\Gamma }_n}^h-{\mathcal{U}}_{\Gamma }^h\Vert }_{H^1\left(\Omega \right)}=0$ ，完成了证明。

现陈述以下关于代价函数凸性的结果。

引理2.3 ${\mathcal{J}}_{\delta }^h$ 在集合${\mathscr{H}}_{ad}$ 上关于范数$L_{sym}^2\left(\Omega \right)\times L^2\left(\Omega \right)\times L^2\left(\partial \Omega \right)$ 是凸且连续的。

证明：函数${\mathcal{J}}_{\delta }^h$ 的连续性直接由引理2.2得出。接下来证明${\mathcal{J}}_{\delta }^h$ 是凸的。

令$\Gamma :=\left(A,f,g\right)\in {\mathscr{H}}_{ad}$ 和$\lambda :=\left(H,l,s\right)\in L_{sym}^{\infty }\left(\Omega \right)\times L^2\left(\Omega \right)\times L^2\left(\partial \Omega \right)$ 得到

${\mathcal{U}}_{\Gamma }^h{}^{\prime }\left(\lambda \right)=\frac{\partial {\mathcal{U}}_{\Gamma }^h}{\partial A}H+\frac{\partial {\mathcal{U}}_{\Gamma }^h}{\partial f}l+\frac{\partial {\mathcal{U}}_{\Gamma }^h}{\partial g}s$ 和${\mathcal{J}}_{\delta }^h{}^{\prime }\left(\Gamma \right)\left(\lambda \right)=\frac{\partial {\mathcal{J}}_{\delta }^h\left(\Gamma \right)}{\partial A}H+\frac{\partial {\mathcal{J}}_{\delta }^h\left(\Gamma \right)}{\partial f}l+\frac{\partial {\mathcal{J}}_{\delta }^h\left(\Gamma \right)}{\partial g}s.$

分别计算最后一个方程右侧的每一项可得，

$\begin{array}{c}{\mathcal{J}}_{\delta }^h{}^{\prime }\left(\Gamma \right)\left(\lambda \right)={\displaystyle {\int }_{\Omega }H\nabla }\left({\mathcal{U}}_{\Gamma }^h-{\Pi }_{z_{\delta }}^h\right)\cdot \nabla \left({\mathcal{U}}_{\Gamma }^h-{\Pi }^h{}_{z_{\delta }}\right)+2{\displaystyle {\int }_{\Omega }A\nabla }\left(\frac{\partial {\mathcal{U}}_{\Gamma }^h}{\partial A}H\right)\cdot \nabla \left({\mathcal{U}}_{\Gamma }^h-{\Pi }_{z_{\delta }}^h\right)\\ +2{\displaystyle {\int }_{\Omega }A\nabla }\left(\frac{\partial {\mathcal{U}}_{\Gamma }^h}{\partial f}l\right)\cdot \nabla \left({\mathcal{U}}_{\Gamma }^h-{\Pi }_{z_{\delta }}^h\right)+2{\displaystyle {\int }_{\Omega }A\nabla }\left(\frac{\partial {\mathcal{U}}_{\Gamma }^h}{\partial g}s\right)\cdot \nabla \left({\mathcal{U}}_{\Gamma }^h-{\Pi }_{z_{\delta }}^h\right)\\ =2{\displaystyle {\int }_{\Omega }A\nabla }\left(\frac{\partial {\mathcal{U}}_{\Gamma }^h}{\partial Q}H+\frac{\partial {\mathcal{U}}_{\Gamma }^h}{\partial f}l+\frac{\partial {\mathcal{U}}_{\Gamma }^h}{\partial g}s\right)\cdot \nabla \left({\mathcal{U}}_{\Gamma }^h-{\Pi }_{z_{\delta }}^h\right)+{\displaystyle {\int }_{\Omega }H\nabla }\left({\mathcal{U}}_{\Gamma }^h-{\Pi }_{z_{\delta }}^h\right)\cdot \nabla \left({\mathcal{U}}_{\Gamma }^h-{\Pi }_{z_{\delta }}^h\right)\\ =2{\displaystyle {\int }_{\Omega }A\nabla {\mathcal{U}}_{\Gamma }^h{}^{\prime }}\left(\lambda \right)\cdot \nabla \left({\mathcal{U}}_{\Gamma }^h-{\Pi }_{z_{\delta }}^h\right)+{\displaystyle {\int }_{\Omega }H\nabla }\left({\mathcal{U}}_{\Gamma }^h-{\Pi }_{z_{\delta }}^h\right)\cdot \nabla \left({\mathcal{U}}_{\Gamma }^h-{\Pi }_{z_{\delta }}^h\right)\\ =2{\displaystyle {\int }_{\Omega }A\nabla {\mathcal{U}}_{\Gamma }^h{}^{\prime }}\left(\lambda \right)\cdot \nabla \left({\mathcal{U}}_{\Gamma }^h-{\overline{\Pi }}_{z_{\delta }}^h\right)+{\displaystyle {\int }_{\Omega }H\nabla }\left({\mathcal{U}}_{\Gamma }^h-{\overline{\Pi }}_{z_{\delta }}^h\right)\cdot \nabla \left({\mathcal{U}}_{\Gamma }^h-{\overline{\Pi }}_{z_{\delta }}^h\right),\end{array}$

其中，

${\overline{\Pi }}_{z_{\delta }}^h:={\Pi }_{z_{\delta }}^h-{\left|\Omega \right|}^{-1}\langle 1,\gamma {\Pi }_{z_{\delta }}^h\rangle \in {\mathcal{V}}_1^h$且$\nabla {\overline{\Pi }}_{z_{\delta }}^h=\nabla {\Pi }_{z_{\delta }}^h.$ (2.12)

根据(2.5)，可以推断出

$\begin{array}{c}{\mathcal{J}}_{\delta }^h{}^{\prime }\left(\Gamma \right)\left(\lambda \right)=-2{\displaystyle {\int }_{\Omega }H\nabla {\mathcal{U}}_{\Gamma }^h\cdot \nabla \left({\mathcal{U}}_{\Gamma }^h-{\overline{\Pi }}_{z_{\delta }}^h\right)}+2\left(l,{\mathcal{U}}_{\Gamma }^h-{\overline{\Pi }}_{z_{\delta }}^h\right)\\ +2\langle s,\gamma \left({\mathcal{U}}_{\Gamma }^h-{\overline{\Pi }}_{z_{\delta }}^h\right)\rangle +{\displaystyle {\int }_{\Omega }H\nabla }\left({\mathcal{U}}_{\Gamma }^h-{\overline{\Pi }}_{z_{\delta }}^h\right)\cdot \nabla \left({\mathcal{U}}_{\Gamma }^h-{\overline{\Pi }}_{z_{\delta }}^h\right)\\ =-{\displaystyle {\int }_{\Omega }H\nabla {\mathcal{U}}_{\Gamma }^h\cdot \nabla {\mathcal{U}}_{\Gamma }^h}+{\displaystyle {\int }_{\Omega }H\nabla {\overline{\Pi }}_{z_{\delta }}^h\cdot \nabla {\overline{\Pi }}_{z_{\delta }}^h}\\ +2\left(l,{\mathcal{U}}_{\Gamma }^h-{\overline{\Pi }}_{z_{\delta }}^h\right)+2\langle s,\gamma \left({\mathcal{U}}_{\Gamma }^h-{\overline{\Pi }}_{z_{\delta }}^h\right)\rangle .\end{array}$ (2.13)

再利用(2.5)得出结论

$\begin{array}{c}{\mathcal{J}}_{\delta }^h{}^{\prime \text{}\prime }\left(\Gamma \right)\left(\lambda ,\lambda \right)=-2{\displaystyle {\int }_{\Omega }H\nabla {\mathcal{U}}_{\Gamma }^h\cdot \nabla {\mathcal{U}}_{\Gamma }^h{}^{\prime }\left(\lambda \right)}+2\left(l,{\mathcal{U}}_{\Gamma }^h{}^{\prime }\left(\lambda \right)\right)+2\langle s,\gamma {\mathcal{U}}_{\Gamma }^h{}^{\prime }\left(\lambda \right)\rangle \\ =2{\displaystyle {\int }_{\Omega }A\nabla {\mathcal{U}}_{\Gamma }^h{}^{\prime }\left(\lambda \right)\cdot \nabla {\mathcal{U}}_{\Gamma }^h{}^{\prime }\left(\lambda \right)}\ge 2C{\Vert {\mathcal{U}}_{\Gamma }^h{}^{\prime }\left(\lambda \right)\Vert }_{H^1\left(\Omega \right)}^2\ge 0,\end{array}$

结合(1.5)可得${\mathcal{J}}_{\delta }^h$ 在集合${\mathscr{H}}_{ad}$ 上关于范数$L_{sym}^2\left(\Omega \right)\times L^2\left(\Omega \right)\times L^2\left(\partial \Omega \right)$ 是凸的，证毕。

现在我们证明本节的重要结论。

定理2.4 严格凸优化问题

$\underset{\left(A,f,g\right)\in {\mathscr{H}}_{ad}}{\min }{F}_{\eta }\left(A,f,g\right):={\mathcal{J}}_{\delta }^h\left(A,f,g\right)+\eta ℛ\left(A,f,g\right)$

存在唯一的极小化解。此外，$\Gamma :=\left(A,f,g\right)\in {\mathscr{H}}_{ad}$ 是的唯一极小化解，当且仅当系统

$A\left(x\right)=P_{\mathcal{K}}\left(\frac{1}{2\rho }\left(\nabla {\mathcal{U}}_{\Gamma }^h\left(x\right)\otimes \nabla {\mathcal{U}}_{\Gamma }^h\left(x\right)-\nabla {\overline{\Pi }}_{z_{\delta }}^h\left(x\right)\otimes \nabla {\overline{\Pi }}_{z_{\delta }}^h\left(x\right)\right)\right),$ (2.14)

$f\left(x\right)=\frac{1}{\rho }\left({\overline{\Pi }}_{z_{\delta }}^h\left(x\right)-{\mathcal{U}}_{\Gamma }^h\left(x\right)\right),$ (2.15)

$g\left(x\right)=\frac{1}{\rho }\Upsilon \left({\overline{\Pi }}_{z_{\delta }}^h\left(x\right)-{\mathcal{U}}_{\Gamma }^h\left(x\right)\right)$ (2.16)

在$\Omega$ 中几乎处处成立，其中${\overline{\Pi }}^h$ 是根据(2.12)从${\Pi }^h$ 生成的。

证明：令${\left({\Gamma }_n\right)}_n:={\left(A_n,f_n,g_n\right)}_n\subset {\mathscr{H}}_{ad}$ 为的一个极小化序列，即$\underset{n\to \infty }{\lim }{\Upsilon }_{\delta }^{\rho ,h}\left({\Gamma }_n\right)=\underset{\left(A,f,g\right)\in {\mathscr{H}}_{ad}}{\inf }{\Upsilon }_{\delta }^{\rho ,h}\left(A,f,g\right)$ 。因此，序列${\left({\Gamma }_n\right)}_n$ 在范数$L_{sym}^2\left(\Omega \right)\times L^2\left(\Omega \right)\times L^2\left(\partial \Omega \right)$ 下是有界的。存在一个子序列(不重新标记)和$\Gamma :=\left(A,f,g\right)\in L_{sym}^2\left(\Omega \right)\times L^2\left(\Omega \right)\times L^2\left(\partial \Omega \right)$ ，使得${\Gamma }_n$ 在$L_{sym}^2\left(\Omega \right)\times L^2\left(\Omega \right)\times L^2\left(\partial \Omega \right)$ 中弱收敛于$\Gamma$ 。另一方面，由于${\mathscr{H}}_{ad}$ 是$L_{sym}^2\left(\Omega \right)\times L^2\left(\Omega \right)\times L^2\left(\partial \Omega \right)$ 中的一个闭凸子集，因此它是弱闭的，这意味着$\Gamma \in {\mathscr{H}}_{ad}$ 。根据引理2.3，${\mathcal{J}}_{\delta }^h$ 和$ℛ$ 在${\mathscr{H}}_{ad}$ 上都是弱下半连续的，因此

${\mathcal{J}}_{\delta }^h\left(\Gamma \right)\le \underset{n\to \infty }{\lim \inf }{\mathcal{J}}_{\delta }^h\left({\Gamma }_n\right)$ 且$ℛ\left(\Gamma \right)\le \underset{n\to \infty }{\lim \inf }ℛ\left({\Gamma }_n\right).$

则

$\begin{array}{c}{\mathcal{J}}_{\delta }^h\left(\Gamma \right)+\eta ℛ\left(\Gamma \right)\le \underset{n\to \infty }{\lim \inf }{\mathcal{J}}_{\delta }^h\left({\Gamma }_n\right)+\underset{n\to \infty }{\lim \inf }\eta ℛ\left({\Gamma }_n\right)\le \underset{n\to \infty }{\lim \inf }\left({\mathcal{J}}_{\delta }^h\left({\Gamma }_n\right)+\eta ℛ\left({\Gamma }_n\right)\right)\\ =\underset{n\to \infty }{\lim }{\Upsilon }_{\delta }^{\rho ,h}\left({\Gamma }_n\right)=\underset{\left(A,f,g\right)\in {\mathscr{H}}_{ad}}{\inf }{\Upsilon }_{\delta }^{\rho ,h}\left(A,f,g\right)\end{array}$

故$\Gamma$ 是的一个极小化解。由于$F_{\eta }$ 是严格凸的，因此这个极小化解是唯一的。另外，元素$\Gamma :=\left(A,f,g\right)\in {\mathscr{H}}_{ad}$ 是的极小化解当且仅当对于所有$\overline{\Gamma }=\left(H,l,s\right)\in {\mathscr{H}}_{ad}$ ，${\Upsilon }_{\delta }^{\rho ,h}{}^{\prime }\left(\Gamma \right)\left(\overline{\Gamma }-\Gamma \right)\ge 0$ 。则根据表达式(2.13)，得到

$\begin{array}{c}0\le {\displaystyle {\int }_{\Omega }\left(H-A\right)\nabla {\overline{\Pi }}_{z_{\delta }}^h\cdot \nabla {\overline{\Pi }}_{z_{\delta }}^h}-{\displaystyle {\int }_{\Omega }\left(H-A\right)\nabla {\mathcal{U}}_{\Gamma }^h\cdot \nabla {\mathcal{U}}_{\Gamma }^h}+2\rho \left(H-A,A\right)\\ +2\left(l-f,{\mathcal{U}}_{\Gamma }^h-{\overline{\Pi }}_{z_{\delta }}^h\right)+2\rho \left(l-f,f\right)+2\langle s-g,\gamma \left({\mathcal{U}}_{\Gamma }^h-{\overline{\Pi }}_{z_{\delta }}^h\right)\rangle +2\rho \langle s-g,g\rangle \\ ={\displaystyle {\int }_{\Omega }\left(H-A\right)\cdot \left(\nabla {\overline{\Pi }}_{z_{\delta }}^h\otimes \nabla {\overline{\Pi }}_{z_{\delta }}^h-\nabla {\mathcal{U}}_{\Gamma }^h\otimes \nabla {\mathcal{U}}_{\Gamma }^h+2\rho A\right)}\\ +2\left(l-f,{\mathcal{U}}_{\Gamma }^h-{\overline{\Pi }}_{z_{\delta }}^h+\rho f\right)+2\langle s-g,\gamma \left({\mathcal{U}}_{\Gamma }^h-{\overline{\Pi }}_{z_{\delta }}^h\right)+\rho g\rangle \end{array}$

对于所有$\overline{\Gamma }=\left(H,l,s\right)\in {\mathscr{H}}_{ad}$ 成立。将${\overline{\Gamma }}_1=\left(H,f,g\right)$ ，${\overline{\Gamma }}_2=\left(A,l,g\right)$ 和${\overline{\Gamma }}_3=\left(A,f,s\right)$ 代入上述不等式即可得到方程(2.14)~(2.16)。证毕。

4. 数值实验

4.1. 投影牛顿算法

根据定理2.4中的必要最优性系统，我们提出了一种投影牛顿算法。第2节中的论证表明，目标函数$F_{\eta }\left(A_h,f_h,g_h^l\right)$ 的方向导数${F^{\prime }}_{\eta }\left(A_h,f_h,g_h^l\right)\left[D,t,s\right]$ 可以表示为

$\begin{array}{l}{F^{\prime }}_{\eta }\left(A_h,f_h,g_h^l\right)\left[D,t,s\right]={\left(D\nabla z_{\delta ,h}^l,\nabla z_{\delta ,h}^l\right)}_{L^2{\left(\Omega \right)}^d}-{\left(D\nabla u_h^l,\nabla u_h^l\right)}_{L^2{\left(\Omega \right)}^d}+2\ast {\left(t,u_h^l-z_{\delta ,h}^l+\eta f_h\right)}_{L^2\left(\Omega \right)}\\ +2\ast {\langle s,u_h^l-z_{\delta ,h}^l+\eta g_h^l\rangle }_{L^2\left(\partial \Omega \right)}+2\ast \eta {\left(D,A_h\right)}_{L^2{\left(\Omega \right)}^{d,d}},\forall \left(D,t,s\right)\in {\mathscr{H}}_h,\end{array}$

其中$u_h^l$ 满足如下方程

${\left(A_h\nabla u_h^l,\nabla v_h\right)}_{L^2{\left(\Omega \right)}^d}-{\left(f_h,v_h\right)}_{L^2\left(\Omega \right)}-{\langle g_h^l,v_h\rangle }_{L^2\left(\partial \Omega \right)}=0,\forall v_h\in V_h,l=\overline{1,\cdots ,L}$ (3.1)

逐点的最优性条件为：

$A_h\left(x\right)=P_K\left(-\frac{1}{2\eta }\left(\nabla z_{\delta ,h}^l\otimes \nabla z_{\delta ,h}^l-\nabla u_h^l\otimes \nabla u_h^l\right)\right)$ (3.2)

在$\Omega$ 中几乎处处成立。

总之，变分形式下的一阶必要最优性系统为

$\{\begin{array}{l}{\left(A_h\nabla u_h^l,\nabla v_h\right)}_{L^2{\left(\Omega \right)}^d}-{\left(f_h,v_h\right)}_{L^2\left(\Omega \right)}-{\langle g_h^l,v_h\rangle }_{L^2\left(\partial \Omega \right)}=0,\forall v_h\in V_h,l=\overline{1,\cdots ,L},\\ {\left(\left(D_h-A_h\right)\nabla z_{\delta ,h}^l,\nabla z_{\delta ,h}^l\right)}_{L^2{\left(\Omega \right)}^d}-{\left(\left(D_h-A_h\right)\nabla u_h^l,\nabla u_h^l\right)}_{L^2{\left(\Omega \right)}^d}+2\ast {\left(t_h-f_h,u_h^l-z_{\delta ,h}^l+\eta f_h\right)}_{L^2\left(\Omega \right)}\\ +2\ast \eta {\left(A_h,D_h-A_h\right)}_{L^2{\left(\Omega \right)}^{d,d}}+2\ast {\langle s_h-g_h^l,u_h^l-z_{\delta ,h}^l+\eta g_h^l\rangle }_{L^2\left(\partial \Omega \right)}\ge 0,\forall \left(D_h,t_h,s_h\right)\in {\mathscr{H}}_h.\end{array}$

为了应用牛顿法，一个关键的步骤是推导未知数$\left(u^l,f,A\right)$ 的牛顿更新$\left({\overline{u}}^l,\overline{f},\overline{A}\right)$ 。这可以通过求解以下方程来实现

$\{\begin{array}{l}{\left(\overline{A}\nabla u_h^{l,n},\nabla v_h\right)}_{L^2{\left(\Omega \right)}^d}+{\left(A_h\nabla {\overline{u}}^l,\nabla v_h\right)}_{L^2{\left(\Omega \right)}^d}-{\left(\overline{f},v_h\right)}_{L^2\left(\Omega \right)}\\ ={\left(A_h^n\nabla u_h^{l,n},\nabla v_h\right)}_{L^2{\left(\Omega \right)}^d}-{\left(f_h^n,v_h\right)}_{L^2\left(\Omega \right)},\forall v_h\in V_h,l=\overline{1,\cdots ,L},\\ 2\ast \eta {\left(\overline{A},D_h\right)}_{L^2{\left(\Omega \right)}^{d,d}}-2\ast {\left(D_h\nabla {\overline{u}}^l,\nabla u_h^{l,n}\right)}_{L^2{\left(\Omega \right)}^d}+2\ast {\left(t_h,{\overline{u}}^l\right)}_{L^2\left(\Omega \right)}+2\ast \eta {\left(t_h,\overline{f}\right)}_{L^2\left(\Omega \right)}\\ ={\left(D_h\nabla z_{\delta ,h}^l,\nabla z_{\delta ,h}^l\right)}_{L^2{\left(\Omega \right)}^d}-{\left(D_h\nabla u_h^{l,n},\nabla u_h^{l,n}\right)}_{L^2{\left(\Omega \right)}^d}+2\ast {\left(t_h,u_h^{l,n}-z_{\delta ,h}^l+\eta f_h^n\right)}_{L^2\left(\Omega \right)}\\ +2\ast \eta {\left(A_h^n,D_h\right)}_{L^2{\left(\Omega \right)}^{d,d}},\forall \left(D_h,t_h\right)\in {\mathcal{A}}_h\times {ℱ}_h,l=\overline{1,\cdots ,L}.\end{array}$(3.3)

耦合系统(3.3)的解记为$\left({\overline{u}}_h^l,{\overline{f}}_h,{\overline{A}}_h\right)$ ，其中$l=\overline{1,\cdots ,L}$ 。在为系统(3.3)制定牛顿法时，我们不直接处理变分不等式(即，导电率$A_h$ 的特征值的区间约束)，即在(3.2)中的投影算子。因此，我们使用牛顿迭代来更新方程(3.1)和(3.2)的解，而不进行逐点投影$P_K$ 。在更新后，将更新的解投影到凸集$K$ 上。这导致了一个投影牛顿算法，其细节详见算法1。请注意，该方法的每次迭代需要求解一个耦合线性系统，包括状态变量${\overline{u}}_h^l$ 、源项变量${\overline{f}}_h$ 和导电率${\overline{A}}_h$ 。第6行的停止条件是相对误差

$\text{max}\left\{\frac{{\Vert u_h^{l,n+1}-u_h^{l,n}\Vert }_2}{{\Vert u_h^{l,n+1}\Vert }_2},\frac{{\Vert f_h^{l,n+1}-f_h^{l,n}\Vert }_2}{{\Vert f_h^{l,n+1}\Vert }_2},\frac{{\Vert A_h^{n+1}-A_h^n\Vert }_2}{{\Vert A_h^{n+1}\Vert }_2},l=1,\cdots ,L\right\}$

降到给定的容忍度$\tau$ 以下，或者最大迭代次数超过预定数。

众所周知，牛顿型算法在提供良好的初始猜测时收敛非常快。正则化参数$\eta$ 的选择对于获得满意的重建结果至关重要。为了选择合适的正则化参数并为牛顿算法提供良好的初始猜测，我们采用了一种路径跟随策略，这在许多应用中取得了成功[23]。具体而言，设定$\rho \in \left(0,1\right)$ ，然后我们应用算法1，令${\eta }_n={\eta }_0\ast {\rho }^n$ ，初始猜测由${\eta }_{n-1}$ 问题的解给出，即：

$\underset{\left(A,f,g\right)\in {\mathscr{H}}_{ad}}{\min }\left\{F_{{\eta }_{n-1}}\left(A,f,g\right)={\displaystyle {\sum }_{l=1}^L{\displaystyle {\int }_{\Omega }A\nabla }}\left(u_{A,f,g}^{l,h}-z_{\delta }^l\right)\cdot \nabla \left(u_{A,f,g}^{l,h}-z_{\delta }^l\right)\text{d}x+\eta \left({\Vert A\Vert }_{L_{sym}^2\left(\Omega \right)}^2+{\Vert f\Vert }_{L^2\left(\Omega \right)}^2+{\Vert g\Vert }_{L^2\left(\partial \Omega \right)}^2\right)\right\}$

这将为牛顿更新提供热启动，并充分利用了牛顿型方法的快速收敛性。最终的正则化参数是通过偏差原则确定的[24] [25]。

算法1. 投影牛顿法

1: 给定初始猜测$\left(A^0,f^0,g^0\right)$ 和容忍度$\tau$ 。 2:for $n=0,1,2,\cdots$ do 3: 求解牛顿系统(3.3)获得$u_h^{l,n},f_h^n,A_h^n$ 的增量${\overline{u}}_h^{l,n},{\overline{f}}_h^n,{\overline{A}}_h^n$ 。 4: 更新$u_h^{l,n},f_h^n,A_h^n$ 通过 $u_h^{l,n+1}=u_h^{l,n}-{\overline{u}}_h^{l,n}$ ，$f_h^{n+1}=f_h^n-{\overline{f}}_h^n$ ，$A_h^{n+1}=A_h^n-{\overline{A}}_h^n$ ，$l=1,\cdots ,L$ 。 5: 将$A_h^{n+1}$ 投影到${\mathcal{A}}_h$ 中。 6: 检查停机准则。 7: end for

总之，该算法包含两个循环。由于内部牛顿迭代的快速局部收敛性以及来自路径跟随策略的适当初始猜测，通常只需一到两个内部迭代即可确保收敛，因此，算法1预计将非常有效，这在下一部分的数值实验中得到了验证。

注：在实际应用中，面对未知真值的数据，可以采用启发式偏差原则(参考文献([26], Rule 2.1))，该原则仅依赖于测量数据就可以获得最终的正则化参数和反演结果。

4.2. 数值模拟

现在展示二维数值实验，以验证算法的准确性和效率。所有实验均在个人笔记本电脑上使用FreeFem++ [27]进行。所考虑的区域$\Omega$ 取为$\left(-1,1\right)\times \left(-1,1\right)$ 。使用一个$\left(N+1\right)\times \left(N+1\right)$ 的均匀方形网格，其中$N=20$ ，网格大小为$h=2/N$ 。用于反演$\left(A,f,g\right)$ 的数据$z^l$ 是通过在更细的网格上求解问题来模拟的。带噪声的数据$z_{\delta }^l$ 是通过对精确数据$z^l$ 进行扰动生成的，即$z_{\delta }^l\left(x\right)=z^l\left(x\right)\left(1+\delta \xi \left(x\right)\right)$ ，其中$\xi$ 是$\left[-1,1\right]$ 上的正态分布，$\delta >0$ 表示相对噪声水平。在数值实验中，使用高阶Galerkin逼近，即将$V_h\left(T\right)$ 选择为$P_2$ 有限元空间，而$X_h\left(T\right)$ 选择为$P_1$ 有限元空间。数值结果表明，这种选择显著优于低阶有限元方法。这里考虑具有对角或非对角各向异性导电率的例子，数据对应于五个不同的Neumann边界条件$\left(g^1,g^2,g^3,g^4,g^5\right):=\left(x_1^2-x_2^2,x_1-x_2,x_2-0.1\ast x_2^3,x_1+x_2,x_1\right)$ ，使用$A^{†}$ (其元素为$A_{i,j},\left(i,j=1,2\right)$ )来表示精确的导电率张量。测试实例及其不同特征列在表1中。为了对算法进行热启动，我们首先设定${\eta }_0=1$ ，并在每次迭代后将其值减少一个因子$\rho =0.6$ 。重建的导电率$\tilde{A}$ 与精确导电率$A^{†}$ 的准确性通过相对误差$e=\frac{{\Vert \tilde{A}-A^{†}\Vert }_2}{{\Vert A^{†}\Vert }_2}$ 来衡量。

Table 1. The test cases for the conductivity tensor and source term

表1. 传导张量和源项的测试情况

例	$A_{11}$	$A_{12}$	$A_{22}$	f
1	$2+\sin \left(\pi x_1\right)\ast \sin \left(\pi x_2\right)$	$0.5\ast \sin \left(2\pi x_1\right)$	$2+x_1^2+x_2^2$	$x_1+x_2$
2	$2+\left|x_1\right|+\left|x_2\right|$	$\left|x_1\ast x_2\right|$	$2+\left|\sin \left(\pi x_1\right)\ast \sin \left(\pi x_2\right)\right|$	$x_1+x_2$
3	$1+1\ast {\chi }_{S_1}$	$1\ast {\chi }_{S_2}$	$2+1\ast {\chi }_{S_3}$	$-0.5+2\ast {\chi }_{S_4}$

图1~3中展示了精确导电率张量、对于精确数据和带噪声数据的重建张量，以及一个横截面。这些图清晰地显示出，即使在不连续情况下，重建结果也相当准确。实际上，图1~3也显示出我们的算法能稳定且准确地恢复具有不同特征的各向异性导电率张量和源项，无论数据是精确的还是带有0.1%噪声的。

Figure 1. The reconstructions for the conductivity tensor and source term for Example 1: (a)~(d), (e)~(h), (i)~(l) and (m)~(p) refer to $A_{11}$ , $A_{12}$ , $A_{22}$ and $f$ , respectively. (a), (e), (i) and (m) are for exact conductivity and source term; (b), (f), (j) and (n) for reconstruction with exact data; (c), (g), (k) and (o) for noisy data with $\delta =0.1\%$ ; and (d), (h), (l) and (p) for the cross section along {$x_2=0$ }

图1. 例1中导电率张量和源项的重建结果：(a)~(d)、(e)~(h)、(i)~(l)和(m)~(p)分别对应于函数$A_{11}$ 、$A_{12}$ 、$A_{22}$ 和源项$f$ 。(a)、(e)、(i)和(m)为精确的导电率张量和源项；(b)、(f)、(j)和(n)为使用精确数据得到的重建结果；(c)、(g)、(k)和(o)为0.1%噪声下的重建；(d)、(h)、(l)和(p)为沿$x_2=0$ 的横截面

Figure 2. The reconstructions for the conductivity tensor and source term for Example 2: (a)~(d), (e)~(h), (i)~(l) and (m)~(p) refer to $A_{11}$ , $A_{12}$ , $A_{22}$ and $f$ , respectively. (a), (e), (i) and (m) are for exact conductivity and source term; (b), (f), (j) and (n) for reconstruction with exact data; (c), (g), (k) and (o) for noisy data with $\delta =0.1\%$ ; and (d), (h), (l) and (p) for the cross section along {$x_2=0$ }

图2. 例2中导电率张量和源项的重建结果：(a)~(d)、(e)~(h)、(i)~(l)和(m)~(p)分别对应于函数$A_{11}$ 、$A_{12}$ 、$A_{22}$ 和源项$f$ 。(a)、(e)、(i)和(m)为精确的导电率张量和源项；(b)、(f)、(j)和(n)为使用精确数据得到的重建结果；(c)、(g)、(k)和(o)为0.1%噪声下的重建；(d)、(h)、(l)和(p)为沿$x_2=0$ 的横截面

接下来简要检查算法的收敛性。图4展示了例3的相对误差和残差随迭代次数的变化。观察到，当噪声水平$\delta$ 从1%降至0%时，重构解的准确性稳步提高。同时从图4中可以看出，算法具有良好的收敛性：无论噪声水平$\delta$ 如何，约40次迭代后即收敛。这些结果验证了投影牛顿算法在导电率张量恢复中的效率。虽然未进一步展示，但这些观察结果在其他例子中同样适用。

Figure 3. The reconstructions for the conductivity tensor and source term for Example 3: (a)~(d), (e)~(h), (i)~(l) and (m)~(p) refer to $A_{11}$ , $A_{12}$ , $A_{22}$ and $f$ , respectively. (a), (e), (i) and (m) are for exact conductivity and source term; (b), (f), (j) and (n) for reconstruction with exact data; (c), (g), (k) and (o) for noisy data with $\delta =0.1\%$ ; and (d), (h), (l) and (p) for the cross section along {$x_2=0$ }

图3. 例3中导电率张量和源项的重建结果：(a)~(d)、(e)~(h)、(i)~(l)和(m)~(p)分别对应于函数$A_{11}$ 、$A_{12}$ 、$A_{22}$ 和源项$f$ 。(a)、(e)、(i)和(m)为精确的导电率张量和源项；(b)、(f)、(j)和(n)为使用精确数据得到的重建结果；(c)、(g)、(k)和(o)为0.1%噪声下的重建；(d)、(h)、(l)和(p)为沿$x_2=0$ 的横截面

Figure 4. The convergence of the algorithm for Example 3 with noisy data. (a) represents the relative error in identifying the conductivity tensor $A$ ; (b) represents the residual error in identifying the conductivity tensor $A$ ; and (c) represents the relative error of the source term

图4. 带有噪声数据的例3中算法的收敛性。(a) 是识别导电率张量$A$ 的相对误差；(b) 是识别导电率张量$A$ 的残差；(c) 是源项$f$ 的相对误差

5. 结论

本文针对椭圆方程的多参数识别问题提出了一种牛顿型算法。基于能量代价泛函和罚项构建了Tikhonov正则化模型，并推导出了该模型的一阶最优性系统。随后，采用投影牛顿算法在有限元空间中求解该最优性系统，实现了对椭圆方程中导电张量、源项和状态函数的同时反演。最后通过实例验证了算法的有效性和收敛性。

基金项目

江苏省大学生创新训练项目(202413573112Y)，金陵科技学院科教融合项目(2024KJRH41)，国家自然科学基金青年项目(No. 12101276)。

NOTES

^*通讯作者。