非线性不适定问题数值求解及在近场光学问题中的应用
|
摘要
本文主要讨论了非线性不适定问题数值求解及在近场光学问题中的应用.第一章是绪论,介绍了散射问题和反散射问题的数学模型和常见方法,近场光学的模型,以及不适定问题的概念和正则化方法.第二章主要介绍谱截断矩量法的相关理论及其在处理严重不适定问题中的应用,特别是在近场光学问题中的应用.通过数值计算可以看到,相比常规的正则化方法,比如Tikhonov正则化,矩量法等,谱截断矩量法处理严重不适定问题更为有效.第三章主要介绍无穷维非线性不适定问题的同伦方法.使用同伦方法主要需要解决两个问题：一是同伦曲线的存在性;二是同伦曲线的数值跟踪方法.该章主要讨论了无穷维不适定问题的两种同伦方法,分别是带Tikhonov正则化的同伦和无导数同伦,并给出了同伦曲线的几种跟踪方法.在跟踪过程中,为提高计算效率,还提出了自适应跟踪技巧.
1.近场光学的数学模型
近场光学是上世纪80年代以来出现的一个新兴学科,突破了传统光学的分辨率极限.近场光学显微镜有三种主要的模型：近场扫描光学显微镜(NSOM)[23],全内反射显微镜(TIRM)[24,25],以及光子扫描隧道显微镜(PSTM)[26,27].
NSOM有两种主要类型：照射模式(illumination mode) NSOM和接受模式(collection mode) NSOM照射模式NSOM中具有细小的光学探针,其尖端的孔径远小于光的波长,探针作为光源在被观察物体(散射体)的近场区域扫描,散射场将被探测并记录.接受模式NSOM的被观察物体由在远场的光源照射,其探针在被观察物体的近场区域探测总场(或散射场).
我们知道,当光线以大于临界角度射入一个棱镜时,光线会被发生全内反射现象,这时将产生倏逝波.全内反射显微镜(TIRM)由全内反射产生倏逝波,被观
图1模型问题的几何结构.散射体q(x)位于紧集D内,探针在r：=｛x=(x1,x2)：x2=b>0,-L 察物体由倏逝波照射,并在远场探测散射场.
光子扫描隧道显微镜(PSTM)可以看作是NSOM和TIRM的结合,被观察物体由倏逝波照射,并在近场探测散射场.
我们主要给出光子扫描隧道显微镜(PSTM)的数学模型.
假设被观察物体(散射体)置于均匀介质基座上,基座足够厚所以我们只考虑基座的一个面.在模型中,只考虑TM(transverse magnetic)极化情形,且不考虑探针的影响,如图1所示.
以基座的面作为空间的分界面,设下半空间的折射率为n0,上半空间在除散射体外的部分折射率为1[19,28].即等价地,可认为背景物质由双层介质组成,其折射率为
其中x=(x1,x2),n0>1.
记散射体为q(x),位于紧集D(?)碾R+2=｛x=(x1,x2)：x2>0}内,其中q(x)/4π为散射体的极化率(susceptibility)[19]以时谐平面波
作为入射波,从下半空间入射,其中
k0为真空中波数.总场u满足Helmholtz方程
记uref为参考场,即无散射体q时的总场,易知其满足Helmholtz方程：
由方程(2)以及界面的连续性条件,我们可以得到参考场的解析表达式
其中
和
分别为透射波和反射波[28-30],其中
易见当|α|>k0时,发生全内反射现象,透射波ut变为倏逝波,在xx方向指数衰减.定义散射场us=u-uref,由(1),(2)可得
同时散射场还要满足辐射条件[28,31]
其中∑R为圆心在原点半径为R的圆周,v是∑R上的单位外法向量.
2.谱截断矩量法及其应用
2.1谱截断矩量法
谱截断矩量法是一种易于实现的方法,形式上看,它是矩量法与谱截断的结合.
设Y=C[a,b],X为Hilbert空间,K：X→Y为有界线性算子,且是的.假设算子方程Kx=y存在唯一解,a≤s1＜s2＜…＜sn≤b为配置点,Xn(?X为n维子空间,则有配置方程
由Riesz定理,存在kj∈X,使得
取Xn=span{kj,j=1,…,n},则形成求解矩量解
的线性代数方程组
其中
之后对系数矩阵进行奇异值分解(SVD)：
其中U,V都是酉矩阵.且有
以及
其中λ1≥λ2≥…≥λn≥0对选定的正整数K(1≤K≤n),记α(δ):=∧2≥λK2,Σ=diag(λ1-1,λ2-1,…,λK-1,0,…,0),则(6)的谱截断解为
记
则xnα(δ)(δ)，δ为问题Kx=y的谱截断矩量解(moment solution with truncated SVD).对谱截断矩量法,我们有如下的误差估计：
定理1.设{s1(n),…,s(n)}(?)[a,b],n∈N为配置点列,子空间
a是问题(6)以β为右端项的解,假设存在z∈Cn,|z|≤E,使得a=A*z,则
设常数c>0,取^=(?)cδ/E有
为最优阶估计,其中
2.2谱截断矩量法在近场光学线性化问题中的应用
我们首先求解近场光学的线性化问题,即Born近似下的近场光学反问题.具体地,我们要从方程
中求解散射体q(x).注意到我们仅仅能在线段Γ={x:x2=b,-L＜x1＜L}上探测到散射场数据,所以我们在{x:x2=b,-∞＜x1＜∞}\Γ上将us补充定义为0.这里我们假设当L充分大时,补充定义的us与真实值的误差充分小,即有‖us-uexacts‖＜δ1,其中δ1依赖于L.由于散射数据与α有关,在反演散射体q(x)时我们需要多个角度入射.
于是,有
将基本解G(x,y)和ut的表达式代入方程(11),注意到us依赖于参数α,方程(11)可写为
在方程(12)的两端同时乘以e-iωx1并关于x1从-∞到∞积分,有
注意到
(13)可化为
注意到
(14)可化简为(15)
当|ω|＜k0以及|α|<‰＜k0时,β1和γ都是实数;当k0<＜|α|＜n0k时,γ是纯虚数;当|ω|>k0时β1也纯虚数.所以,积分方程(15)的求解既涉及到有限区间上的Fourier逆变换,又涉及到有限区间上的Laplace逆变换.而有限区间上的Fourier逆变换和Laplace逆变换都是不适定的,特别地,Laplace逆变换是指数度不适定的(参见2.1.2节与文献[53]).
我们使用谱截断矩量法,就在一定程度上避免了计算机最小精度的限制.
我们将二维的积分方程(15)分解为两个一维问题.这样计算效率会显著提高,计算中所占用内存也会大幅度下降.在计算其中第一个一维问题时,我们采用谱截断矩量法(moment method with truncated SVD),这样我们可以在一定程度上减少对计算机最小精度的依赖,并获得更好的分辨率.
首先,对每个ω-α=c(k),k=1,2,,积分方程(15)化为一维问题,对每个k,我们可以求解得到q(c(k),y2);
需要注意的是,在求解每个q(c(kl),y2),k=1,2,时,对不同的k,要选取不同的ω.如果我们有n个不同入射角,对应于αj,j=1,…,n,则对k=1,2,选取ωj=αj+c(k),j=1,…,n我们设所有的ωj包含于集合[-W,W]中.
第二步,通过有限区间上的Fourier逆变换,由q(c(k),y2)得到q(y1,y2)
经过以上两步计算,得到线性化问题(10)的数值解.以上两步的计算均是不适定的,特别第一步的不适定性更严重,一方面是Laplace逆变换的不适定,另一方面是探测数据的不完整(只在Γ上有限个点做探测),我们将数据缺失作为一种扰动来处理.
下面我们给出求解积分方程(15)的具体方法.首先,我们对每个k,k1,2,…,求解算子方程
其中(17)
记
由于k(ω,y2)∈L.([-W,W]×[0,b]),显然K:L2[0,b]→L2[-W,W]为紧算子设Xn(?)L2[0,b]为有限维子空间,其维数dim Xn=n设0≤ω1ωn≤b为配置点,其中ωj=αj+c(k),j=1,…,n于是配置方程为
其中
则qn∈Xn为问题(18)的配置解.
选择Xn的基底为{Φj,j=1,…,n},则qn(c(k),y2)可表为
则系数aj为线性代数方程组
的解,其中
对于不适定问题(19),对右端数据β很小的误差扰动都会使解不稳定,而我们只能在线段Γ探测数据,我们将数据的缺失也作为扰动来处理,即设有‖us-uexacts‖＜δ1为方便记,我们设带有噪声的探测数据为βδ,并设|β-βδ|<δ,其中δ包含了δ1的影响,这里的范数|·|为Cn空间的Euclidean范数.
我们使用谱截断矩量法求解问题(18).则Riesz定理,
其中-表示复共轭.如果kj,j=1,…,n线性无关,即系数矩阵Ai,j=(kj,ki)的行列式不为0,则问题(18)存在唯一的矩量解qn
设A的奇异值分解为
其中λ1≥…≥λn≥0,取正整数K(1≤K≤n),记α(δ):=∧2≥λk2, Σ=diag(λ1-1,λ2-1,…,λK-1,0,…,0),则
于是(18)的谱截断矩量解为
于是由定理1,可以得到关于近场光学线性化问题的误差估计：
定理2.设{ω1(n),ω2(n),…,ωn(n)}(?)[0,b],n∈N,为配置点列,
假设存在z∈Cn,|z|≤E,使得a=A*z.qn～δ是问题(18)的谱截断矩量解,则有如下的误差估计
如果对常数c>0,取∧=cδ/E,则有
其中
由谱截断法,可以选取合适的^,使得δ/Λ有界.除了选取∧=cδ/E外,我们还可以针对不同问题,更灵活地加以选择,比如可以选取∧=δα,0＜α＜1,这时有δ/∧+∧|z|～δ1-α+cδα
由上面的方法求解得到q(c(k),y2),k=1,2,后,可以利用Fourier逆变换公式,通过数值积分得到q(y1,y2)由Parseval公式,有‖qn～δ-q‖L22(D)=‖qn～0-q‖L22(D)/(2π)
数值实验表明该方法快捷,可靠,即使在噪声数据下也可以获较为满意的分辨率.
3.无穷维非线性问题的同伦方法及其应用
3.1带Tikhonov正则化的同伦
我们将无穷维空间的同伦方法与Tikhonov正则化结合到一起,构造一种解决无穷维不适定问题的同伦方法.
设X,Y为Hilbert空间,考虑非线性不适定问题
其中F:D(F)(?)X→Y为非线性不适定(紧)算子.构造极小化泛函
其中x0∈D(F),A(λ):X→X为依赖于参数λ∈[0,1]的算子.
注1.若取A(λ)=1-(1-α)λI,易见λ=1时,(25)恰为Tikhonov泛函.若同伦曲线存在,我们将得到问题(24)的Tikhonov正则化解.
假设F于D(F)内对任何x都Frechet可微,则对任何h∈X,有
当
时,J将取得极小.于是,我们得到极小化泛函(25)的法方程：
注2.对A(λ)=√1-(1-α)λI的情形,有
恰为求解(24)的Tikhonov正则化解的不动点同伦.
对于A(λ)=√1-(1-α)λI的情形,我们得到带Tikhonov正则化的不动点同伦：
定理3.设D(F)为有界连通开集,F(x)∈C2(D(F),Y),F弱闭且存在一个元素x∈D(F),使得F(x)=y.对x0∈D(F),设h(x,t)如(26)所示,其中A(λ)=1-(1-α)λI,如果对(x,λ)∈(?)D(F)×[0,1),有h(x,λ)≠0,则h(x,λ)=0的解集合中至少存在一条解曲线,连接x0到(24)的一个正则化解.
在同伦曲线的跟踪中,可以采用Seidel技巧,用最新得到的x(λ)代替同伦曲线(27)中的x0,以获得更好的计算效果.
3.2无导数同伦
带Tikhonov正则化的同伦中含有非线性映射F的一阶导数,在同伦曲线的跟踪过程中,可能还需要计算F的二阶导数,在理论和计算上都有不小的难度,本小节介绍一种不含F的导数的同伦方法,它主要源于这样一个事实：
我们对问题(24)直接构造“不动点”同伦：
其中0<ρ<1.显然对λ∈[0,1-ρ],(29)中h(x,λ)=0的计算均是适定的.与定理3的证明类似,我们可以得到
定理4.设D(F)为有界连通开集,F(x)∈C2(D(F),Y),对x0∈D(F),设h(x,t)如(29)所示,如果对(x,λ)∈(?)D(F)×[0,1),有h(x,λ)≠0,则h(x,λ)=0的解集合中至少存在一条解曲线,连接x0到x1
如果h(x,1-ρ)=0的解x1。与F(x)=y的一个解充分接近,我们可以x1。作为初始近似,采用某种正则化方法求解(24).
3.3同伦曲线的跟踪
对于同伦曲线,我们要使用合适的数值方法去跟踪求解.假设x可由λ参数化,同伦曲线为h(x(λ),λ)=0,其中h(x,0)=0的解x(0)已知,我们要对0=λ0＜λ1＜…＜λn=1,依次求解h(x,λj)=0,j=1,…,n的解x(j),j=1,…,n.常见的方法有Stefenssen方法,牛顿法,欧拉法等.
Stefenssen方法
Stefenssen方法是从j=1,…,n顺次求解h(x,-)=0的解x(j),所以总可以x(j-1)作为初始近似来求解.具体地,我们采用一种类似求解非线性代数方程组的Stefenssen方法.对l=0,1,…,令
每一步迭代中的初始近似x(j-1)与x(j)足够接近是收敛的必要条件,而对于不适定问题,要想使x(j-1)成为h(x,λj)=0的足够接近的初值,就必须要λj-λj-1足够小,但这样必然会带来较大的计算量,我们会采用一种自适应的方法来处理这个问题.
牛顿迭代法
与简单迭代法的思想类似,假设h(x,λj-1)=0的解x(j-1)已经求得,我们将其作为一个初始近似,采用具有局部收敛的牛顿法求解h(x,λj)=0有
特别地,对于同伦(27),我们有
于是,从x(j-1)计算到h(x,λj)=0的解x(j)的牛顿法为：
当‖Δl(j)‖充分小时设定x(j)-xl(j),并进入从x(j)到x(j+1)的计算过程.可见牛顿迭代法中每一步迭代需要解一个关于△l(j)的线性算子方程.
对于同伦(29),有
于是,从x(j-1)计算到h(x,λj)=0的解x(j)的牛顿法为：
‖Δl(j)‖充分小时设定x(j)=xl(j),并进入从x(j)到x(j+1)的计算过程.
从λ=1-δ到λ=1的跟踪过程中,需要求解的线性算子方程变为
是一个不适定线性算子方程,我们可以采用某种正则化策略予以求解,如Tikhonov正则化,谱截断法等.
欧拉法
欧拉法是首先对h(x(λ),λ)=0关于λ求导,得到一个形式上的常微分方程组,
其中x0已知,再对其采用欧拉折线法进行数值求解.特别地,对同伦(27),有
对于给定的点x,上述方程是关于dx的线性算子方程.对λ选取合适的剖分步长,即对0=λ0<λ1<…<λn=1,有
对于同伦(29),有
对于给定的点x,上述方程是关于dx的线性算子方程.对λ选取合适的剖分步长,即对0=λ0<λ1＜…＜λn-1=1-δ,有
从λ=1-δ到λ=1的跟踪过程中,需要求解的线性算子方程变为
是一个不适定线性算子方程,我们采用某种正则化策略予以求解,如Tikhonov正则化,谱截断法等.
在从λ=λj-1到λ=λj的迭代步中,每一步迭代中的初始近似x(j-1)与x(j)都要足够接近才能取得较好的计算效果,这就要λj-λj1足够小,但这样必然会带来较大的计算量,这里给出一种自适应的方法.
第1步：选定阈值ε>0,N为正整数.对λ∈[0,1].给出初始粗网格0=λ0＜λ1＜…＜λn=1,并令j=0;
第2步：如果j=n,终止计算;否则进入第3步;
第3步：计算h(x(λj),λj),如果h(x(λj),λj)＜∈,则令j=j-＞next,返回第2步;否则将[λj-1,λj]再次剖分,将-1/2加入计算网格,形成新的较细的网格
令j=j-＞former,返回第2步;若细分次数大于给定阈值N,终止运算,重新选择初值x0与正则化参数α.
3.4同伦方法在近场光学问题中的应用
本节将同伦方法应用于近场光学问题.我们已经给出近场扫描隧道光学显微镜(PSTM)的数学模型：
以及辐射条件
其中∑R为圆心在原点半径为R的圆周,v是∑R上的单位外法向量.
如果记(32)(33)的解us为S(q),则显然S(q)：L2(D)→L2(R2)是关于q的非线性映射.我们在前一部分针对线性化问题的求解进行了讨论,主要克服其严重不适定性.本节我们讨论处理其非线性的同伦方法.同伦的构造和同伦曲线的跟踪方法已在前节说明,本节主要给出在计算过程中所需要的各阶导数的计算过程.
首先,我们推导S’(q).由于us=S(q)满足(32)(33),故对h∈D(?)Ω(?)R+2, us(q+h)=S(q+h)满足
(34)-(32)并省去高阶项,有
记v=S’(q)h,则v满足
以及辐射条件.上式中右端将n2省略是因为于R+2中n=1(注意h∈D(?)R+2).
由此,可以继续推导二阶导数S”(q).首先假设S”(q)存在.设V1=S’(q+h1)h,其中h1,h∈D(?)Ω(?)R+2,则v1满足
以及辐射条件.(37)-(36)得
利用关系式
舍去高阶项,并令w=S”(q)h1h,有
其中v=S’(q)h1满足
其中w和v均满足辐射条件.
在同伦方法中,需要求解S'(q)*(S(q)-y),其中y是给定测量值.仍记v=S’(q)h1,则其满足(41).在(41)两端乘以(?),并在R2上积分,利用Green公式及辐射条件,有
由关系式
和S'(q)h1=v可知,对给定S(q)-y,求(?)满足
和辐射条件,则
这样,对于带Tikhonov正则化的同伦,即
我们得到h(q(λ),λ)的计算方法：由于(43)(44)式已经给出S'(q)*(S(g)-yδ)的表达式,我们需要按(43)和辐射条件求解微分方程正问题,得到(?),再将(?)代入(44)得到S'(q)*(S(g)-yδ),之后直接运算即得到h(q(λ),λ).
接下来,我们按照自适应跟踪技巧,使用Stefenssen方法进行曲线的跟踪.
第1步：选定阈值∈,∈',N,N',初始网格为0=λ0<λ1<…<λ。=1,并令j=0;
第2步：如果j=n,终止计算;否则进入第3步;
第3步：以qj→former为初值,使用Stefenssen方法计算q(λj).具体地,对l=0,1,…,令
若‖ql+1(j)-ql(j)‖＜∈'且迭代次数l<=N’,则令q(λj)=ql+1(j),否则认为h(q(λj),λj)＞
如果h(q(λj),λj)＜∈,则令j=j→＞next,返回第2步;否则将[λj-1,λj]再次剖分,将λj1/2加入计算网格,形成新的较细的网格
令j=j→＞former,返回第2步;若细分次数大于给定阈值N,终止运算,重新选择初值q0与正则化参数α.
这样,我们得到近场光学问题的同伦算法.
In this paper, the numerical method for solving nonlinear ill-posed problems isproposed, analyzed and applied to the near-field optics. Chapter1is an introduction toscattering problems and inverse scattering problems, the model of near-field optics, theconcept of ill-posed problems and regularization method. Chapter2introduces the mo-ment method with truncated SVD and its applications to severely ill-posed problems,especially in near-field optics. From the numerical experiments, we can conclude thatthemomentmethodwithtruncatedSVDcandealwithseverelyill-posedproblemsmoreeffectively, comparedtotheconventionalregularizationmethodssuchasTikhonovreg-ularization or moments method. Chapter3introduces infinite dimensional homotopymethod to solving nonlinear ill-posed problems. There are two main aspects in ho-motopy method: one is the existence of the homotopy path, the other is the numericalmethod for tracking the homotopy path. In this chapter, two homotopy methods forsolving ill-posed problems are discussed, namely homotopy with Tikhonov regular-ization and homotopy without derivative, and several tracking methods are given forhomotopy path. In the tracking process, adaptive tracking skills are given in order toimprove computing efficiency.
1. Scattering Model of near-field optics
Near-field optics is a new subject since the1980s, breaking the traditional limit ofoptical resolution. There are three important modalities in near-field optics: near-field scanning optical microscopy(NSOM)[23], total internal reflection microscopy(TIRM)[24,25], and photon scanning tunneling microscopy (PSTM)[26,27].
The two principal genres of NSOM are illumination mode and collection mode. In illumination mode NSOM, a probe is scanned over the scatterer in the near zone and the scattering field is measured in the far zone of scatterer. In collection mode NSOM, the scatterer is illuminated by a source located in the far zone of the scatterer and the total field is collected in the near zone of the scatterer through the small aperture in the probe as the probe is scanned over the scatterer. In TIRM, the scatterer is illuminated by an evanescent wave that is generated by total internal reflection, and the scattering field is measured in the far zone of the scatterer. In the PSTM, the scatterer is illuminated by an evanescent wave and the scattering field is measured in the near zone.
In all of the three modalities, reconstruction of the scatterer from the measured data is desirable.
We describe the scattering model mathematically as follows. Assume that the ma-terials are nonmagnetic in the problem considered in this paper. Throughout, we only consider the problem which can be reduced to a two-dimensional model in the case of transverse magnetic(TM) polarization, as is shown in Figure1.
The index of refraction in the lower half-space has a constant value n0. The index of refraction in the upper half-space varies within the domain of the sample but other-wise has a value of unity[19,28]. That is, the background media consists of two layers with the index of refraction where x=(x1,x2), n0＞1.
We denote the scatterer by q(x) within a compact support in D(?)R2+:={x
{x1,x2):x2＞0}, with q(x)/4π the susceptibility of the scatterer [19]. The incident wave comes from the lower half space, which is a time-harmonic plane wave where
Figure1The geometry of the model problem. The scatterer q(x) is in a compact support in D, the measurements are taken on Γ, and is a compact set containing D and Γ. and θ∈(-π/2,π/2), k0is the free space wavenumber. The geometry of the model problem is shown in Figure1. The measurements are taken on Γ:={x=(x1,x2): x2=b＞0,-L＜xi＜L}, and Ω(?)R2is a compact set containing D and Γ.
It is easy to check that the total field u satisfies the Helmholtz equation
Let uref be the reference field, which satisfies the Helmholtz equation without the scat-terer:
From Eq.(2) and the continuity condition on the interface we can analytically obtain
that
where
and are the transmitted and reflected waves respectively[28-30], and
It can be readily seen that the transmitted wave ut becomes evanescent wave when|α|＞k0. That is, ut decays exponentially in the x2direction.
Define the scattering field us=u-uref. Then from (1),(2), we have
And we require the radiation condition [28,31]
where ΣR is the sphere centered in the origin with radius R, and v is the unit outward normal to ΣR.
2. Moment method with truncated SVD and applications
2.1Moment method with truncated SVD
Moment method with truncated SVD is easy to implement. Formally, it is a com-bination of moment method and truncated SVD.
Let Y=C[a, b] and X be Hilbert spaces, the operator K:X→Y be bounded and one-to-one. Let a≤s1＜s2＜…＜sn≤b be the collocation points, and Xn C X be finite-dimensional subspaces with dimXn=n, then the collocation equation is
By Riesz Theorem, there exists kj E X, such that
Kz(sj)=(z,kj),(?)z∈X, j=1,…,n.
Let Xn=span{kj,j=1,…,n}, then the moment solution is given by
where a solves the linear system
with
Then the singular value decomposition(SVD)of A is made by
where U,V are unitary with
and
where λ1≥λ2≥…≥λn≥0.For given K(1≤K≤n),denote α(δ):=∧2≥λK2, Σ=diag(λ1-1,λ2-1,…,λK-1,0,…,0),then the solution of(6)by SVD is
Denote
then xnα(δ),δ is the moment solution with truncated SVD of equation Kx=y.
Now an error estimate of the Moment method with truncated SVD is given as follows:
Theorem1.Let{sn(n),…,sn(n))(?)[a,b],n∈N be a sequence of collocation points, and
Suppose that a solves(6)with β the right-hand side,and assume there exists z∈Cn|z|≤E,such that a=A*z,then
For constant c＞0,choose∧=cδ/E,then
which is optimal,where
2.2Application to linearized near-field optics problem
The linearized near-field optics problem,i.e. inverse problem concerning Born approximation,is to solve the scatterer q(x)from
Note that the only known data lies on the segment Γ={x:x2=b,-L＜x1＜L},so we extend us on Γ to the whole line{x:x2=b,-∞＜x1＜∞}by defining us=0on the remaining interval of Γ.We assume that the error between the expanded data us and the "exact" data is small when L is sufficiently large.Namely,we assume that‖us-uexacts‖＜δ1,where δ1depends on L.The scattering data us depends on α,and we need several different incident waves,i.e.,we need[α1,…,αn]∈[-n0k0,n0k0]to reconstruct the scatterer q(x).
Thus we have
Substituting the fundamental solution G(x,y)and ul into(11),and realizing the depen-dence of us on α,we rewrite(11)as
Multiplying both sides of(12)by e-iωx1and integrating from-∞to∞with respect to x1,
and noting that
(13) can be written as
Realizing that
(14) can be reduced to
When|ω|＜k0and|α|＜k0, both β1and γ are real; and for k0＜|α|＜nok,γ is pure imaginary, and so is β1with|ω|＞k0. It is obvious that, the inversion of (15) involves both inverse Fourier transform and inverse Laplace transform. In fact, both inverse Fourier transform and inverse Laplace transform defined on finite interval are ill posed, especially the inverse Laplace transform is exponentially ill posed(see section2.1.2and [53]).
Here, we split the two-dimensional problem (15) into two one-dimensional prob-lems. So the computation runs much faster and uses less memory. Moreover, we apply the moment method with truncated SVD to get a better resolution where the limitation of numerical precisions is circumvented in a certain extent.
Step1:For ω-α=c(k),k=1,2,…,obtain q(c(k),y2) from a one-dimensional integral equation for every k;
So for every calculation of q(c(k),y2), k=1,2,…, we need special choices of ω That is, if we have n incident angles, then we set ωj=αj+c(k),j=1,…,n, for k=1,2,…. We suppose all the ωj's lie in [-W, W].
Step2:Obtain q(y1,y2) from q(c(k),y2) by inverse Fourier transform on finite intervals.
After the two steps, an approximate solution of the linearized inverse problem (10) is obtained. Obviously, the computation of step1is quite difficult for its severe ill-posedness,especially on a finite data set.As is discussed above,we treat it as a kind of perturbation of the measured data.
The first step is to solve the operator equation
for every specified k,k=1,2,…,where Denote then K:L2[0,b]→L2[-W,W]is a compact operator since k(ω,y2)∈L2([-W,W]×
[0,b]).
Let Xn(?)L2[0,b]be finite-dimensional subspaces with dim Xn=n,and0≤ω1＜…＜ωn≤b be the collocation points,which are chosen as ωj=αj+c(k),j=1,…,n.Then the collocation equation is where
then qn∈Xn is a collocation solution of(18).
Choosing a basis{Φj,j=1,…,n)of Xn,then qn(c(k),y2)can be represented as
where the coefficients aj's are the solution of the system
with
For ill-posed problems such as(19),even small perturbation of data β can change the solution a lot.Recall that we can only obtain the scattering data on Γ,and we have assumed that‖us-uexacts‖＜δ1.Now we denote the perturbed data,or the measured data,by βδ,and let|β-βδ|＜δ,where δtakes in the influence of the error bound δ1, and|·|denotes Euclidean norm in Cn
We will find the moment solution with truncated SVD to(18).By Riesz Theorem,
where is conjugate.If kj,j=1,…,n are linearly independent,i.e.,the determinant of the collocation matrix∧Ai,j=(kj,ki)is not equal to0,by choosing Φj=kj,there exists one and only one moment solution qn of(18).
Using SVD,we can write
where λ1≥…≥λn≥0,let α(δ):=∧2≥λK2for some K(1≤K≤n), Σ=diag(λ1-1,λ2-1,…,λK-1,0,…,0),then
So the moment solution with truncated SVD to(18)is given by
By using Theorem1,we have the following estimate for linearized near-field optics problem:
Theorem2.Let{ω1(n),ω2(n),…,ωn(n)}(?)[0,b],n∈N,be a sequence of collocation points and
And assume that there exist z∈Cn,|z|≤E,such that a=A*z.Let qn～δ be the moment solution with truncated SVD to(18),then the following error estimate holds:
For some constant c＞0,choose∧=cδ/E,then
where
By our truncated SVD method,an appropriate∧could be chosen so that δ/∧is bounded.In practice,we could also choose∧=δα,0＜α＜1,then δ/∧+∧｜z δ1-α+cδα is independent of n. ically by inverse Fourier transform. By Parseval formula,it is obvious that‖qn～δ q‖L22(D)＝‖qn～δ-q‖L22(D)/(2π).
Numerical experiments show that this method computes fast,stable,and can get satisfying resolution even under the perturbed data.
3.Homotopy method for infinite dimensional nonlinear ill-posed problems and applications
3.1Homotopy with Tikhonov regularization
We combine homotopy method in infinite dimensional spaces with Tikhonov reg-ularization to solve infinite dimensional ill-posed problems.
Let X,Y be Hilbert spaces,consider nonlinear ill-posed problem
where F:D(F)(?)X→Y is nonlinear ill-posed(compact)operator.Construct minimization funcitional
where x0∈D(F),A(λ):X→X is an operator dependent to λ∈[0,1].
Remark1.If we set A(λ)=1-(1-α)λI,the when λ=1,(25)is nothing but Tikhonov functional.If the homotopy path exists,the Tikhonov regularization solution to(24)is obtained.
Assume that for every x∈D(F),F is Frechet differentiable,then for every h∈X,
J is minimized when Re(AF'(x)*(F(x)-yδ)+A*(λ)A(λ)(x-x0))=0.
Thus,the normal equation of(25)is given by:
Remark2.For the case when A(λ)=1-(1-α)λI, which is just the fixed-point homotopy with Tikhonov regularization of(24).
For the case when A(λ)=1-(1-α)λI,,the fixed-point homotopy with Tikhonov regularization of(24)is obtained:
Theorem3.LET D(F)be a bounded connected open set,F(x)∈C2(D(F),Y),F weakly closed, and there exists an element x∈D(F),such that F(x)=y.For x0∈D(F),denote h(x,t)as in(26),where A(λ)=1-(1-α)λI.If h(x,λ)≠0when(x,λ)∈(?)D(F)×[0,1),then there exists at least one path in the solution set of h(x,λ)=0,connecting x0to one regularization solution to(24).
In the tracking of the homotopy path,the Seidel skill can be adopted,i.e.,using the latest x(λ)instead of x0in(27)in order to obtain better solution.
3.2Homotopy without derivative
Homotopy with Tikhonov regularization contains the derivative of nonlinear map-ping F,moreover, in the tracking process,the second order derivative of F may have to be calculated,which makes the computing rather difficult. So homotopy without derivative is desirable.
Construct formally the fixed-point homotopy of(24):
wheye0＜ρ＜1.Obviously,for λ∈[0,1-ρ],the problem of solving h(x,λ)=0in
(29)is well posed.
We can obtain the following result similar to Theorem3:
Theorem4.Let D(F)be abounded connected open set,F(x)∈C2(D(F),Y).For x0∈D(F),denote h(x,t)as in(29).If h(x,λ)≠0when(x,λ)∈(?)D(F)×[0,1), then there exists at least one path in the solution set of h(x,λ)＝0,connecting x0to x1-ρ.
If the solution x1-ρ of h(x,1-ρ)=0is sufficiently close to the solution of F(x)=y,we can take x1-ρ as the initial approximation,and using some regularization method to solve(24).
3.3Tracking the homotopy path
Numerically,we need appropriate method to track the homotopy path.Assume that x can be parameterized by λ,i.e.,the homotopy path can be written as h(x(λ),λ)=0, with x(0)the solution of h(x,0)=0,which is given.For0=λ0＜λ1＜…＜λn=1, we need to solve h(x,λj)=0,j=1,…,n one by one,whose solutions are denoted by x(j),j=1,…,n.We introduce some methods such as Stefenssen method,Newton method,and Euler method.
Stefenssen method
Stefenssen method means solving h(x,λj)=0sequentially by j=1,…,n, whose solutions are denoted by x(j),so x(j-1)could be an initial approximation.Namely, for l=0,1,…,let
In each iteration,the initial approximation x(j-1)should be close to x(j)enough, but for ill-posed problems,we have to make λj-λj1small enough to make it the case,which brings over a large amount of computation,so we will adopt an adaptive approach to deal with this problem.
Adaptive tracking skills
As is discussed,we have to make λj-λj1small enough to make the initial ap-proximation x(j-1)close to x(j),but this brings over a large amount of computation.We introduce the adaptive tracking skills as follows.
Step1:Select threshold∈＞0,N a positive integer.For λ∈[0,1],give initial coarse grid0=λ0＜λ1＜…＜λn=1,and let j=0;
Step2:When j=n,terminate the calculation;otherwise go to Step3;
Step3:Calculate h(x(λj),λj),if h(x(λj),λj)＜∈,then let j=j-＞next,return to step2;otherwise split[λj-1,λj],and put λj-1/2into the new computing grid,i.e.,a new finer grid is formed as
then let j=j-＞former,return to step2;If the subdivision number is greater than the given threshold N,terminate the calculation,and re-select the initial x0and regularization parameter α.
Newton Method
Similar to Stefenssen method,we assume that the solution to h(x,λj-1)=0,i.e., x(j-1)has been obtained,and use it as an initial approximation to solve h(x,λj)=0. Here we use the local convergent Newton's method.Let
In particular, for the homotopy(27),
Thus,the Newton's method for solving h(x,λj)=0from the initial approximation
x(j-1)can be written as:
When‖△l(j)‖is sufficiently small we set x(j)=xl(j),and compute x(j+1)from the initial approximation x(j).It is readily that every iteration in Newton's method needs to solve a linear operator equation with respect to△l(j)
For homotopy(29),
Thus,the Newton's method for solving h(x,λj)=0from the initial approximation x(j-1)can be written as:
When‖△l(j)‖is sufficiently small we set x(j)=xl(j), and compute x(j+1)from the initial approxim ation x(j).
From λ=1-δ to λ=1,the linear operator equation becomes
which are ill-posed linear operator equations,so we can solve them by some kind of regularization method,such as Tikhonov regularization,SVD method.
Euler method
By Euler method,we derivative h(x(λ),λ)=0with respect to λ,and obtain a
where x0is given. Then solve it by Euler method numerically.
In particular, for the homotopy (27), we have
For a fixed point x, the above equation is a linear operator equation with respect to dx. For0=λ0＜λ1＜…＜λn=1, we have
For the homotopy (29), we have
For a fixed point x, the above equation is linear operator equation with respect to dx. For0=λ0＜λ1＜…＜λn=1, we have
From λ=1-δ to λ=1, the linear operator equation becomes
which is an ill-posed linear operator equation, so we can solve it by some kind of regu-larization method, such as Tikhonov regularization, SVD method.
3.4Application to near-field optics
In this section we will apply the homotopy method to near-field optics. The model of photon scanning tunneling microscopy (PSTM) has been given by
with radiation condition
where ΣR is the sphere centered in the origin with radius R,and v is the unit outward normal to∑R.
Denote the solution to(32)(33),i.e.us by S(q),then it is readily seen that S(q) L2(D)→L2(R2)is nonlinear in q.
In the previous part we solved the problem for the linearized case,mainly to over-come its severe ill-posedness.In this section we discuss the homotopy method to deal with the nonlinearity.Homotopy method and homotopy path have been given in the previous section,so in this section we calculate the derivatives in need.
First,we derive S'(q).
since us=S(q)satisfies(32)(33),so for h∈D(?)Ω(?)R+2,us(q+h)=S(q+h) satisfies
(34)-(32)and omit higher order terms,we get
△(S'(q)h)+n22k02(1+q)(S'(q)h)+n2k02hS(q)=+k02hut,in R+2.(35) Denote v=S'(q)h,then v satisfies
△v+n2k02(1+q)v=-k02hut-n2k2hS(q)=-k02h(ut+us(q)),in R2+,(36)
and radiation condition.
Thus,we can continue to derive the second derivative.Assume that S"(q)exists. Let v1=S'(q+h1)h,where h1,h∈D(?)Ω(?)R+2,then v1satisfies
and radiation condition.(37)-(36)and we get
Using the relationship
and omitting higher order terms, and denoting ω=S"(q)h1h, we have
where v=S'(q)h1satisfies
where both w and v satisfy radiation condition.
In homotopy method, we need to solve S'(q)*(S(q)-y), where y is measured data. Again, denote v=S'(q)h1, then it satisfies (41). Multiply0on both sides of (41), and integrate on R2, by Green's formula and radiation condition,
By the relationship
and S'(q)h1=v, we can get that, for given S(q)-y, solve Φ that satisfies
and radiation condition, then
Thus, for homotopy with Tikhonov regularization, i.e.,
the numerical method for computing h(q(λ),λ) has been obtained:since the formula of S'(q)*(S(q)-yδ) has been given in (43) and (44), we solve0from PDE (43) and radiation condition, then substitute0into (44) and get S'(q)*(S(q)-yδ), and then h(q(λ),λ).
Next, by the adaptive tracking skills, we use Stefenssen method to track the homo-topy path.
Step1:Select threshold∈,∈',and N,N' positive integers,and the initial coarse grid0=λ0＞λ1＜…＜λn=1,and let j=0;
Step2:When j=n,terminate the calculation;otherwise go to Step3;
Step3:Calculate q(λj)by Stefenssen method,using qj-＞former as initial guess, namely,for l=0,1,…,let
if‖ql+1(j)-ql(j)‖＜∈' and the iteration steps l＜=N',then let q(λj)=ql-1(j),otherwise let h(q(λj),λj)＞∈.
If h(q(λj),λj)＜∈,then let j=j-＞nexl,and return to step2;otherwise split [λj-1,λj],and put λj-1/2into the new computing grid,i.e.,a new finer grid is formed as
then let j=j-＞former,and return to step2;If the subdivision number is greater than the given threshold N,terminate the calculation,and re-select the initial q0and regularization parameter α.
In this way,the homotopy algorithm for near-field optics is obtained.
1.近场光学的数学模型
近场光学是上世纪80年代以来出现的一个新兴学科,突破了传统光学的分辨率极限.近场光学显微镜有三种主要的模型：近场扫描光学显微镜(NSOM)[23],全内反射显微镜(TIRM)[24,25],以及光子扫描隧道显微镜(PSTM)[26,27].
NSOM有两种主要类型：照射模式(illumination mode) NSOM和接受模式(collection mode) NSOM照射模式NSOM中具有细小的光学探针,其尖端的孔径远小于光的波长,探针作为光源在被观察物体(散射体)的近场区域扫描,散射场将被探测并记录.接受模式NSOM的被观察物体由在远场的光源照射,其探针在被观察物体的近场区域探测总场(或散射场).
我们知道,当光线以大于临界角度射入一个棱镜时,光线会被发生全内反射现象,这时将产生倏逝波.全内反射显微镜(TIRM)由全内反射产生倏逝波,被观
图1模型问题的几何结构.散射体q(x)位于紧集D内,探针在r：=｛x=(x1,x2)：x2=b>0,-L
光子扫描隧道显微镜(PSTM)可以看作是NSOM和TIRM的结合,被观察物体由倏逝波照射,并在近场探测散射场.
我们主要给出光子扫描隧道显微镜(PSTM)的数学模型.
假设被观察物体(散射体)置于均匀介质基座上,基座足够厚所以我们只考虑基座的一个面.在模型中,只考虑TM(transverse magnetic)极化情形,且不考虑探针的影响,如图1所示.
以基座的面作为空间的分界面,设下半空间的折射率为n0,上半空间在除散射体外的部分折射率为1[19,28].即等价地,可认为背景物质由双层介质组成,其折射率为
其中x=(x1,x2),n0>1.
记散射体为q(x),位于紧集D(?)碾R+2=｛x=(x1,x2)：x2>0}内,其中q(x)/4π为散射体的极化率(susceptibility)[19]以时谐平面波
作为入射波,从下半空间入射,其中
k0为真空中波数.总场u满足Helmholtz方程
记uref为参考场,即无散射体q时的总场,易知其满足Helmholtz方程：
由方程(2)以及界面的连续性条件,我们可以得到参考场的解析表达式
其中
和
分别为透射波和反射波[28-30],其中
易见当|α|>k0时,发生全内反射现象,透射波ut变为倏逝波,在xx方向指数衰减.定义散射场us=u-uref,由(1),(2)可得
同时散射场还要满足辐射条件[28,31]
其中∑R为圆心在原点半径为R的圆周,v是∑R上的单位外法向量.
2.谱截断矩量法及其应用
2.1谱截断矩量法
谱截断矩量法是一种易于实现的方法,形式上看,它是矩量法与谱截断的结合.
设Y=C[a,b],X为Hilbert空间,K：X→Y为有界线性算子,且是的.假设算子方程Kx=y存在唯一解,a≤s1＜s2＜…＜sn≤b为配置点,Xn(?X为n维子空间,则有配置方程
由Riesz定理,存在kj∈X,使得
取Xn=span{kj,j=1,…,n},则形成求解矩量解
的线性代数方程组
其中
之后对系数矩阵进行奇异值分解(SVD)：
其中U,V都是酉矩阵.且有
以及
其中λ1≥λ2≥…≥λn≥0对选定的正整数K(1≤K≤n),记α(δ):=∧2≥λK2,Σ=diag(λ1-1,λ2-1,…,λK-1,0,…,0),则(6)的谱截断解为
记
则xnα(δ)(δ)，δ为问题Kx=y的谱截断矩量解(moment solution with truncated SVD).对谱截断矩量法,我们有如下的误差估计：
定理1.设{s1(n),…,s(n)}(?)[a,b],n∈N为配置点列,子空间
a是问题(6)以β为右端项的解,假设存在z∈Cn,|z|≤E,使得a=A*z,则
设常数c>0,取^=(?)cδ/E有
为最优阶估计,其中
2.2谱截断矩量法在近场光学线性化问题中的应用
我们首先求解近场光学的线性化问题,即Born近似下的近场光学反问题.具体地,我们要从方程
中求解散射体q(x).注意到我们仅仅能在线段Γ={x:x2=b,-L＜x1＜L}上探测到散射场数据,所以我们在{x:x2=b,-∞＜x1＜∞}\Γ上将us补充定义为0.这里我们假设当L充分大时,补充定义的us与真实值的误差充分小,即有‖us-uexacts‖＜δ1,其中δ1依赖于L.由于散射数据与α有关,在反演散射体q(x)时我们需要多个角度入射.
于是,有
将基本解G(x,y)和ut的表达式代入方程(11),注意到us依赖于参数α,方程(11)可写为
在方程(12)的两端同时乘以e-iωx1并关于x1从-∞到∞积分,有
注意到
(13)可化为
注意到
(14)可化简为(15)
当|ω|＜k0以及|α|<‰＜k0时,β1和γ都是实数;当k0<＜|α|＜n0k时,γ是纯虚数;当|ω|>k0时β1也纯虚数.所以,积分方程(15)的求解既涉及到有限区间上的Fourier逆变换,又涉及到有限区间上的Laplace逆变换.而有限区间上的Fourier逆变换和Laplace逆变换都是不适定的,特别地,Laplace逆变换是指数度不适定的(参见2.1.2节与文献[53]).
我们使用谱截断矩量法,就在一定程度上避免了计算机最小精度的限制.
我们将二维的积分方程(15)分解为两个一维问题.这样计算效率会显著提高,计算中所占用内存也会大幅度下降.在计算其中第一个一维问题时,我们采用谱截断矩量法(moment method with truncated SVD),这样我们可以在一定程度上减少对计算机最小精度的依赖,并获得更好的分辨率.
首先,对每个ω-α=c(k),k=1,2,,积分方程(15)化为一维问题,对每个k,我们可以求解得到q(c(k),y2);
需要注意的是,在求解每个q(c(kl),y2),k=1,2,时,对不同的k,要选取不同的ω.如果我们有n个不同入射角,对应于αj,j=1,…,n,则对k=1,2,选取ωj=αj+c(k),j=1,…,n我们设所有的ωj包含于集合[-W,W]中.
第二步,通过有限区间上的Fourier逆变换,由q(c(k),y2)得到q(y1,y2)
经过以上两步计算,得到线性化问题(10)的数值解.以上两步的计算均是不适定的,特别第一步的不适定性更严重,一方面是Laplace逆变换的不适定,另一方面是探测数据的不完整(只在Γ上有限个点做探测),我们将数据缺失作为一种扰动来处理.
下面我们给出求解积分方程(15)的具体方法.首先,我们对每个k,k1,2,…,求解算子方程
其中(17)
记
由于k(ω,y2)∈L.([-W,W]×[0,b]),显然K:L2[0,b]→L2[-W,W]为紧算子设Xn(?)L2[0,b]为有限维子空间,其维数dim Xn=n设0≤ω1ωn≤b为配置点,其中ωj=αj+c(k),j=1,…,n于是配置方程为
其中
则qn∈Xn为问题(18)的配置解.
选择Xn的基底为{Φj,j=1,…,n},则qn(c(k),y2)可表为
则系数aj为线性代数方程组
的解,其中
对于不适定问题(19),对右端数据β很小的误差扰动都会使解不稳定,而我们只能在线段Γ探测数据,我们将数据的缺失也作为扰动来处理,即设有‖us-uexacts‖＜δ1为方便记,我们设带有噪声的探测数据为βδ,并设|β-βδ|<δ,其中δ包含了δ1的影响,这里的范数|·|为Cn空间的Euclidean范数.
我们使用谱截断矩量法求解问题(18).则Riesz定理,
其中-表示复共轭.如果kj,j=1,…,n线性无关,即系数矩阵Ai,j=(kj,ki)的行列式不为0,则问题(18)存在唯一的矩量解qn
设A的奇异值分解为
其中λ1≥…≥λn≥0,取正整数K(1≤K≤n),记α(δ):=∧2≥λk2, Σ=diag(λ1-1,λ2-1,…,λK-1,0,…,0),则
于是(18)的谱截断矩量解为
于是由定理1,可以得到关于近场光学线性化问题的误差估计：
定理2.设{ω1(n),ω2(n),…,ωn(n)}(?)[0,b],n∈N,为配置点列,
假设存在z∈Cn,|z|≤E,使得a=A*z.qn～δ是问题(18)的谱截断矩量解,则有如下的误差估计
如果对常数c>0,取∧=cδ/E,则有
其中
由谱截断法,可以选取合适的^,使得δ/Λ有界.除了选取∧=cδ/E外,我们还可以针对不同问题,更灵活地加以选择,比如可以选取∧=δα,0＜α＜1,这时有δ/∧+∧|z|～δ1-α+cδα
由上面的方法求解得到q(c(k),y2),k=1,2,后,可以利用Fourier逆变换公式,通过数值积分得到q(y1,y2)由Parseval公式,有‖qn～δ-q‖L22(D)=‖qn～0-q‖L22(D)/(2π)
数值实验表明该方法快捷,可靠,即使在噪声数据下也可以获较为满意的分辨率.
3.无穷维非线性问题的同伦方法及其应用
3.1带Tikhonov正则化的同伦
我们将无穷维空间的同伦方法与Tikhonov正则化结合到一起,构造一种解决无穷维不适定问题的同伦方法.
设X,Y为Hilbert空间,考虑非线性不适定问题
其中F:D(F)(?)X→Y为非线性不适定(紧)算子.构造极小化泛函
其中x0∈D(F),A(λ):X→X为依赖于参数λ∈[0,1]的算子.
注1.若取A(λ)=1-(1-α)λI,易见λ=1时,(25)恰为Tikhonov泛函.若同伦曲线存在,我们将得到问题(24)的Tikhonov正则化解.
假设F于D(F)内对任何x都Frechet可微,则对任何h∈X,有
当
时,J将取得极小.于是,我们得到极小化泛函(25)的法方程：
注2.对A(λ)=√1-(1-α)λI的情形,有
恰为求解(24)的Tikhonov正则化解的不动点同伦.
对于A(λ)=√1-(1-α)λI的情形,我们得到带Tikhonov正则化的不动点同伦：
定理3.设D(F)为有界连通开集,F(x)∈C2(D(F),Y),F弱闭且存在一个元素x∈D(F),使得F(x)=y.对x0∈D(F),设h(x,t)如(26)所示,其中A(λ)=1-(1-α)λI,如果对(x,λ)∈(?)D(F)×[0,1),有h(x,λ)≠0,则h(x,λ)=0的解集合中至少存在一条解曲线,连接x0到(24)的一个正则化解.
在同伦曲线的跟踪中,可以采用Seidel技巧,用最新得到的x(λ)代替同伦曲线(27)中的x0,以获得更好的计算效果.
3.2无导数同伦
带Tikhonov正则化的同伦中含有非线性映射F的一阶导数,在同伦曲线的跟踪过程中,可能还需要计算F的二阶导数,在理论和计算上都有不小的难度,本小节介绍一种不含F的导数的同伦方法,它主要源于这样一个事实：
我们对问题(24)直接构造“不动点”同伦：
其中0<ρ<1.显然对λ∈[0,1-ρ],(29)中h(x,λ)=0的计算均是适定的.与定理3的证明类似,我们可以得到
定理4.设D(F)为有界连通开集,F(x)∈C2(D(F),Y),对x0∈D(F),设h(x,t)如(29)所示,如果对(x,λ)∈(?)D(F)×[0,1),有h(x,λ)≠0,则h(x,λ)=0的解集合中至少存在一条解曲线,连接x0到x1
如果h(x,1-ρ)=0的解x1。与F(x)=y的一个解充分接近,我们可以x1。作为初始近似,采用某种正则化方法求解(24).
3.3同伦曲线的跟踪
对于同伦曲线,我们要使用合适的数值方法去跟踪求解.假设x可由λ参数化,同伦曲线为h(x(λ),λ)=0,其中h(x,0)=0的解x(0)已知,我们要对0=λ0＜λ1＜…＜λn=1,依次求解h(x,λj)=0,j=1,…,n的解x(j),j=1,…,n.常见的方法有Stefenssen方法,牛顿法,欧拉法等.
Stefenssen方法
Stefenssen方法是从j=1,…,n顺次求解h(x,-)=0的解x(j),所以总可以x(j-1)作为初始近似来求解.具体地,我们采用一种类似求解非线性代数方程组的Stefenssen方法.对l=0,1,…,令
每一步迭代中的初始近似x(j-1)与x(j)足够接近是收敛的必要条件,而对于不适定问题,要想使x(j-1)成为h(x,λj)=0的足够接近的初值,就必须要λj-λj-1足够小,但这样必然会带来较大的计算量,我们会采用一种自适应的方法来处理这个问题.
牛顿迭代法
与简单迭代法的思想类似,假设h(x,λj-1)=0的解x(j-1)已经求得,我们将其作为一个初始近似,采用具有局部收敛的牛顿法求解h(x,λj)=0有
特别地,对于同伦(27),我们有
于是,从x(j-1)计算到h(x,λj)=0的解x(j)的牛顿法为：
当‖Δl(j)‖充分小时设定x(j)-xl(j),并进入从x(j)到x(j+1)的计算过程.可见牛顿迭代法中每一步迭代需要解一个关于△l(j)的线性算子方程.
对于同伦(29),有
于是,从x(j-1)计算到h(x,λj)=0的解x(j)的牛顿法为：
‖Δl(j)‖充分小时设定x(j)=xl(j),并进入从x(j)到x(j+1)的计算过程.
从λ=1-δ到λ=1的跟踪过程中,需要求解的线性算子方程变为
是一个不适定线性算子方程,我们可以采用某种正则化策略予以求解,如Tikhonov正则化,谱截断法等.
欧拉法
欧拉法是首先对h(x(λ),λ)=0关于λ求导,得到一个形式上的常微分方程组,
其中x0已知,再对其采用欧拉折线法进行数值求解.特别地,对同伦(27),有
对于给定的点x,上述方程是关于dx的线性算子方程.对λ选取合适的剖分步长,即对0=λ0<λ1<…<λn=1,有
对于同伦(29),有
对于给定的点x,上述方程是关于dx的线性算子方程.对λ选取合适的剖分步长,即对0=λ0<λ1＜…＜λn-1=1-δ,有
从λ=1-δ到λ=1的跟踪过程中,需要求解的线性算子方程变为
是一个不适定线性算子方程,我们采用某种正则化策略予以求解,如Tikhonov正则化,谱截断法等.
在从λ=λj-1到λ=λj的迭代步中,每一步迭代中的初始近似x(j-1)与x(j)都要足够接近才能取得较好的计算效果,这就要λj-λj1足够小,但这样必然会带来较大的计算量,这里给出一种自适应的方法.
第1步：选定阈值ε>0,N为正整数.对λ∈[0,1].给出初始粗网格0=λ0＜λ1＜…＜λn=1,并令j=0;
第2步：如果j=n,终止计算;否则进入第3步;
第3步：计算h(x(λj),λj),如果h(x(λj),λj)＜∈,则令j=j-＞next,返回第2步;否则将[λj-1,λj]再次剖分,将-1/2加入计算网格,形成新的较细的网格
令j=j-＞former,返回第2步;若细分次数大于给定阈值N,终止运算,重新选择初值x0与正则化参数α.
3.4同伦方法在近场光学问题中的应用
本节将同伦方法应用于近场光学问题.我们已经给出近场扫描隧道光学显微镜(PSTM)的数学模型：
以及辐射条件
其中∑R为圆心在原点半径为R的圆周,v是∑R上的单位外法向量.
如果记(32)(33)的解us为S(q),则显然S(q)：L2(D)→L2(R2)是关于q的非线性映射.我们在前一部分针对线性化问题的求解进行了讨论,主要克服其严重不适定性.本节我们讨论处理其非线性的同伦方法.同伦的构造和同伦曲线的跟踪方法已在前节说明,本节主要给出在计算过程中所需要的各阶导数的计算过程.
首先,我们推导S’(q).由于us=S(q)满足(32)(33),故对h∈D(?)Ω(?)R+2, us(q+h)=S(q+h)满足
(34)-(32)并省去高阶项,有
记v=S’(q)h,则v满足
以及辐射条件.上式中右端将n2省略是因为于R+2中n=1(注意h∈D(?)R+2).
由此,可以继续推导二阶导数S”(q).首先假设S”(q)存在.设V1=S’(q+h1)h,其中h1,h∈D(?)Ω(?)R+2,则v1满足
以及辐射条件.(37)-(36)得
利用关系式
舍去高阶项,并令w=S”(q)h1h,有
其中v=S’(q)h1满足
其中w和v均满足辐射条件.
在同伦方法中,需要求解S'(q)*(S(q)-y),其中y是给定测量值.仍记v=S’(q)h1,则其满足(41).在(41)两端乘以(?),并在R2上积分,利用Green公式及辐射条件,有
由关系式
和S'(q)h1=v可知,对给定S(q)-y,求(?)满足
和辐射条件,则
这样,对于带Tikhonov正则化的同伦,即
我们得到h(q(λ),λ)的计算方法：由于(43)(44)式已经给出S'(q)*(S(g)-yδ)的表达式,我们需要按(43)和辐射条件求解微分方程正问题,得到(?),再将(?)代入(44)得到S'(q)*(S(g)-yδ),之后直接运算即得到h(q(λ),λ).
接下来,我们按照自适应跟踪技巧,使用Stefenssen方法进行曲线的跟踪.
第1步：选定阈值∈,∈',N,N',初始网格为0=λ0<λ1<…<λ。=1,并令j=0;
第2步：如果j=n,终止计算;否则进入第3步;
第3步：以qj→former为初值,使用Stefenssen方法计算q(λj).具体地,对l=0,1,…,令
若‖ql+1(j)-ql(j)‖＜∈'且迭代次数l<=N’,则令q(λj)=ql+1(j),否则认为h(q(λj),λj)＞
如果h(q(λj),λj)＜∈,则令j=j→＞next,返回第2步;否则将[λj-1,λj]再次剖分,将λj1/2加入计算网格,形成新的较细的网格
令j=j→＞former,返回第2步;若细分次数大于给定阈值N,终止运算,重新选择初值q0与正则化参数α.
这样,我们得到近场光学问题的同伦算法.
In this paper, the numerical method for solving nonlinear ill-posed problems isproposed, analyzed and applied to the near-field optics. Chapter1is an introduction toscattering problems and inverse scattering problems, the model of near-field optics, theconcept of ill-posed problems and regularization method. Chapter2introduces the mo-ment method with truncated SVD and its applications to severely ill-posed problems,especially in near-field optics. From the numerical experiments, we can conclude thatthemomentmethodwithtruncatedSVDcandealwithseverelyill-posedproblemsmoreeffectively, comparedtotheconventionalregularizationmethodssuchasTikhonovreg-ularization or moments method. Chapter3introduces infinite dimensional homotopymethod to solving nonlinear ill-posed problems. There are two main aspects in ho-motopy method: one is the existence of the homotopy path, the other is the numericalmethod for tracking the homotopy path. In this chapter, two homotopy methods forsolving ill-posed problems are discussed, namely homotopy with Tikhonov regular-ization and homotopy without derivative, and several tracking methods are given forhomotopy path. In the tracking process, adaptive tracking skills are given in order toimprove computing efficiency.
1. Scattering Model of near-field optics
Near-field optics is a new subject since the1980s, breaking the traditional limit ofoptical resolution. There are three important modalities in near-field optics: near-field scanning optical microscopy(NSOM)[23], total internal reflection microscopy(TIRM)[24,25], and photon scanning tunneling microscopy (PSTM)[26,27].
The two principal genres of NSOM are illumination mode and collection mode. In illumination mode NSOM, a probe is scanned over the scatterer in the near zone and the scattering field is measured in the far zone of scatterer. In collection mode NSOM, the scatterer is illuminated by a source located in the far zone of the scatterer and the total field is collected in the near zone of the scatterer through the small aperture in the probe as the probe is scanned over the scatterer. In TIRM, the scatterer is illuminated by an evanescent wave that is generated by total internal reflection, and the scattering field is measured in the far zone of the scatterer. In the PSTM, the scatterer is illuminated by an evanescent wave and the scattering field is measured in the near zone.
In all of the three modalities, reconstruction of the scatterer from the measured data is desirable.
We describe the scattering model mathematically as follows. Assume that the ma-terials are nonmagnetic in the problem considered in this paper. Throughout, we only consider the problem which can be reduced to a two-dimensional model in the case of transverse magnetic(TM) polarization, as is shown in Figure1.
The index of refraction in the lower half-space has a constant value n0. The index of refraction in the upper half-space varies within the domain of the sample but other-wise has a value of unity[19,28]. That is, the background media consists of two layers with the index of refraction where x=(x1,x2), n0＞1.
We denote the scatterer by q(x) within a compact support in D(?)R2+:={x
{x1,x2):x2＞0}, with q(x)/4π the susceptibility of the scatterer [19]. The incident wave comes from the lower half space, which is a time-harmonic plane wave where
Figure1The geometry of the model problem. The scatterer q(x) is in a compact support in D, the measurements are taken on Γ, and is a compact set containing D and Γ. and θ∈(-π/2,π/2), k0is the free space wavenumber. The geometry of the model problem is shown in Figure1. The measurements are taken on Γ:={x=(x1,x2): x2=b＞0,-L＜xi＜L}, and Ω(?)R2is a compact set containing D and Γ.
It is easy to check that the total field u satisfies the Helmholtz equation
Let uref be the reference field, which satisfies the Helmholtz equation without the scat-terer:
From Eq.(2) and the continuity condition on the interface we can analytically obtain
that
where
and are the transmitted and reflected waves respectively[28-30], and
It can be readily seen that the transmitted wave ut becomes evanescent wave when|α|＞k0. That is, ut decays exponentially in the x2direction.
Define the scattering field us=u-uref. Then from (1),(2), we have
And we require the radiation condition [28,31]
where ΣR is the sphere centered in the origin with radius R, and v is the unit outward normal to ΣR.
2. Moment method with truncated SVD and applications
2.1Moment method with truncated SVD
Moment method with truncated SVD is easy to implement. Formally, it is a com-bination of moment method and truncated SVD.
Let Y=C[a, b] and X be Hilbert spaces, the operator K:X→Y be bounded and one-to-one. Let a≤s1＜s2＜…＜sn≤b be the collocation points, and Xn C X be finite-dimensional subspaces with dimXn=n, then the collocation equation is
By Riesz Theorem, there exists kj E X, such that
Kz(sj)=(z,kj),(?)z∈X, j=1,…,n.
Let Xn=span{kj,j=1,…,n}, then the moment solution is given by
where a solves the linear system
with
Then the singular value decomposition(SVD)of A is made by
where U,V are unitary with
and
where λ1≥λ2≥…≥λn≥0.For given K(1≤K≤n),denote α(δ):=∧2≥λK2, Σ=diag(λ1-1,λ2-1,…,λK-1,0,…,0),then the solution of(6)by SVD is
Denote
then xnα(δ),δ is the moment solution with truncated SVD of equation Kx=y.
Now an error estimate of the Moment method with truncated SVD is given as follows:
Theorem1.Let{sn(n),…,sn(n))(?)[a,b],n∈N be a sequence of collocation points, and
Suppose that a solves(6)with β the right-hand side,and assume there exists z∈Cn|z|≤E,such that a=A*z,then
For constant c＞0,choose∧=cδ/E,then
which is optimal,where
2.2Application to linearized near-field optics problem
The linearized near-field optics problem,i.e. inverse problem concerning Born approximation,is to solve the scatterer q(x)from
Note that the only known data lies on the segment Γ={x:x2=b,-L＜x1＜L},so we extend us on Γ to the whole line{x:x2=b,-∞＜x1＜∞}by defining us=0on the remaining interval of Γ.We assume that the error between the expanded data us and the "exact" data is small when L is sufficiently large.Namely,we assume that‖us-uexacts‖＜δ1,where δ1depends on L.The scattering data us depends on α,and we need several different incident waves,i.e.,we need[α1,…,αn]∈[-n0k0,n0k0]to reconstruct the scatterer q(x).
Thus we have
Substituting the fundamental solution G(x,y)and ul into(11),and realizing the depen-dence of us on α,we rewrite(11)as
Multiplying both sides of(12)by e-iωx1and integrating from-∞to∞with respect to x1,
and noting that
(13) can be written as
Realizing that
(14) can be reduced to
When|ω|＜k0and|α|＜k0, both β1and γ are real; and for k0＜|α|＜nok,γ is pure imaginary, and so is β1with|ω|＞k0. It is obvious that, the inversion of (15) involves both inverse Fourier transform and inverse Laplace transform. In fact, both inverse Fourier transform and inverse Laplace transform defined on finite interval are ill posed, especially the inverse Laplace transform is exponentially ill posed(see section2.1.2and [53]).
Here, we split the two-dimensional problem (15) into two one-dimensional prob-lems. So the computation runs much faster and uses less memory. Moreover, we apply the moment method with truncated SVD to get a better resolution where the limitation of numerical precisions is circumvented in a certain extent.
Step1:For ω-α=c(k),k=1,2,…,obtain q(c(k),y2) from a one-dimensional integral equation for every k;
So for every calculation of q(c(k),y2), k=1,2,…, we need special choices of ω That is, if we have n incident angles, then we set ωj=αj+c(k),j=1,…,n, for k=1,2,…. We suppose all the ωj's lie in [-W, W].
Step2:Obtain q(y1,y2) from q(c(k),y2) by inverse Fourier transform on finite intervals.
After the two steps, an approximate solution of the linearized inverse problem (10) is obtained. Obviously, the computation of step1is quite difficult for its severe ill-posedness,especially on a finite data set.As is discussed above,we treat it as a kind of perturbation of the measured data.
The first step is to solve the operator equation
for every specified k,k=1,2,…,where Denote then K:L2[0,b]→L2[-W,W]is a compact operator since k(ω,y2)∈L2([-W,W]×
[0,b]).
Let Xn(?)L2[0,b]be finite-dimensional subspaces with dim Xn=n,and0≤ω1＜…＜ωn≤b be the collocation points,which are chosen as ωj=αj+c(k),j=1,…,n.Then the collocation equation is where
then qn∈Xn is a collocation solution of(18).
Choosing a basis{Φj,j=1,…,n)of Xn,then qn(c(k),y2)can be represented as
where the coefficients aj's are the solution of the system
with
For ill-posed problems such as(19),even small perturbation of data β can change the solution a lot.Recall that we can only obtain the scattering data on Γ,and we have assumed that‖us-uexacts‖＜δ1.Now we denote the perturbed data,or the measured data,by βδ,and let|β-βδ|＜δ,where δtakes in the influence of the error bound δ1, and|·|denotes Euclidean norm in Cn
We will find the moment solution with truncated SVD to(18).By Riesz Theorem,
where is conjugate.If kj,j=1,…,n are linearly independent,i.e.,the determinant of the collocation matrix∧Ai,j=(kj,ki)is not equal to0,by choosing Φj=kj,there exists one and only one moment solution qn of(18).
Using SVD,we can write
where λ1≥…≥λn≥0,let α(δ):=∧2≥λK2for some K(1≤K≤n), Σ=diag(λ1-1,λ2-1,…,λK-1,0,…,0),then
So the moment solution with truncated SVD to(18)is given by
By using Theorem1,we have the following estimate for linearized near-field optics problem:
Theorem2.Let{ω1(n),ω2(n),…,ωn(n)}(?)[0,b],n∈N,be a sequence of collocation points and
And assume that there exist z∈Cn,|z|≤E,such that a=A*z.Let qn～δ be the moment solution with truncated SVD to(18),then the following error estimate holds:
For some constant c＞0,choose∧=cδ/E,then
where
By our truncated SVD method,an appropriate∧could be chosen so that δ/∧is bounded.In practice,we could also choose∧=δα,0＜α＜1,then δ/∧+∧｜z δ1-α+cδα is independent of n. ically by inverse Fourier transform. By Parseval formula,it is obvious that‖qn～δ q‖L22(D)＝‖qn～δ-q‖L22(D)/(2π).
Numerical experiments show that this method computes fast,stable,and can get satisfying resolution even under the perturbed data.
3.Homotopy method for infinite dimensional nonlinear ill-posed problems and applications
3.1Homotopy with Tikhonov regularization
We combine homotopy method in infinite dimensional spaces with Tikhonov reg-ularization to solve infinite dimensional ill-posed problems.
Let X,Y be Hilbert spaces,consider nonlinear ill-posed problem
where F:D(F)(?)X→Y is nonlinear ill-posed(compact)operator.Construct minimization funcitional
where x0∈D(F),A(λ):X→X is an operator dependent to λ∈[0,1].
Remark1.If we set A(λ)=1-(1-α)λI,the when λ=1,(25)is nothing but Tikhonov functional.If the homotopy path exists,the Tikhonov regularization solution to(24)is obtained.
Assume that for every x∈D(F),F is Frechet differentiable,then for every h∈X,
J is minimized when Re(AF'(x)*(F(x)-yδ)+A*(λ)A(λ)(x-x0))=0.
Thus,the normal equation of(25)is given by:
Remark2.For the case when A(λ)=1-(1-α)λI, which is just the fixed-point homotopy with Tikhonov regularization of(24).
For the case when A(λ)=1-(1-α)λI,,the fixed-point homotopy with Tikhonov regularization of(24)is obtained:
Theorem3.LET D(F)be a bounded connected open set,F(x)∈C2(D(F),Y),F weakly closed, and there exists an element x∈D(F),such that F(x)=y.For x0∈D(F),denote h(x,t)as in(26),where A(λ)=1-(1-α)λI.If h(x,λ)≠0when(x,λ)∈(?)D(F)×[0,1),then there exists at least one path in the solution set of h(x,λ)=0,connecting x0to one regularization solution to(24).
In the tracking of the homotopy path,the Seidel skill can be adopted,i.e.,using the latest x(λ)instead of x0in(27)in order to obtain better solution.
3.2Homotopy without derivative
Homotopy with Tikhonov regularization contains the derivative of nonlinear map-ping F,moreover, in the tracking process,the second order derivative of F may have to be calculated,which makes the computing rather difficult. So homotopy without derivative is desirable.
Construct formally the fixed-point homotopy of(24):
wheye0＜ρ＜1.Obviously,for λ∈[0,1-ρ],the problem of solving h(x,λ)=0in
(29)is well posed.
We can obtain the following result similar to Theorem3:
Theorem4.Let D(F)be abounded connected open set,F(x)∈C2(D(F),Y).For x0∈D(F),denote h(x,t)as in(29).If h(x,λ)≠0when(x,λ)∈(?)D(F)×[0,1), then there exists at least one path in the solution set of h(x,λ)＝0,connecting x0to x1-ρ.
If the solution x1-ρ of h(x,1-ρ)=0is sufficiently close to the solution of F(x)=y,we can take x1-ρ as the initial approximation,and using some regularization method to solve(24).
3.3Tracking the homotopy path
Numerically,we need appropriate method to track the homotopy path.Assume that x can be parameterized by λ,i.e.,the homotopy path can be written as h(x(λ),λ)=0, with x(0)the solution of h(x,0)=0,which is given.For0=λ0＜λ1＜…＜λn=1, we need to solve h(x,λj)=0,j=1,…,n one by one,whose solutions are denoted by x(j),j=1,…,n.We introduce some methods such as Stefenssen method,Newton method,and Euler method.
Stefenssen method
Stefenssen method means solving h(x,λj)=0sequentially by j=1,…,n, whose solutions are denoted by x(j),so x(j-1)could be an initial approximation.Namely, for l=0,1,…,let
In each iteration,the initial approximation x(j-1)should be close to x(j)enough, but for ill-posed problems,we have to make λj-λj1small enough to make it the case,which brings over a large amount of computation,so we will adopt an adaptive approach to deal with this problem.
Adaptive tracking skills
As is discussed,we have to make λj-λj1small enough to make the initial ap-proximation x(j-1)close to x(j),but this brings over a large amount of computation.We introduce the adaptive tracking skills as follows.
Step1:Select threshold∈＞0,N a positive integer.For λ∈[0,1],give initial coarse grid0=λ0＜λ1＜…＜λn=1,and let j=0;
Step2:When j=n,terminate the calculation;otherwise go to Step3;
Step3:Calculate h(x(λj),λj),if h(x(λj),λj)＜∈,then let j=j-＞next,return to step2;otherwise split[λj-1,λj],and put λj-1/2into the new computing grid,i.e.,a new finer grid is formed as
then let j=j-＞former,return to step2;If the subdivision number is greater than the given threshold N,terminate the calculation,and re-select the initial x0and regularization parameter α.
Newton Method
Similar to Stefenssen method,we assume that the solution to h(x,λj-1)=0,i.e., x(j-1)has been obtained,and use it as an initial approximation to solve h(x,λj)=0. Here we use the local convergent Newton's method.Let
In particular, for the homotopy(27),
Thus,the Newton's method for solving h(x,λj)=0from the initial approximation
x(j-1)can be written as:
When‖△l(j)‖is sufficiently small we set x(j)=xl(j),and compute x(j+1)from the initial approximation x(j).It is readily that every iteration in Newton's method needs to solve a linear operator equation with respect to△l(j)
For homotopy(29),
Thus,the Newton's method for solving h(x,λj)=0from the initial approximation x(j-1)can be written as:
When‖△l(j)‖is sufficiently small we set x(j)=xl(j), and compute x(j+1)from the initial approxim ation x(j).
From λ=1-δ to λ=1,the linear operator equation becomes
which are ill-posed linear operator equations,so we can solve them by some kind of regularization method,such as Tikhonov regularization,SVD method.
Euler method
By Euler method,we derivative h(x(λ),λ)=0with respect to λ,and obtain a
where x0is given. Then solve it by Euler method numerically.
In particular, for the homotopy (27), we have
For a fixed point x, the above equation is a linear operator equation with respect to dx. For0=λ0＜λ1＜…＜λn=1, we have
For the homotopy (29), we have
For a fixed point x, the above equation is linear operator equation with respect to dx. For0=λ0＜λ1＜…＜λn=1, we have
From λ=1-δ to λ=1, the linear operator equation becomes
which is an ill-posed linear operator equation, so we can solve it by some kind of regu-larization method, such as Tikhonov regularization, SVD method.
3.4Application to near-field optics
In this section we will apply the homotopy method to near-field optics. The model of photon scanning tunneling microscopy (PSTM) has been given by
with radiation condition
where ΣR is the sphere centered in the origin with radius R,and v is the unit outward normal to∑R.
Denote the solution to(32)(33),i.e.us by S(q),then it is readily seen that S(q) L2(D)→L2(R2)is nonlinear in q.
In the previous part we solved the problem for the linearized case,mainly to over-come its severe ill-posedness.In this section we discuss the homotopy method to deal with the nonlinearity.Homotopy method and homotopy path have been given in the previous section,so in this section we calculate the derivatives in need.
First,we derive S'(q).
since us=S(q)satisfies(32)(33),so for h∈D(?)Ω(?)R+2,us(q+h)=S(q+h) satisfies
(34)-(32)and omit higher order terms,we get
△(S'(q)h)+n22k02(1+q)(S'(q)h)+n2k02hS(q)=+k02hut,in R+2.(35) Denote v=S'(q)h,then v satisfies
△v+n2k02(1+q)v=-k02hut-n2k2hS(q)=-k02h(ut+us(q)),in R2+,(36)
and radiation condition.
Thus,we can continue to derive the second derivative.Assume that S"(q)exists. Let v1=S'(q+h1)h,where h1,h∈D(?)Ω(?)R+2,then v1satisfies
and radiation condition.(37)-(36)and we get
Using the relationship
and omitting higher order terms, and denoting ω=S"(q)h1h, we have
where v=S'(q)h1satisfies
where both w and v satisfy radiation condition.
In homotopy method, we need to solve S'(q)*(S(q)-y), where y is measured data. Again, denote v=S'(q)h1, then it satisfies (41). Multiply0on both sides of (41), and integrate on R2, by Green's formula and radiation condition,
By the relationship
and S'(q)h1=v, we can get that, for given S(q)-y, solve Φ that satisfies
and radiation condition, then
Thus, for homotopy with Tikhonov regularization, i.e.,
the numerical method for computing h(q(λ),λ) has been obtained:since the formula of S'(q)*(S(q)-yδ) has been given in (43) and (44), we solve0from PDE (43) and radiation condition, then substitute0into (44) and get S'(q)*(S(q)-yδ), and then h(q(λ),λ).
Next, by the adaptive tracking skills, we use Stefenssen method to track the homo-topy path.
Step1:Select threshold∈,∈',and N,N' positive integers,and the initial coarse grid0=λ0＞λ1＜…＜λn=1,and let j=0;
Step2:When j=n,terminate the calculation;otherwise go to Step3;
Step3:Calculate q(λj)by Stefenssen method,using qj-＞former as initial guess, namely,for l=0,1,…,let
if‖ql+1(j)-ql(j)‖＜∈' and the iteration steps l＜=N',then let q(λj)=ql-1(j),otherwise let h(q(λj),λj)＞∈.
If h(q(λj),λj)＜∈,then let j=j-＞nexl,and return to step2;otherwise split [λj-1,λj],and put λj-1/2into the new computing grid,i.e.,a new finer grid is formed as
then let j=j-＞former,and return to step2;If the subdivision number is greater than the given threshold N,terminate the calculation,and re-select the initial q0and regularization parameter α.
In this way,the homotopy algorithm for near-field optics is obtained.
引文
[1] COLTOND,KRESSR.AppliedMathematicalSciences, Vol93: Inverseacousticand electromagnetic scattering theory[M/OL]. Third.[S.l.]: Springer, New York,2013: xiv+405. http://dx.doi.org/10.1007/978-1-4614-4942-3.
[2] LI P. Numerical simulations of global approach for photon scanning tunnelingmicroscopy: coupling of finite-element and boundary integral methods[J]. JOSAA,2008,25(8):1929–1936.
[3] LI P. Coupling of finite element and boundary integral methods for electromag-netic scattering in a two-layered medium[J]. Journal of Computational Physics,2010,229(2):481–497.
[4] BERENGER J-P. A perfectly matched layer for the absorption of electromagneticwaves[J]. Journal of computational physics,1994,114(2):185–200.
[5] TURKEL E, YEFET A. Absorbing PML boundary layers for wave-like equa-tions[J/OL]. Appl. Numer. Math.,1998,27(4):533–557. http://dx.doi.org/10.1016/S0168-9274(98)00026-9.
[6] RAPPAPORTCM,WINTONS.UsingthePMLABCforair/soilwaveinteractionmodeling in the time and frequency domains[J]. Subsurface Sensing Technologiesand Applications,2000,1(3):289–303.
[7] TSYNKOVS,TURKELE.ACartesianperfectlymatchedlayerfortheHelmholtzequation[G]//ANON. Absorbing boundaries and layers, domain decompositionmethods. Huntington, NY: Nova Sci. Publ.,2001:279–309.
[8] CESSENAT O, DESPRES B. Application of an ultra weak variational for-mulation of elliptic PDEs to the two-dimensional Helmholtz problem[J/OL].SIAM J. Numer. Anal.,1998,35(1):255–299. http://dx.doi.org/10.1137/S0036142995285873.
[9] CESSENAT O, DESPRéS B. Using plane waves as base functions for solv-ing time harmonic equations with the ultra weak variational formulation[J/OL].J. Comput. Acoust.,2003,11(2):227–238. http://dx.doi.org/10.1142/S0218396X03001912.
[10] HUTTUNEN T, MONK P, KAIPIO J P. Computational aspects of the ultra-weakvariational formulation[J/OL]. J. Comput. Phys.,2002,182(1):27–46. http://dx.doi.org/10.1006/jcph.2002.7148.
[11] HUTTUNEN T, KAIPIO J P, MONK P. The perfectly matched layer for the ul-tra weak variational formulation of the3D Helmholtz equation[J/OL]. Internat. J.Numer. Methods Engrg.,2004,61(7):1072–1092. http://dx.doi.org/10.1002/nme.1105.
[12] BUFFA A, MONK P. Error estimates for the ultra weak variational formulation ofthe Helmholtz equation[J/OL]. M2AN Math. Model. Numer. Anal.,2008,42(6):925–940. http://dx.doi.org/10.1051/m2an:2008033.
[13] LUOSTARIT, HUTTUNENT, MONKP.Improvementsfortheultraweakvaria-tional formulation[J/OL]. Internat. J. Numer. Methods Engrg.,2013,94(6):598–624. http://dx.doi.org/10.1002/nme.4469.
[14] LUOSTARI T, HUTTUNEN T, MONK P. The ultra weak variational formulationusing Bessel basis functions[J/OL]. Commun. Comput. Phys.,2012,11(2):400–414. http://dx.doi.org/10.4208/cicp.121209.040111s.
[15] COLTOND,KIRSCHA.Asimplemethodforsolvinginversescatteringproblemsin the resonance region[J/OL]. Inverse Problems,1996,12(4):383–393. http://dx.doi.org/10.1088/0266-5611/12/4/003.
[16] CAKONI F, COLTON D, MONK P. CBMS-NSF Regional Conference Series inApplied Mathematics, Vol80: The linear sampling method in inverse electro-magnetic scattering[M/OL].[S.l.]: Society for Industrial and Applied Mathemat-ics (SIAM), Philadelphia, PA,2011: x+138. http://dx.doi.org/10.1137/1.9780898719406.
[17]KIRSCH A, GRINBERG N. Oxford Lecture Series in Mathematics and its Appli-cations, Vol36:The factorization method for inverse problems[M].[S.1.]:Oxford University Press, Oxford,2008:xiv+201.
[18]KIRSCH A. The factorization method for Maxwell's equations[J/OL]. In-verse Problems,2004,20(6):S117-S134. http://dx.doi.org/10.1088/0266-5611/20/6/S08.
[19]CARNEY P S, SCHOTLAND J C. Near-field tomography[J]. Inside out:inverse problems and applications,2003,47:133-168.
[20]COURJON D, BAINIER C. Near field microscopy and near field optics[J]. Re-ports on progress in Physics,1994,57(10):989.
[21]朱星.近场光学与近场光学显微镜[J].北京大学学报(自然科学版),1997,33(3)：394-407.
[22]祝生祥.传统光学显微镜与近场光学显微镜[J].光学仪器,2000,22(6)：34-41.
[23]CARNEY P S, SCHOTLAND J C. Inverse scattering for near-field microscopy [J]. Applied Physics Letters,2000,77(18):2798-2800.
[24]CARNEY P S, SCHOTLAND J C. Three-dimensional total internal reflection mi-croscopy [J]. Optics Letters,2001,26(14):1072-1074.
[25]CARNEY P S, SCHOTLAND J C. Theory of total-internal-reflection tomogra-phy [J]. JOSA A,2003,20(3):542-547.
[26]CARNEY P S, SCHOTLAND J C. Determination of three-dimensional structure in photon scanning tunnelling microscopy [J]. Journal of Optics A:Pure and Ap-plied Optics,2002,4(5):S140.
[27]COURJON D, SARAYEDDINE K, SPAJER M. Scanning tunneling optical mi-croscopy [J]. Optics communications,1989,71(1):23-28.
[28] BAO G, LI P. Inverse medium scattering problems in near-field optics[J]. Journalof Computational Mathematics-International Edition-,2007,25(3):252.
[29] AKDUMAN I, ALKUMRU A. A generalized ART algorithm for inverse scat-tering problems related to buried cylindrical bodies[J]. Inverse problems,1999,11(6):1125.
[30] IDEMEN M, AKDUMAN I. Two-dimensional inverse scattering problems con-nected with bodies buried in a slab[J]. Inverse problems,1999,6(5):749.
[31] XU Y. Radiation condition and scattering problem for time-harmonic acousticwaves in a stratified medium with a nonstratified inhomogeneity[J]. IMA journalof applied mathematics,1995,54(1):9–29.
[32] KRISTENSSON G. A uniqueness theorem for Helmholtz’equation: penetrablemedia with an infinite interface[J]. SIAM Journal on Mathematical Analysis,1980,11(6):1104–1117.
[33] ZHANGB, CHANDLER-WILDESN.Acousticscatteringbyaninhomogeneouslayer on a rigid plate[J]. SIAM Journal on Applied Mathematics,1998,58(6):1931–1950.
[34] CUTZACH P-M, HAZARD C. Existence, uniqueness and analyticity propertiesforelectromagneticscatteringinatwo-layeredmedium[J].Mathematicalmethodsin the applied sciences,1998,21(5):433–461.
[35] DURáN M, MUGA I, NéDéLEC J-C. The Helmholtz equation with impedancein a half-plane[J]. Comptes Rendus Mathematique,2005,340(7):483–488.
[36] CHENY. Inversescatteringvia Heisenberg’s uncertaintyprinciple[J/OL]. InverseProblems,1997,13(2):253–282.http://dx.doi.org/10.1088/0266-5611/13/2/005.
[37] BAO G, LI P. Inverse medium scattering for the Helmholtz equation at fixed fre-quency[J]. Inverse problems,2005,21(5):1621.
[38] BAO G, LI P. Inverse medium scattering problems for electromagneticwaves[J/OL].SIAMJ.Appl.Math.,2005,65(6):2049–2066(electronic).http://dx.doi.org/10.1137/040607435.
[39] BAO G, CHEN Y, MA F. Regularity and stability for the scattering map of alinearized inverse medium problem[J]. Journal of mathematical analysis and ap-plications,2000,247(1):255–271.
[40] SCHOTLAND J C, MARKEL V A. Inverse scattering with diffusingwaves[J/OL].J.Opt.Soc.Amer.A,2001,18(11):2767–2777.http://dx.doi.org/10.1364/JOSAA.18.002767.
[41] MOSKOW S, SCHOTLAND J C. Convergence and stability of the inverse scat-tering series for diffuse waves[J/OL]. Inverse Problems,2008,24(6):065005,16.http://dx.doi.org/10.1088/0266-5611/24/6/065005.
[42] ARRIDGE S, MOSKOW S, SCHOTLAND J C. Inverse Born series for theCalderon problem[J/OL]. Inverse Problems,2012,28(3):035003,16. http://dx.doi.org/10.1088/0266-5611/28/3/035003.
[43] MARKEL V, MITAL V, SCHOTLAND J. Inverse problem in optical diffusiontomography. III. Inversion formulas and singular-value decomposition[J]. JOSAA,2003,20(5):890–902.
[44] HOFMANN B, KINDERMANN S, OTHERS. On the degree of ill-posedness forlinearproblemswithnoncompactoperators[J].MethodsandApplicationsofAnal-ysis,2010,17(4):445–462.
[45] HOFMANN B. Mathematik inverser Probleme[M].[S.l.]: Teubner Stuttgart,1999.
[46] KIRSCH A. Applied Mathematical Sciences, Vol120: An introduction to themathematical theory of inverse problems[M/OL]. New York: Springer-Verlag,1996: x+282. http://dx.doi.org/10.1007/978-1-4612-5338-9.
[47]MOROZO V V A. Regularization methods for ill-posed problems[M]. Boca Raton, FL:CRC Press,1993:vi+257.
[48]BAKUSHINSKY A, GONCHARSKY A. Mathematics and its Applications, Vol301:Ill-posed problems:theory and applications[M/OL]. Dordrecht:Kluwer Academic Publishers Group,1994:x+256. http://dx.doi.org/10.1007/978-94-011-1026-6.
[49]刘继军.不适定问题的正则化方法及应用[M].[S.1.]：北京：科学出版社,2005.
[50]韩波,李莉.非线性不适定问题的求解方法及其应用[M].[S.1.]：北京：科学出版社,2011.
[51]LU S, PEREVERZEV S V, RAMLAU R. An analysis of Tikhonov regular-ization for nonlinear ill-posed problems under a general smoothness assump-tion[J/OL]. Inverse Problems,2007,23(1):217-230. http://dx.doi. org/10.1088/0266-5611/23/1/011.
[52]ENGL H W, KUNISCH K, NEUB AUER A. Convergence rates for Tikhonov reg-ularisation of non-linear ill-posed problems[J]. Inverse problems,1989,5(4):523.
[53]EPSTEIN C L, SCHOTLAND J. The bad truth about Laplace's transform[J]. SIAM review,2008,50(3):504-520.
[54]D'AMORE L, LACCETTI G, MURLI A. An implementation of a Fourier series method for the numerical inversion of the Laplace transform [J]. ACM Transac-tions on Mathematical Software (TOMS),1999,25(3):279-305.
[55]ABATE J, VALKO P. Multi-precision Laplace transform inversion[J]. Interna-tional Journal for Numerical Methods in Engineering,2004,60(5):979-993.
[56]ABATE J, CHOUDHURY G L, WHITT W. On the Laguerre method for numer-ically inverting Laplace transforms [J]. INFORMS Journal on Computing,1996,8(4):413-427.
[57]VALKO P P, ABATE J. Comparison of sequence accelerators forthe Gaver method of numerical Laplace transform inversion[J]. Computers&Mathematics with Ap-plications,2004,48(3):629-636.
[58]WIMP J, WIMP J. Sequence transformations and their applications[M].[S.1.]: Academic Press New York,1981.
[59]STEHFEST H. Algorithm368:Numerical inversion of Laplace transforms [D5][J]. Communications of the ACM,1970,13(1):47-49.
[60]THAMBAN NAIR M, PEREVERZEV S. Regularized collocation method for Fredholm integral equations of the first kind[J]. Journal of Complexity,2007,23(4):454-467.
[61]江泽坚,孙善利.泛函分析[M].[S.1.]：北京：高等教育出版社,1994.
[62]张恭庆,林源渠.泛函分析讲义(上册)[M].[S.1.]：北京大学出版社,1987.
[63]CHANG K-C. Springer Monographs in Mathematics:Methods in nonlinear anal-ysis[M]. Berlin:Springer-Verlag,2005:x+439.
[64]ZEIDLER E. Nonlinear functional analysis and its applications. Ⅰ[M/OL]. New York:Springer-Verlag,1986:xxi+897. http://dx.doi.org/10.1007/978-1-4612-4838-5.
[65]ZEIDLER E. Nonlinear functional analysis and its applications. Ⅱ/A[M/OL]. New York:Springer-Verlag,1990:xviii+467. http://dx.doi.org/10.1007/978-1-4612-0985-0.
[66]ZEIDLER E. Nonlinear functional analysis and its applications. Ⅱ/B[M/OL]. New York:Springer-Verlag,1990:i-1202. http://dx.doi.org/10.1007/978-1-4612-0985-0.
[67]ZEIDLER E. Nonlinear functional analysis and its applications. Ⅲ[M]. New York:Springer-Verlag,1985:xxii+662.
[68]ZEIDLER E. Nonlinear functional analysis and its applications. IV[M]. New York:Springer-Verlag,1988:xxiv+975.
[69]DEIMLING K. Nonlinear functional analysis[M]. Berlin:Springer-Verlag,1985: xiv+450.
[70]SMALE S. An infinite dimensional version of Sard's theorem[J]. American Jour-nal of Mathematics,1965,87(4):861-866.
[71]SCHROTER T. On a special class of nonlinear Fredholm integral equations of the first kind[J]. Computing,1997,58(3):259-279.
[72]钟承奎,范先令,陈文塬.非线性泛函分析引论[M].[S.1.]：兰州：兰州大学出版社,2004.
[73]HRYCAK T, ISAKOV V. Increased stability in the continuation of solutions to the Helmholtz equation[J/OL]. Inverse Problems,2004,20(3):697-712. http://dx.doi.org/10.1088/0266-5611/20/3/004.
[74]MAKINO M, OISHI S. Constructive analysis for infinite-dimensional nonlin-ear systems-Infinite-dimensional version of homotopy method[J]. Electronics and Communications in Japan (Part III:Fundamental Electronic Science),1991,74(8):1-10.
[75]DEIMLING K. Nonlinear functional analysis[M].[S.1.]:Courier Dover Publica-tions,2011.
[76]HAN B, FU H, LIU H. A homotopy method for well-log constraint waveform inversion[J]. Geophysics,2006,72(1):R1-R7.
[77]ZHAO J, LIU T, LIU S. An adaptive homotopy method for permeability estima-tion of a nonlinear diffusion equation[J]. Inverse Problems in Science and Engi-neering,2013,21(4):585-604.
[2] LI P. Numerical simulations of global approach for photon scanning tunnelingmicroscopy: coupling of finite-element and boundary integral methods[J]. JOSAA,2008,25(8):1929–1936.
[3] LI P. Coupling of finite element and boundary integral methods for electromag-netic scattering in a two-layered medium[J]. Journal of Computational Physics,2010,229(2):481–497.
[4] BERENGER J-P. A perfectly matched layer for the absorption of electromagneticwaves[J]. Journal of computational physics,1994,114(2):185–200.
[5] TURKEL E, YEFET A. Absorbing PML boundary layers for wave-like equa-tions[J/OL]. Appl. Numer. Math.,1998,27(4):533–557. http://dx.doi.org/10.1016/S0168-9274(98)00026-9.
[6] RAPPAPORTCM,WINTONS.UsingthePMLABCforair/soilwaveinteractionmodeling in the time and frequency domains[J]. Subsurface Sensing Technologiesand Applications,2000,1(3):289–303.
[7] TSYNKOVS,TURKELE.ACartesianperfectlymatchedlayerfortheHelmholtzequation[G]//ANON. Absorbing boundaries and layers, domain decompositionmethods. Huntington, NY: Nova Sci. Publ.,2001:279–309.
[8] CESSENAT O, DESPRES B. Application of an ultra weak variational for-mulation of elliptic PDEs to the two-dimensional Helmholtz problem[J/OL].SIAM J. Numer. Anal.,1998,35(1):255–299. http://dx.doi.org/10.1137/S0036142995285873.
[9] CESSENAT O, DESPRéS B. Using plane waves as base functions for solv-ing time harmonic equations with the ultra weak variational formulation[J/OL].J. Comput. Acoust.,2003,11(2):227–238. http://dx.doi.org/10.1142/S0218396X03001912.
[10] HUTTUNEN T, MONK P, KAIPIO J P. Computational aspects of the ultra-weakvariational formulation[J/OL]. J. Comput. Phys.,2002,182(1):27–46. http://dx.doi.org/10.1006/jcph.2002.7148.
[11] HUTTUNEN T, KAIPIO J P, MONK P. The perfectly matched layer for the ul-tra weak variational formulation of the3D Helmholtz equation[J/OL]. Internat. J.Numer. Methods Engrg.,2004,61(7):1072–1092. http://dx.doi.org/10.1002/nme.1105.
[12] BUFFA A, MONK P. Error estimates for the ultra weak variational formulation ofthe Helmholtz equation[J/OL]. M2AN Math. Model. Numer. Anal.,2008,42(6):925–940. http://dx.doi.org/10.1051/m2an:2008033.
[13] LUOSTARIT, HUTTUNENT, MONKP.Improvementsfortheultraweakvaria-tional formulation[J/OL]. Internat. J. Numer. Methods Engrg.,2013,94(6):598–624. http://dx.doi.org/10.1002/nme.4469.
[14] LUOSTARI T, HUTTUNEN T, MONK P. The ultra weak variational formulationusing Bessel basis functions[J/OL]. Commun. Comput. Phys.,2012,11(2):400–414. http://dx.doi.org/10.4208/cicp.121209.040111s.
[15] COLTOND,KIRSCHA.Asimplemethodforsolvinginversescatteringproblemsin the resonance region[J/OL]. Inverse Problems,1996,12(4):383–393. http://dx.doi.org/10.1088/0266-5611/12/4/003.
[16] CAKONI F, COLTON D, MONK P. CBMS-NSF Regional Conference Series inApplied Mathematics, Vol80: The linear sampling method in inverse electro-magnetic scattering[M/OL].[S.l.]: Society for Industrial and Applied Mathemat-ics (SIAM), Philadelphia, PA,2011: x+138. http://dx.doi.org/10.1137/1.9780898719406.
[17]KIRSCH A, GRINBERG N. Oxford Lecture Series in Mathematics and its Appli-cations, Vol36:The factorization method for inverse problems[M].[S.1.]:Oxford University Press, Oxford,2008:xiv+201.
[18]KIRSCH A. The factorization method for Maxwell's equations[J/OL]. In-verse Problems,2004,20(6):S117-S134. http://dx.doi.org/10.1088/0266-5611/20/6/S08.
[19]CARNEY P S, SCHOTLAND J C. Near-field tomography[J]. Inside out:inverse problems and applications,2003,47:133-168.
[20]COURJON D, BAINIER C. Near field microscopy and near field optics[J]. Re-ports on progress in Physics,1994,57(10):989.
[21]朱星.近场光学与近场光学显微镜[J].北京大学学报(自然科学版),1997,33(3)：394-407.
[22]祝生祥.传统光学显微镜与近场光学显微镜[J].光学仪器,2000,22(6)：34-41.
[23]CARNEY P S, SCHOTLAND J C. Inverse scattering for near-field microscopy [J]. Applied Physics Letters,2000,77(18):2798-2800.
[24]CARNEY P S, SCHOTLAND J C. Three-dimensional total internal reflection mi-croscopy [J]. Optics Letters,2001,26(14):1072-1074.
[25]CARNEY P S, SCHOTLAND J C. Theory of total-internal-reflection tomogra-phy [J]. JOSA A,2003,20(3):542-547.
[26]CARNEY P S, SCHOTLAND J C. Determination of three-dimensional structure in photon scanning tunnelling microscopy [J]. Journal of Optics A:Pure and Ap-plied Optics,2002,4(5):S140.
[27]COURJON D, SARAYEDDINE K, SPAJER M. Scanning tunneling optical mi-croscopy [J]. Optics communications,1989,71(1):23-28.
[28] BAO G, LI P. Inverse medium scattering problems in near-field optics[J]. Journalof Computational Mathematics-International Edition-,2007,25(3):252.
[29] AKDUMAN I, ALKUMRU A. A generalized ART algorithm for inverse scat-tering problems related to buried cylindrical bodies[J]. Inverse problems,1999,11(6):1125.
[30] IDEMEN M, AKDUMAN I. Two-dimensional inverse scattering problems con-nected with bodies buried in a slab[J]. Inverse problems,1999,6(5):749.
[31] XU Y. Radiation condition and scattering problem for time-harmonic acousticwaves in a stratified medium with a nonstratified inhomogeneity[J]. IMA journalof applied mathematics,1995,54(1):9–29.
[32] KRISTENSSON G. A uniqueness theorem for Helmholtz’equation: penetrablemedia with an infinite interface[J]. SIAM Journal on Mathematical Analysis,1980,11(6):1104–1117.
[33] ZHANGB, CHANDLER-WILDESN.Acousticscatteringbyaninhomogeneouslayer on a rigid plate[J]. SIAM Journal on Applied Mathematics,1998,58(6):1931–1950.
[34] CUTZACH P-M, HAZARD C. Existence, uniqueness and analyticity propertiesforelectromagneticscatteringinatwo-layeredmedium[J].Mathematicalmethodsin the applied sciences,1998,21(5):433–461.
[35] DURáN M, MUGA I, NéDéLEC J-C. The Helmholtz equation with impedancein a half-plane[J]. Comptes Rendus Mathematique,2005,340(7):483–488.
[36] CHENY. Inversescatteringvia Heisenberg’s uncertaintyprinciple[J/OL]. InverseProblems,1997,13(2):253–282.http://dx.doi.org/10.1088/0266-5611/13/2/005.
[37] BAO G, LI P. Inverse medium scattering for the Helmholtz equation at fixed fre-quency[J]. Inverse problems,2005,21(5):1621.
[38] BAO G, LI P. Inverse medium scattering problems for electromagneticwaves[J/OL].SIAMJ.Appl.Math.,2005,65(6):2049–2066(electronic).http://dx.doi.org/10.1137/040607435.
[39] BAO G, CHEN Y, MA F. Regularity and stability for the scattering map of alinearized inverse medium problem[J]. Journal of mathematical analysis and ap-plications,2000,247(1):255–271.
[40] SCHOTLAND J C, MARKEL V A. Inverse scattering with diffusingwaves[J/OL].J.Opt.Soc.Amer.A,2001,18(11):2767–2777.http://dx.doi.org/10.1364/JOSAA.18.002767.
[41] MOSKOW S, SCHOTLAND J C. Convergence and stability of the inverse scat-tering series for diffuse waves[J/OL]. Inverse Problems,2008,24(6):065005,16.http://dx.doi.org/10.1088/0266-5611/24/6/065005.
[42] ARRIDGE S, MOSKOW S, SCHOTLAND J C. Inverse Born series for theCalderon problem[J/OL]. Inverse Problems,2012,28(3):035003,16. http://dx.doi.org/10.1088/0266-5611/28/3/035003.
[43] MARKEL V, MITAL V, SCHOTLAND J. Inverse problem in optical diffusiontomography. III. Inversion formulas and singular-value decomposition[J]. JOSAA,2003,20(5):890–902.
[44] HOFMANN B, KINDERMANN S, OTHERS. On the degree of ill-posedness forlinearproblemswithnoncompactoperators[J].MethodsandApplicationsofAnal-ysis,2010,17(4):445–462.
[45] HOFMANN B. Mathematik inverser Probleme[M].[S.l.]: Teubner Stuttgart,1999.
[46] KIRSCH A. Applied Mathematical Sciences, Vol120: An introduction to themathematical theory of inverse problems[M/OL]. New York: Springer-Verlag,1996: x+282. http://dx.doi.org/10.1007/978-1-4612-5338-9.
[47]MOROZO V V A. Regularization methods for ill-posed problems[M]. Boca Raton, FL:CRC Press,1993:vi+257.
[48]BAKUSHINSKY A, GONCHARSKY A. Mathematics and its Applications, Vol301:Ill-posed problems:theory and applications[M/OL]. Dordrecht:Kluwer Academic Publishers Group,1994:x+256. http://dx.doi.org/10.1007/978-94-011-1026-6.
[49]刘继军.不适定问题的正则化方法及应用[M].[S.1.]：北京：科学出版社,2005.
[50]韩波,李莉.非线性不适定问题的求解方法及其应用[M].[S.1.]：北京：科学出版社,2011.
[51]LU S, PEREVERZEV S V, RAMLAU R. An analysis of Tikhonov regular-ization for nonlinear ill-posed problems under a general smoothness assump-tion[J/OL]. Inverse Problems,2007,23(1):217-230. http://dx.doi. org/10.1088/0266-5611/23/1/011.
[52]ENGL H W, KUNISCH K, NEUB AUER A. Convergence rates for Tikhonov reg-ularisation of non-linear ill-posed problems[J]. Inverse problems,1989,5(4):523.
[53]EPSTEIN C L, SCHOTLAND J. The bad truth about Laplace's transform[J]. SIAM review,2008,50(3):504-520.
[54]D'AMORE L, LACCETTI G, MURLI A. An implementation of a Fourier series method for the numerical inversion of the Laplace transform [J]. ACM Transac-tions on Mathematical Software (TOMS),1999,25(3):279-305.
[55]ABATE J, VALKO P. Multi-precision Laplace transform inversion[J]. Interna-tional Journal for Numerical Methods in Engineering,2004,60(5):979-993.
[56]ABATE J, CHOUDHURY G L, WHITT W. On the Laguerre method for numer-ically inverting Laplace transforms [J]. INFORMS Journal on Computing,1996,8(4):413-427.
[57]VALKO P P, ABATE J. Comparison of sequence accelerators forthe Gaver method of numerical Laplace transform inversion[J]. Computers&Mathematics with Ap-plications,2004,48(3):629-636.
[58]WIMP J, WIMP J. Sequence transformations and their applications[M].[S.1.]: Academic Press New York,1981.
[59]STEHFEST H. Algorithm368:Numerical inversion of Laplace transforms [D5][J]. Communications of the ACM,1970,13(1):47-49.
[60]THAMBAN NAIR M, PEREVERZEV S. Regularized collocation method for Fredholm integral equations of the first kind[J]. Journal of Complexity,2007,23(4):454-467.
[61]江泽坚,孙善利.泛函分析[M].[S.1.]：北京：高等教育出版社,1994.
[62]张恭庆,林源渠.泛函分析讲义(上册)[M].[S.1.]：北京大学出版社,1987.
[63]CHANG K-C. Springer Monographs in Mathematics:Methods in nonlinear anal-ysis[M]. Berlin:Springer-Verlag,2005:x+439.
[64]ZEIDLER E. Nonlinear functional analysis and its applications. Ⅰ[M/OL]. New York:Springer-Verlag,1986:xxi+897. http://dx.doi.org/10.1007/978-1-4612-4838-5.
[65]ZEIDLER E. Nonlinear functional analysis and its applications. Ⅱ/A[M/OL]. New York:Springer-Verlag,1990:xviii+467. http://dx.doi.org/10.1007/978-1-4612-0985-0.
[66]ZEIDLER E. Nonlinear functional analysis and its applications. Ⅱ/B[M/OL]. New York:Springer-Verlag,1990:i-1202. http://dx.doi.org/10.1007/978-1-4612-0985-0.
[67]ZEIDLER E. Nonlinear functional analysis and its applications. Ⅲ[M]. New York:Springer-Verlag,1985:xxii+662.
[68]ZEIDLER E. Nonlinear functional analysis and its applications. IV[M]. New York:Springer-Verlag,1988:xxiv+975.
[69]DEIMLING K. Nonlinear functional analysis[M]. Berlin:Springer-Verlag,1985: xiv+450.
[70]SMALE S. An infinite dimensional version of Sard's theorem[J]. American Jour-nal of Mathematics,1965,87(4):861-866.
[71]SCHROTER T. On a special class of nonlinear Fredholm integral equations of the first kind[J]. Computing,1997,58(3):259-279.
[72]钟承奎,范先令,陈文塬.非线性泛函分析引论[M].[S.1.]：兰州：兰州大学出版社,2004.
[73]HRYCAK T, ISAKOV V. Increased stability in the continuation of solutions to the Helmholtz equation[J/OL]. Inverse Problems,2004,20(3):697-712. http://dx.doi.org/10.1088/0266-5611/20/3/004.
[74]MAKINO M, OISHI S. Constructive analysis for infinite-dimensional nonlin-ear systems-Infinite-dimensional version of homotopy method[J]. Electronics and Communications in Japan (Part III:Fundamental Electronic Science),1991,74(8):1-10.
[75]DEIMLING K. Nonlinear functional analysis[M].[S.1.]:Courier Dover Publica-tions,2011.
[76]HAN B, FU H, LIU H. A homotopy method for well-log constraint waveform inversion[J]. Geophysics,2006,72(1):R1-R7.
[77]ZHAO J, LIU T, LIU S. An adaptive homotopy method for permeability estima-tion of a nonlinear diffusion equation[J]. Inverse Problems in Science and Engi-neering,2013,21(4):585-604.