Numerical Method For Nonlinear Ill-posed Problems With Applications In Near-field Optics

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.