Crack Detecting via Newton’s Method

Crack detecting is an important issue in many areas of science and engineering. Often it is dangerous or it is only laborious to find the cracks directly, e.g., the detection of cracks within the nuclear power plants or of the cracks buried beneath the earth. One possible way out is to make use of the phenomenon of wave scattering. Sending an incident wave into the area in which we are interested, the information about the crack is hidden in the scattered wave that can be measured at some convenient place distance away. Thus, crack detection is equivalent to the extraction of the hidden information from the measured data. This is what we want to deal with here.

about the crack from the measured data. Later we will see that both for the direct and the inverse problems, the same equations are used almost all the time. This is beneficial both in the apprehension and in the computation.

Direct scattering problem
The scattering problem is mathematically modelled by an exterior boundary value problem governed by Helmholtz equation with some prescribed boundary conditions. The two dimensional problem arises when the obstacle is an infintely long cylinder. Thus, a planar crack can be seen as the intersction of a plane with an infintely long thin cylinder which runs in the direction normal to the plane. Mathematically, a planar crack can be given by a regular non-intersecting C 3 -smooth open arc Γ ⊂ R 2 which can be described as 1], z ∈ C 3 [−1, 1] and |z (s)| = 0, ∀s ∈ [−1, 1]}.
The two end points of the crack are denoted by z * −1 , z * 1 respectively. The left hand side and the right hand side of the crack are written by Γ + and Γ − respectively. The unit normal to Γ + is denoted by ν. Further we set Γ 0 := Γ \ {z * −1 , z * 1 }. The direct scattering problem that we are considering is as follows:

Problem 1. (DP)
Given an incident planar wave u i (x, d) := e ik<x,d> with a wave number k > 0 and a unit vector d giving the direction of propagation, find a solution u s ∈ C 2 (R 2 \ Γ) ∩ C(R 2 \ Γ 0 ) to the Helmholtz equation which satisfies the Neumann boundary conditions on both sides of the crack and the Sommerfeld radiation condition lim r→∞ √ r ∂u s ∂r − iku s = 0, r := |x| (3) uniformly for all directionsx := x |x| .
In (2), the limits are required to exist in the sense of locally uniform convergence.
Note that the boundary conditions (2) can be reformulated as homogeneous Neumann conditions for the total field u := u i + u s , i.e., ∂u ± ∂ν = 0. This means that the normal component of the velocity of the total wave vainishes on the crack, i.e., the crack is sound-hard. Using boundary integral equations, this direct problem can be solved via the layer approach, see for example [2] and [3]. Indeed, in terms of the fundamental solution to the Helmholtz equation in R 2 Φ(x, y) := i 4 H (1) 0 (k|x − y|), x = y we have the following theorem which ensures the unique solvability of the direct scattering problem.
Theorem 1. The direct Neumann problem 1 has a unique solution given by where ϕ ∈ C 1,α, * (Γ) is the (unique) solution to the following integral equation and for 0 < α < 1, the function space C 1,α, * (Γ) is defined by Applying the Green's integral theorem, the uniqueness is ensured by the Rellich's lemma and the radiation condition. Using the potential ansatz, the solvability of the boundary value problem is then converted to the solvability of the induced boundary integral equation (5) which can be determined by the Riesz theory. For details we refer to [4].
At this place we want to point out that in the scattering problem, one is particularly interested in the far-field pattern u ∞ of the scattered field u s . The far-field pattern describes the behavior of the scattered wave at infinity uniformly for all directionsx ∈ Ω := {x ∈ R 2 ||x| = 1}. The one-to-one correspondence between radiating waves and their far field patterns is established by the Rellich's lemma. In the case of a sound-hard crack, the far-field pattern is declared via with the density function ϕ given by theorem 1 and the constant ρ = k 8π e −i ß 4 .
For further treatment, we would like to transform our integral equation (5) into operator form. For this purpose, the following integral operators may be defined: To reduce the hypersingularity of the operator T 0 , we use the Maue's identity to split it into two milder parts (for a proof of this splitting, see for example Theorem 7.29 in [3]) where ϑ is the unit tangent vector. The integral equation (5) now becomes where we set g = − ∂u i ∂ν .

