Identification of far field electric field based on near field distributed measurements.
In this note, we consider the evaluation of an incoming directed energy at parts of the domain based on distributed measurements collected at other parts of the domain. The eventual goal is to be able to detect incoming laser beams. The mathematical model involves the propagation of the electromagnetic waves through the atmosphere. Writing the Maxwell's equation in terms of the electric field, E(x), leads to the following vector wave equation 
[[nabla].sup.2]E + [nabla] [E.[nabla](log([epsilon]))] = [mu] [[partial derivative].sup.2] ([epsilon]E) / [partial derivative][t.sup.2],
where, [mu] is the perz meability and [epsilon](x) is the dielectric constant. Assuming a periodic incident wave, E(t, x) = U(x)cos([omega]t), one arrives at
[[nabla].sup.2] U + [[omega].sup.2] [mu][epsilon]U = - [nabla] [U.[nabla](log([epsilon]))]
for the amplitude of the electric field, U(x). It can be argued that, for a wide range of applications, the right hand side in the above equation can be neglected in comparison to the other two terms  (page 232). Assuming that the dielectric constant is uniform, i.e., [mu][epsilon](x) = 1, it is then clear that one is led to deal with the Helmholtz equation.
Under the above conditions, the problem of identifying an incoming laser beam is equivalent to identifying the boundary conditions at the part of the domain (outer) that is not accessible. At the part of the domain that is accessible, it is possible to obtain distributed measurements. This situation calls for the solution of an elliptic system, i.e., Helmholtz equation, for which a complete set of boundary conditions are not available. This problem is often refered to as Cauchy problem and it is well known that these problems are ill-posed [2, 3]. A similar situation can also arise in a number of other physical applications including corrosion detection , soil venting , bioelectric , and eddy-current modeling .
A number of investigators have considered this topic and developed various methods including Quasi-reversibility , Tikhonov regularizations [9, 10], quasi-boundary-value method , Carleman estimate , mollification regularization , and Bakus-Gilbert method .
The purpose of this note is to develop two algorithms for the solution of a Cauchy problem for Helmholtz equation. The first algorithm is the application of a recently developed method to the problem at hand. This method is quite versatile and can be applied to various systems. It was originally developed for parabolic systems [15, 16]. The second algorithm is specific to the mathematical situation at hand. It is more effective in recovering the unknown boundary condition. Both algorithms are iterative in nature and share the same structure. In what follows, we first describe the two algorithms in the next three sections. In section 5, we use a number of numerical examples to investigate the applicability of the algorithms.
2. Mathematical Formulation
Let S [subset] [R.sup.2] be a closed bounded set. Consider a 2-D Helmholtz equation given by
[DELTA]u + [k.sup.2]u = 0, x [member of] S [subset] [R.sup.2]. (1)
It is often the case that the condition on the part of the boundary, [partial derivative]][S.sub.1], is not known. In addition, the remaining part of the boundary [partial derivative][S.sub.2],with [partial derivative][S.sub.1][union][partial derivative][S.sub.2] = [partial derivative]S and [partial derivative][S.sub.1][intersection][partial derivative][S.sub.2] = 0, is accessible. It is possible to obtain both Dirichlet and Neuman condition on the part of the boundary that is accessible, i.e.,
u(x) = [tau](x), [partial derivative]u / [partial derivative]n (x) = [sigma] (x), x [member of] [partial derivative] [S.sub.2] [subset] [partial derivative] S. (2)
The problem of interest is now to solve for u(x) [member of] S based on the information on parts of the boundary only.
The algorithms presented in this note are iterative in nature. They assume a value for the missing boundary condition. At every iteration, they obtain corrections to the assumed value. This is the main part of these algorithms. It is possible to develop a number of different methods to accomplish this task which are the main contributions of this work. For clarity, both algorithms consists of the following three steps:
1: Assume a value for the unknown boundary condition and obtain the error field.
2: Use the additional boundary condition for the error field and obtain the correction to the assumed value in step 1.
3: Update the assumed value and go to step 1.
3. Algorithm 1
For simplicity, consider a square domain shown in figure 1.
The first algorithm consists of the following steps.
I. Assume an initial guess: In the above figure, the top portion of the domain cannot be reached and the boundary condition is unknown. It is often the case that, using the known boundary conditions, an initial guess for the unknown part of the boundary can be obtained. Let [??] (x) denote the initial guess for the missing boundary condition. The initial guess together with the given Dirichlet boundary conditions lead to a background field, [??](x), which satisfies the forward problem given by
[DELTA][??] + [k.sup.2][??]]u = 0, x [member of] S [subset] [R.sup.2]. (3)
This field satisfies the Dirichlet conditions at all four sides. Since the initial guess is almost always not the actual function, the computed field is not the actual field. Let e(x) denote the error, i.e., e = u - [??], and subtract the two equations (Eqn. (1) and Eqn. (3)) to obtain an equation for the error given by
[DELTA]e + [k.sup.2]e = 0. (4)
[FIGURE 1 OMITTED]
The error function e(x), shown in figure 2, is required to satisfy two boundary conditions on the part of the boundary at which the data is collected.
[FIGURE 2 OMITTED]
II. Obtain the error field e(x):
This is the main step in the algorithm. The above problem for the error field is also ill-posed. In the following section we develop a number of methods to compute the error field.
III. Update the assumed value [??]:
The assumed value for the missing boundary condition can be updated according to
[??] (x) = [??] (x) + e(x, 1) (5)
In the rest of this section, we develop a number of methods to accomplish step II of the above algorithm.
3.1. Differential formulation
The problem that one needs to solve for the error field e(x, y), shown in figure 2, is still ill-posed. It is possible to consider two well-posed problems for the error field. This approach was originally developed for ill-posed parabolic systems [17, 16]. The two well-posed elliptic problems for e(x) are given by
[FIGURE 3 OMITTED]
Note that g(x) is the collected data, and [[??].sub.y](x) can be computed form the background field.
The same symbol, e(x), is being used here to denote the error in both of these problems. In what follows, it will specifically be clear which problem is being considered. Both of these problems are well-posed since we are treating h(x) as a known quantity. Using a second-order accurate central differencing for the Helmholtz operator, these two equations lead to
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], Problem I (6)
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], Problem II (7)
where now, the matrices [A.sub.0] and [A.sub.1] are known. We are using a bar, i.e., [bar.*], to denote vector representations of variables. If the domain is divided into [n.sub.e] equal intervals in each direction, i.e., [DELTA]x = [DELTA]y = 1 / [n.sub.e], then the size of the unknown vector [bar.h] is n = [n.sub.e] + 1, and [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] with [n.sup.2] = n x n. The condition [e.sub.y] = g(x) - [[??].sub.y] shows up in the vector [bar.f] in Eqn. (7). Also, the boundary condition at y = 1, which is represented by [bar.h], is not known. This quantity appears in the vectors [bar.p] and [bar.f]. In essence, we have two equations and two unknowns, namely, e and h. Eliminating the unknown vector [bar.e] leads to
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (8)
Note that both problems in figure 3 are well-posed and the matrix [A.sub.0] can be readily inverted. It is now possible to obtain an equation for the unknown boundary condition, i.e. [bar.h]. The coefficient matrix in Eqn. (8) can be partitioned according to
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].
This is due to the fact that the first ([n.sup.2] - n) entries in the column vector on the left hand side of Eqn. (8) are zeros. After denoting the last n columns of the coefficient matrix by the ([n.sup.2] x n) nonsquare matrix [PHI], then the above equation leads to
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (9)
The above equation has the unknown vector [bar.h] in the right hand side as well. Moving the unknown vector [bar.h] in the last n equation from the left hand side to the right hand side leads to
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (10)
where, the matrix [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is equal to the matrix [PHI] after subtracting an (n x n) identity matrix to account for the above rearrangement. The above system can now be solved for the unknown vector [bar.h]. It is noted that the coefficient matrix is not square. In addition, as is often the case with ill-posed problems, the coefficient matrix is highly singular. Tikhonov regularization can now be used to stabilize the inversion according to
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (11)
where, the matrix [THETA] [member of] [R.sup.(nxn)] is a full rank, first-order differential operator, and [beta] is a constant coefficient. Note that we are also defining new variables for simplicity. The least-square solution to the above over-determined system is now given by
[bar.h] = [[[[GAMMA].sup.T][GAMMA]].sup,-1] [[GAMMA].sup.T] [bar.[sigma]. (12)
This completes the third step of the algorithm. It is now possible to update the assumed value according to Eqn. (5). We next proceed to develop a second formulation.
4. Algorithm 2
The second algorithm is similar to the first algorithm, in the sense that, it is also an iterative algorithm. However, it computes the correction to the assumed value, i.e., step III, differently.
I: Assume an initial guess.
II: Obtain the error field e(x).
III: Use algorithm 1 to obtain a value for the missing boundary condition, i.e., h(x) = e(x, 1).
IV: Use a semi-analytical method to compute a better estimate of the missing boundary condition.
V: Use the value just computed in step IV as the regularization matrix and update the computed value for the correction to the unknown value. Repeat steps I-V until convergence.
The first two steps in this algorithm are similar to the first two steps in algorithm 1. The formulation for the error field makes it possible to improve the estimate for the missing boundary condition. Consider the error field given in figure 2. The procedure that is to follow is very similar to the variation of parameter's formulation for the solution of nonhomogeneous transient boundary value problems  (page 103). If one assumes a solution of the form
e(x, y) = [[infinity].summation over (j = 1)] [B.sub.j](y)sin(j[pi]x) (13)
for the error, then using orthogonality condition, it is possible to obtain relations for the coefficients according to
[B.sub.j](y) = 2 [[integral].sup.1.sub.0] e(x, y)sin(j[pi]x)dx,for, j = 1,2, ... [infinity]. (14)
Differentiating the above equations with respect to y, twice leads to
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (15)
Rearranging the above equations leads to
[d.sup.2][B.sub.j(y) / [dy.sub.2] + [k.sup.2][B.sub.j] = -2 [[integral.sup.1.sub.0] [[e.sub.xx]sin[j[pi]]x)dx, j = 1,2, ... [infinity]. (16)
Integrating the right hand side by parts twice leads to
[d.sup.2][B.sub.j(y) / [dy.sub.2] + [k.sup.2][B.sub.j] = -2[(j[pi]).sup.2] [[integral.sup.1.sub.0] [[e.sub.xx]sin[j[pi]]x)dx, j = 1,2, ... [infinity]. (17)
It is possible to treat the first few terms in the series differently. If k> jn then,
[d.sup.2][B.sub.j(y) / [dy.sub.2] + ([k.sup.2] - [(j[pi]).sup.2])[B.sub.j](y) = 0 if, k > j[pi] (18)
[d.sup.2][B.sub.j(y) / [dy.sub.2] + [k.sup.2] [B.sub.j](y) = -2 [(j[pi]).sup.2], [[integral.sup.1.sub.0] e(x,y) sin (j[pi]x)dx = [[mu].sub.j] (y), if, k < j[pi], (19)
where, for simplicity, we are also defining a new variable, [[mu].sub.j](y). It is then possible to obtain exact solutions for the terms for which k > j[pi]. For the cases where k < j[pi], the exact solution leads to exponential functions which will lead to instability. Therefore, for the remaining terms, it is possible to set up an iteration. It is possible to use the first algorithm, presented in section 3, to obtain an initial error field e(x, y). Then, using this error field, the integral terms on the right hand side can be computed. It is then possible to obtain solutions for the above differential equations. The required boundary conditions are given by
[B.sub.j](0) = 0, [dB.sub.j](0) / dy = [[integral].sup.1.sub.0] [e.sub.y](x,0)sin(j[pi]x)dx = [v.sub.j] j = 1,2, ... [infinity], (20)
where, for simplicity, we are defining new scalars [v.sub.j]. The exact solutions are given by
[B.sub.j] (y) = [v.sub.j] / [square root of [k.sup.2] - [(j[pi]).sup.2]] sin ([square root of [k.sup.2] - [(j[pi]).sup.2]y), k > j[pi](21)
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (22)
Then the error at the top can be updated according to
h(x) = e (x,1) = [[infinity].summation over (j = 1)] [B.sub.j] (1) sin (j[pi]x). (23)
For clarity, we state the second algorithm again.
1: Assume a value for the unknown boundary condition.
2: Obtain the background field and, in turn, the bounday conditions for the error field, i.e. figure 2.
3: Use the first algorithm to obtain an initial value for the error field.
4: Use Eqn. (13) and Eqns. (20-23) to obtain a better estimate of the error field , e(x, y) and the unknown boundary condition h(x) in Eqn. (23).
5: Use this value of the unknown boundary condition as the regularization term in Eqn. (11) and obtain a new value for the unknown boundary condition for the error, i.e., h(x) = e(x,1)
6: Update the assumed value for the missing boundary condition and go to step 2.
Step 5 in the present algorithm uses this value of h(x) as the regularization matrix. This process simply Eqn. (11) with
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (24)
where, [??] is the vector representation of the function given in Eqn. (23). Note that, since I is an n x n identity matrix, with nonzero scalar [[beta].sub.2], the coefficient matrix has a full rank. The above linear system represents an over-determined system, and least-square method can be used to obtain the solution.
5. Numerical implementation and examples
In this section, we use a number of examples to investigate the applicability of the proposed algorithms. The above algorithms involves the numerical solution of the Helmholtz equation. It is well known that, for high frequencies, the numerical solution of Helmholtz equation remains challenging [19,20]. For simplicity, consider a 2-D domain and assume that the known boundary conditions are given by
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (25)
[FIGURE 4 OMITTED]
It is also possible to collect data at the lower boundary in the form of [u.sub.y](x, 0). The condition at the top boundary is unknown.
Example 1: (Algorithm 1)
First Consider recovering an unknown boundary condition given by
u (x,1) = 1 + exp [- [(x - 0.58).sup.4] / 0.00003]. (26)
If one divides the domain into [n.sub.e] = 70 equal intervals in both x and y directions, the forward problem is composed of a linear system with dimension [n.sup.2] = [71.sup.2]. This mesh is sufficient to produce accurate numerical results for the examples discussed in this note with k = 10. It is often the case that the collected data is noisy. Therefore, after generating the data numerically, the data is contaminated with a zero mean Gaussian white noise. Using a second-order accurate finite difference method, it is now possible to obtain various matrices in the above algorithms.
Consider the first algorithm. The first algorithm calls for the least-square solution of Eqns. (11-12). Consider the above least-square problem before the inclusion of the Tikhonov regularizations. It is instructive to look at the singular values of the coefficient matrix [GAMMA] (for [beta] = 0), or the eigenvalues of the symmetric non-negative matrix [[[GAMMA].sup.T][GAMMA]]. Figure 4 shows the eigenvalues of the symmetric matrix [[[GAMMA].sup.T][GAMMA]] for k = 17. The eigenvalues are normalized with respect to their highest value. The figure shows that the coefficient matrix has only a number of non-zero eigenvalues and Tikhonov regularization is necessary for stable inversion.
We next proceed to recover the unknown function. Figure 5 shows the recovered function, using the first algorithm. The noise to signal ratio is %4. The regularization parameter is [beta] = 45 and the algorithm is able to reduce the error by four orders of magnitude Error : [10.sup.4 [right arrow] 1. After four iterations, there are no significant reduction in the error. The error is the difference between the measured data and the computed value according to
Error = [n.summation over (l = 1) [[g.sub.l] - [[[??].sub.l].sup.2] = [n.summation over (l = 1) [[[gamma].sub.l].sup.2]. (27)
[FIGURE 5 OMITTED]
Example 2: (Algorithm 1)
We next consider recovering an unknown boundary condition with two targets. Figure 6 compares the recovered function and the actual function. The noise to signal ration is also %4. The regularization parameter is [beta] = 45. The algorithm can recover the unknown function with two targets. Similar to the first example, It can reduce the error by four orders of magnitude Error : [10.sup.4] [right arrow] 1. After four iterations, there are no significant reduction in the error, or any improvement in the recovered function.
Example 3: (Algorithm 2)
We next proceed to apply the second algorithm to the Cauchy problem at hand. For this case, we use the first algorithm to obtain an initial value for the error field in the third step of the algorithm. Figure 7 shows the recovered function when the second algorithm is used. The same amount of noise is added to the data. For this case, the Tikhonov regularization parameter for the differential formulation is set at [beta] = 45. The solution obtained at step III, is also used as a regularization matrix with [[beta].sub.2] = 0.1. The amount of error is reduced by four orders of magnitudes (Error: [10.sup.4] [right arrow] 1). After two iterations, There are no signifant improvement in the recovered function.
Example 4: (Algorithm 2)
We next proceed to recover an unknown function with three targets. Figure 8 compares the actual function with the recovered unknown boundary condition. The Tikhonov regularization parameter for the differential formulation is set at [beta] = 45. The solution obtained at step III, is also used as a regularization matrix with [[beta].sub.2] = 0.1. The amount of error is reduced by four orders of magnitudes (Error: [10.sup.4] [right arrow] 1). Similar to example 3, after two iterations, there are no significant improvement in the recovered function.
[FIGURE 6 OMITTED]
[FIGURE 7 OMITTED]
[FIGURE 8 OMITTED]
In this note, we studied the application of two new algorithms to a Cauchy problem for the Helmholtz equation. The algorithms seek to indentify the missing boundary conditions at the part of the domain that is not accessible. The model under study is suitable to indentify the amplitude of an incoming directed energy. Future studies will seek to identify more specific information on the incoming directed energy.
Both methods show good robustness to the noise. The second algorithm was able to recover a better estimate of the unknown boundary condition. The results presented here are for low range frequencies. This is entirely a numerical issue. Numerical results can also be greatly improved by using a finer mesh and/or, more accurate finite-difference approximations.
 J.W. Strohbehn (Ed.). Laser Beam Propagation in the Atmosphere. Springer-Verlag, New York, 1978.
 G. Alessandrini, L. Rondi, E. Rosset, and S. Vessella. The stability for the cauchy problem for elliptic equations. Inverse Problems, 25:123004(47pp), 2009.
 L.E. Payne. Improperly Posed Problems in Partial Differential Equations. SIAM, Philadelphia, PA, (1975).
 G. Inglese. An inverse problem in corrosion detection. Inverse Problems, 13(4):977-994, 1997.
 M. Slodicka and H. De. Schepper. On an inverse problem of pressure recovery arising from soil venting facilities. Applied Mathematics and Computation, 129 (2-3):469-480, 2002.
 C.R. Johnson. Computational and numerical methods for bioelectric field problems. Critical reviews in Biomedical Engineering, 25:1-81, 1997.
 Marian Slodicka and Valdemar Melcher. An iterative algorithm for a Cauchy problem in eddy-current modelling. Applied Mathematics and Computation, 217:237-246, 2010.
 R. Lattes and J. L. Lions. The method of Quasi-Reversibility. American Elsevier, New York, (1969).
 Houde Han, Leevan Ling, and Tomoya Takeuchi. An energy regularization for Cauchy problems of Laplace equation in annulus domain. Communications in computational physics, 9(4):878-896, 2011.
 X.T. Xiong. A regularization method for a Cauchy problem of the Helmholtz equation. Journal of computational and Applied Mathematics, 233:1723-1732, 2010.
 Xiao-Li Feng, Lars Elden, and Chu-Li Fu. A quasi-boundary-value method for the cauchy problem for elliptic equations with nonhomogeneous neumann data. Journal of inverse and ill-posed problems, 18:617-645, 2010.
 Hui Cao, Michael V. Klibanov, and Sergei V Pereverzev. A Carleman estimate and the balancing principle in the Quasi-reversibility method for solving the cauchy problem for the laplace equation. Inverse Problems, 25(3):035005(21-pp), 2009.
 Hao Cheng, Xiao-Li Feng, and Chu-Li Fu. A mollification regularization method for the cauchy problem of an elliptic equation in a multi-dimensional case. Inverse Problems in Science and Engineering, 18(7):971-982, 2010.
 Y.C. Hon and T. Wei. Backus-gilbert algorithm for the cauchy problem of the laplace equation. Inverse Problems, 17(2):261-271, 2001.
 M. Tadi. Evaluation of diffusion coefficient based on multiple forward solutions. Applied Mathematics and Computations, 216(12):3707-17, 2010.
 M. Tadi. A computational method for an inverse problem in a parabolic systems. Discrete and Continuous Dynamical System-B, 12(1):205-218, 2009.
 M. Tadi. An iterative method for the solution of ill-posed parabolic systems. Applied Mathematics and Computations, 201(1-2):843-851, 2008.
 Glen E. Myers. Analytical Methods in Conduction Heat Transfer. Genium, NY, (1987).
 M. Navabi, M.H. Kamran Siddiqui, and J. Daragahi. A new 9-point sixth-order accurate compact finite difference method for the Helmholtz equation. Journal of Sound and Vibration, 307:972-982, 2007.
 G. Sutmann. Compact finite difference schemes of sixth order for the Helmholz equation. Journal of Computational and Applied Mahematics, 203:15-31, 2007.
M. Tadi1 (1)
Department of Mechanical Engineering, University of Colorado Denver, Campus Box 112, P.O.Box 173364, Denver, CO 80217- 3364
Naval Post-Graduate School, Monterey, CA 93943
(1) This work was performed while the author was a National Research Council Senior Research Fellow at the Naval Postgraduate School.
(2) This work was supported in parts by the Counter-Directed Energy Weapon program (C-DEW) of the Office of the Naval Research (ONR).
|Printer friendly Cite/link Email Feedback|
|Author:||Tadi, M.; Sritharan, S.S.|
|Publication:||International Journal of Computational and Applied Mathematics|
|Date:||May 1, 2012|
|Previous Article:||A modified regularization scheme for classification problems.|
|Next Article:||One step hybrid multistep method with two off-step points for solution of second order ordinary differential equations.|