【作者】 刘明辉；

【导师】 马富明；

【作者基本信息】 吉林大学 ， 计算数学， 2014， 博士

【摘要】 本文主要讨论了非线性不适定问题数值求解及在近场光学问题中的应用.第一章是绪论,介绍了散射问题和反散射问题的数学模型和常见方法,近场光学的模型,以及不适定问题的概念和正则化方法.第二章主要介绍谱截断矩量法的相关理论及其在处理严重不适定问题中的应用,特别是在近场光学问题中的应用.通过数值计算可以看到,相比常规的正则化方法,比如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<x1<L}上探测,Ω(?)R12为包含D和r的紧集.察物体由倏逝波照射,并在远场探测散射场.光子扫描隧道显微镜(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与正则化参数α.这样,我们得到近场光学问题的同伦算法.更多还原

【Abstract】 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 opticsNear-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 equationLet 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 obtainthatwhereand 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 haveAnd 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 applications2.1Moment method with truncated SVDMoment 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 isBy Riesz Theorem, there exists kj E X, such thatKz(sj)=(z,kj),(?)z∈X, j=1,…,n.Let Xn=span{kj,j=1,…,n}, then the moment solution is given bywhere a solves the linear system withThen the singular value decomposition(SVD)of A is made bywhere U,V are unitary withandwhere λ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 isDenotethen 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, andSuppose that a solves(6)with β the right-hand side,and assume there exists z∈Cn|z|≤E,such that a=A*z,thenFor constant c＞0,choose∧=cδ/E,then which is optimal,where2.2Application to linearized near-field optics problemThe linearized near-field optics problem,i.e. inverse problem concerning Born approximation,is to solve the scatterer q(x)fromNote 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 haveSubstituting the fundamental solution G(x,y)and ul into(11),and realizing the depen-dence of us on α,we rewrite(11)asMultiplying both sides of(12)by e-iωx1and integrating from-∞to∞with respect to x1, and noting that(13) can be written asRealizing that(14) can be reduced toWhen|ω|＜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 equationfor 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 wherethen 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 aswhere the coefficients aj’s are the solution of the systemwith 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 CnWe 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 writewhere λ1≥…≥λn≥0,let α(δ):=∧2≥λK2for some K(1≤K≤n), Σ=diag(λ1-1,λ2-1,…,λK-1,0,…,0),thenSo the moment solution with truncated SVD to(18)is given byBy 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 andAnd 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,thenwhereBy 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 applications3.1Homotopy with Tikhonov regularizationWe 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 problemwhere F:D(F)(?)X→Y is nonlinear ill-posed(compact)operator.Construct minimization funcitionalwhere 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 derivativeHomotopy 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 pathNumerically,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 methodStefenssen 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 skillsAs 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 asthen 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 MethodSimilar 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.LetIn 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 becomeswhich are ill-posed linear operator equations,so we can solve them by some kind of regularization method,such as Tikhonov regularization,SVD method.Euler methodBy 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 haveFor a fixed point x, the above equation is a linear operator equation with respect to dx. For0=λ0＜λ1＜…＜λn=1, we haveFor the homotopy (29), we haveFor a fixed point x, the above equation is linear operator equation with respect to dx. For0=λ0＜λ1＜…＜λn=1, we haveFrom λ=1-δ to λ=1, the linear operator equation becomeswhich 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 opticsIn 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 conditionwhere Σ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 v1satisfiesand radiation condition.(37)-(36)and we getUsing the relationship and omitting higher order terms, and denoting ω=S"(q)h1h, we havewhere v=S’(q)h1satisfieswhere 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 relationshipand S’(q)h1=v, we can get that, for given S(q)-y, solve Φ that satisfiesand radiation condition, thenThus, 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,…,letif‖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 asthen 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.更多还原