Inverse scattering problem
After introducing the direct scattering problem in the last section, we consider the following inverse problem:

Problem 2. (IP)
Determine the crack Γ if the far-field pattern u ∞ is known for one incident wave.
About this inverse problem, the first thing to ask is the uniqueness, that is, the identifiability of the crack. Unfortunately, there exists no theoretical result for our setting of inverse problem. However, if all the far-field patterns from all possible incident directions are measured, the crack can be uniquely identified. We quote the following result from [5].
Theorem 2. Assume that Γ 1 and Γ 2 are two sound-hard cracks with the property that for a fixed wave number k > 0, the far-field patterns u 1,∞ and u 2,∞ coincide for all incident directions d. Then we have Our aim is to find and to reconstruct the crack from the measured far-field data resulted from just one single incident wave. Despite of the lack of a theoretical proof, our setting of inverse problem arises naturally from the viewpoint of the real applications. It is always desirable to find the crack with less effort.
For the detection of the crack, the only available information is the measured far-field data. Equation (6) relates the far-field pattern u ∞ with the crack Γ and is therefore suitable as the starting point for the reconstruction. For future reference, we rewrite it in the more compact form and call it the far-field equation as in the literature. The inverse problem is then equivalent to the solving of this equation. However, the sloving of (8) is not as straight forward as it appears. Mathematically this equation is not even solvable because of the nature of the operator F. Being a compact operator in infinite dimensional function space, F cannot have a bounded inverse. This means that (8) cannot be solved in any reasonable normed space.
Therefore, certain regularization techniques must be introduced at this place. For the sake of completeness, we briefly discuss regularization for linear compact operators in the next section. We note that the operator F will be differently interpreted for different methods in the following sectons. But the right hand side of (8) will remain the same.

Regularization
Definition 3 (Regularization). Assume that X, Y are normed spaces. Let the operator A : X → Y be linear, bounded and injective. A family of bounded linear operators R α : Y → X, α > 0 is called a regularization scheme for Aϕ = f (9) if it satisfies the following pointwise convergence In this case, the parameter α is called the regularization parameter.
From the above definition, it is easy to see that if the operator A has a bounded inverse, then we can simply choose R α = A −1 , ∀α > 0. In the case where the inverse of the injective operator A is not bounded, for example when A is compact and the dimension of X is infinite, then difficulty arises since we are trying to approximate an unbounded operator with a sequence of bounded operators. The task of regularization is to find a stable approximate solution. The following theorem, which is proved in [3] (theorem 15.6), shows the limitations of the regularization.
Theorem 3. Let X and Y be normed spaces with dimX = ∞ and let A : X → Y be linear and compact. Then for a regularization scheme the operators R α cannot be uniformly bounded with respect to α, and the operators R α A cannot be norm convergent as α → 0.
Because of the lack of uniform convergence, the choice of the regularization parameter α is apparently crucial. To have a vivid impression of this, let's examine the aproximation error in a more pratical setting where measurement errors or noises are present. Let's denote the contaminated data by f δ assuming an error level δ, i.e., f − f δ ≤ δ. For α > 0, we wirte the regularized approximation as ϕ δ α := R α f δ The total approximation error of the regularization scheme can be written as From this, we have the following error estimate The right-hand side of (10) reveals the difficulty of the approximation. As α tends to 0, the second term of the right-hand side becomes small because of the regularization. At the same time, the error δ will be amplified by the factor R α which tends to infinity. Therefore, the total error is dominated by the first term and is much larger than the datat error. For larger α, the first term on the right-hand side is smaller but the second term will be larger. Hence certain strategy for choosing α must be developed to balance the effects of the two factors. The discrepancy principle proposed by Morozov ([6], [7]) which based on the consideration that a problem cannot be solved more accurate than the data error, is a natural a posteriori strategy for determining the parameter α. More precisely, one looks for an α > 0 such that for some prescribed γ ≥ 1. It can be shown that there always exists a smallest α > 0 for which the approxiamte solution f δ satisfies the above condition. (see [3]) Although the existence for a "best" α is theoretically assured, there is no obvious way to select it in practice. In most cases, the only way to choose a good parameter is by trial and error. For our purpose, it is sufficient to consider regularization in the Hilbert space setting. We denote by A * the adjoint operator of A. Without going into details, we adapt the following regularization scheme.
Applying the Tikhonov Regularization scheme to find an approximate solution to (9), we actually solve the follwing equation for ϕ α .

Classical Newton's method
As noted earlier, (8) connects the far-field pattern with the crack and is therefore suitable for the inverse task of finding the crack. It is therefore natural to solve it for the unknown crack when one has the measured far-field data at hand. For this nonlinear problem, Newton's method provides a simple way to handle it. Utilizing the first Frechét derivative which existence is proved in [5], Newton's method approximates the nonlinear equation (8) with a linear one. This means that instead of solving (8), we endeavor to solve the linear equation for an update q. We remark here that the derivative F is still a compact operator, this means that the above equation can't be solved directly. However, we can apply the regularization technique from the last section. Thus, the following equation has to be solved Numerically, starting from an initial guess γ of the unknown crack, (13) finds an update q. The same equation is then solved iteratively with γ replaced by γ + q until some convergence criterion is met.
This scheme is conceptually very simple, yet the computation of the Fréchet derivative F is not straightforward. To calculate F , one must first solve another direct problem since F is characterized as the solution of that problem (see [5]). Thus additional computational cost arises.
At this place, we want to discuss the uniqueness of the equation (12). As pointed out in [8], the Fréchet derivative F is not injective. The null space of this operator is given by This reflects the fact that different parametrizations of the arc leading to the same set of points which in turn give the same far-field pattern. We can avoid this ambiguity by limiting our solution space to the set of arcs representable as the graph of a function as suggested in [8].
Once this restriction is made, the operator F is then an injective linear compact operator. The same restriction will aslo be made in the next two sections.

Nonlinear integral equations method
Based on the reciprocity gap principle, Kress and Rundell [9] proposed a method using nonlinear integral equations for solving an inverse conducting problem avoiding the computation of a forward problem at each iteration step. The idea behind this approach is quite simple: in a doubly connected bounded domain to which the second Green's theorem applies, if one knows the Cauchy data of a harmonic function on the outer boundary and the type of the boundary condition on the interior boundary, then one can determine the Cauchy data on the interior boundary and also the interior boundary itself. This method is extended to other boundary conditions and cracks in [10] and to obstacle scattering in [11]. Based on a different system of two nonlinear integral equations, a method with far-field data is proposed in [12] without the need of an auxiliary curve around the obstacle.
To understand this method, let's first reformulate our equations for the direct problem. In terms of Γ and ϕ, we can define the operators B : and respectively. Now we can consider the following system of operator equations Note that the operator B is just the boundary operator which maps the unknowns to the Neumann boundary data. The operator F calculates the far-field pattern. From the viewpoint of the direct problem, this system is just the same as the equations system (5) and (6) if Γ is given in advance.
To solve the corresponding inverse problem, the idea of [12] is to use this same system (16). The reasoning is very simple. If Γ solves the inverse scattering problem, then it follows directly from the solution theory of the direct problem that system (16) is satisfied. Conversely, if the pair (Γ, ϕ) solves (16), the first equation of (16) ensures that the total field u defined in theorem 1 satisfies the homogeneous Neumann boundary conditions on Γ. The second equation in (16) then ascertains the correct far field pattern for the scattered field u s . From the uniqueness theorem it follows that Γ is the solution of the inverse problem. We have thus the following main theorem (see [12]).

Theorem 5. Γ is the solution of the inverse problem if and only if Γ, ϕ solve the system of nonlinear integral equations (16).
We note that, mathematically, the system (16) illustrates the spirit of "inverse problems". Indeed, giving the type of the boundary conditions, (16) is just our system of equations for the direct problem which solves the far field pattern of the scattered field if one knows the crack. On the other hand, this very same system is used to find the crack from the knowledge of the far field pattern which is just the solution of the corresponding direct problem. This aspect also demonstrate one of the important features of our method. No new equations are created. What we use to solve the inverse problem is just what we have from the direct problem.
According to theorem 5, our task now is to solve the system (16). After parametrising the above system with the cosine substitution t = cos τ (see [13]) incorporated to take care of the square root singularities of the solution u s at the crack tips, the system (16) can be written in the following form and γ(τ) := (z • cos)(τ), ψ(τ) = sign(π − τ)((ϕ • cos)(τ)), n = (ν • γ) · |γ |. Note that we make an odd extension of our system (16) to the closed interval [0, 2π] which is not necessary but helpful in the theoretical treatment [4].
It is noted that all the operators are nonlinear with respect to γ. Since our numerical method will be based on the Newton's method, we need also the Fréchet derivatives of the operators w.r.t. the boundary γ. The Fréchet derivatives of the integral operators A s are simply given by the Fréchet derivatives of their kernels (see [14]). For brevity, we set a ⊥ = (a 2 , −a 1 ) t if a = (a 1 , a 2 ) t . Also we write Δγ = γ(τ) − γ(σ) and Δq = q(τ) − q(σ) . Using the differentiation rules of the Hankel functions H 0 (z) = −H 1 (z) and (zH 1 (z)) = zH 0 (z), we have Therefore, to solve the parametrized system (17), Newton's method suggests the solving of the following system of linear approximations Numerically this system is solved iteratively: with a current approximation (ψ, γ), the system (18) has to be solved for the pair (χ, q). The updated data are then given by ψ + χ for the density function of the integral operators and γ + q for the unknown crack. For brevity, we rewrite the system (18) in the operator form: Since this linear equation of the first kind is still ill-posed, we have to incorporate some regularization scheme. Instead of solving (19) directly, we solve the following regularized equation with two to be determined regularization parameters α and β.
The major advantage of this method as compared to the classical Newton's method is that the Fréchet derivatives of the integral operators can be easily computed. They are obtained simply by solving the above system which is nothing more than the Gaussian elimination.

Two-steps method
The idea of the nonlinear integral equations method is based on the observation that the inverse problem can be formulated in an equivalent two-by-two system (16) of nonlinear equations. This system is then treated as a coupled system of two variables (γ, ψ). In the regularized version (20), two parameters (α, β) are therefore needed at the same time. Since these parameters are in general selected by trial and error, it is rather complicated. However, the system (16) can be solved in another way. Inspired from the solution theory of the direct problem, This system can be seen as the whole procedure for the calculation of the far-field pattern of the scattered field from a given crack resulted from an incident planar wave in two steps. Indeed, the first equation of (16) solves the scattered field from the crack and the second one calculates the far-field pattern from the solution of the first equation. This motivates our inverse scheme: the inverse solver should be inverse to the direct solver. There are several possibilities to solve (16) in two steps separately. One can solve the first equation first and the second one next, or the other way round. One has further choice concerning which quantity is to be solved at each individual step. Since our method is based on Newton's iteration method, an initial guess for γ or for ψ or for both is unavoidable. It is therefore natural to start with a initial guess for the crack and solve the first equation for ψ. This is beneficial since this is just the direct problem defined by problem 1. Even better is, theorem 1 ensures the unique solvability of the equation for any crack. Next we can solve the second equation for an update of the crack. Note that this step is nonlinear and ill-posed. With slightly abused notations, we have thus splitted the system (16) into the following two separate steps: Note that this system is formally the same as system (17) except that the first variable is fixed in each single equation. That is, γ is fixed in (21a) while ψ is fixed in (21b). In the actual numerical computation, we have to solve a regularized version of (21b). We therefore solve instead of (21b).
The algorithm for the two-steps method can be summarized as follows: 1. Given an initial guess γ 0 for the unknown crack.

Stopping criterion: Discrepancy principle
At this place we comment that from the numerical point of view this scheme is very attractive. As a variant of the nonlinear integral equations method, the Fréchet derivative of the integral operator is computed directly by solving the integral equation (23). Besides, the derivative is very simple as compared to those in the last section, since only the far-field operator which has a smooth kernel is to be differentiated. Numerically it can be easily solved via the rectangular rule, for example. Another advantage of this method is that it splits the problem into two smaller parts and thus makes the computation cheaper.

Numerical results
In this section we will demonstrate the applicability of the two-steps method via some examples. One can compare the results with those from the other methods in [5] and [12]. Using just one incident wave, we reconstruct the unknown crack from the measured far-field pattern at 16 points which are evenly distributed on the unit circle. For the direct problem, the forward solver is applied with 63 collocation points. To avoid committing an inverse crime, the number of collocation points used in the inverse problem is chosen to be different from that of the forward solver. We choose 31 collocation points in the inverse algorithm. In all our examples, we take the incident plane wave coming from the direction d = 1 For the reconstruction, we take the k-th Chebyshev monomials T k (x) = cos(k cos −1 x), k = 1, . . . , m as our basis functions. The selection of the Chebyshev polynomials is based on the fact that they can take care of the square root singularity in their own right. The updates for the crack in each iterative step can be written in the form with to be determined coefficients b 1 k , b 2 k , k = 0, . . . , m − 1. We test our scheme for two different wave numbers, k = 1, 3. The dimension of the test space is taken to be 5, i.e., m = 5. The starting curve, that is, the initial guess for the regularized Newton's method, is taken to be the straight line y = 0 in all examples. According to the discrepancy principle, the stopping criterion for the iterative scheme is given by the relative error which is taken to be 0.001 in case of exact data and 0.03 in the case where 3% data error are present. In all our figures below, the dotted line (blue) represents the initial guess. We denote by the dashed line (red) the true solution and by the solid line (black) the reconstruction.

Example 1.
For the first example, we take the arc which is a polynomial. The numerical results are given in the figures 1, 2 for k = 1 and k = 3, respectively. From the numerical results, we see that the reconstructions for exact data are very good. In the case where random errors are present, the reconstructions are not bad at all.

Example 2.
To demonstrate the benefit of our method, we choose a curve which does not belong to our solution space. To this end, we take the arc  We see from the figures 3, 4 that the reconstructions are very good, even in the case of erroneous data.  We see that in the case k = 1, the reconstructions are very good both with exact and erroneous data (Fig 5). In the case of more oscillations (k = 3, Fig 6), the result is quite good with respect to the bad initial guess for the iteration. We also note that although in this case, the Fréchet derivative of the far field operator is not injective, the actual reconstruction process is still running without any problem.

Conclusion
In this article, we try to detect a planar crack via iterative Newton's method with the knowledge of the measured far-field pattern. By the one to one correspondence of the scattered field and its far-field, the measurements can be taken anywhere outside the crack. The choice of boundary integral equations method reduces the dimension of the problem by one which makes the computation cheaper. The conceptual simplicity and its numerical accuracy make the Newton's method attractive. Besides, the two variants of the Newton's method (nonlinear integral equations method and two-steps method) on one hand simplifiy the calculation of the Fréchet derivatives and on the other hand help clearifying the idea of the inverse problem.

Author details
Kuo-Ming Lee Department of Mathematics, National Cheng Kung University, Taiwan