Numerical analysis of eikonal equation

I Introduction

In this article, we consider the approach to transform the eikonal equations to the ODE system. The first part of the article describes in detail all the mathematical calculations. In the second part we briefly describe Maxwell and Luneberg lenses, and explain the approach to their numerical modeling, which allows to obtain ray trajectories and wave fronts from sources of different shapes.

Ii Application of the characteristics method to the eikonal equation solution

ii.1 The eikonal equation

The eikonal equation can be obtained from Maxwell’s equations, written for the regions free of currents and charges, and under the condition of a time-changing harmonic electromagnetic field in a nonconducting isotropic medium born-wolf:principles_optics::en; stratton:1948::en; ll:2::en; ll:8::en

. In general, the eikonal equation is written as a partial differential equation of the first order:

{\begin{matrix} | \nabla u (r) |^{2} = n^{2} (r), r \in R^{3}, u (r) = φ (r), x \in Γ \subset R^{3} . \end{matrix}

where $r = (x, y, z)^{T}$

is radius-vector,

φ (r)

is the boundary condition,

n (r)

is the refractive index of the medium. The function

u (r)

is the real scalar function with a physical meaning of time. It is also often called the eikonal function bruns:1895; klein:1901:eikonal.

For visualization of lens modeling results, we will consider their projection on the $O x y$ plane. In this case, the eikonal equation is reduced to the following two-dimensional form:

⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ \begin{matrix} {(\frac{\partial u (x, y)}{\partial x})}^{2} + {(\frac{\partial u (x, y)}{\partial y})}^{2} = n^{2} (x, y), (x, y) \in R^{2}, u (x, y) = φ (x, y), (x, y) \in Γ \subset R^{2} . \end{matrix}

(1)

Using the method of characteristics, the eikonal equation can be transformed into an ODE system that can be solved by standard numerical methods.

ii.2 Characteristics of the eikonal equation

Let us briefly describe the method of characteristics jeong:2007:eikonal_parallel; kimmel:1998:geodesic; zhao:2004:fsm; beliakov:1996:numerical_luneburg; gremaud:2006:fsm_eikonal; bak:2010:fast-sweeping-method and the application of this method to the eikonal equation.

The partial differential equation of the following form is considered:

a_{1} (x, y) \frac{\partial u (x, y)}{\partial x} + a_{2} (x, y) \frac{\partial u (x, y)}{\partial y} = f (x, y),

(2)

where $a_{1} (x, y)$ , $a_{2} (x, y)$ , $u (x, y)$ and $f (x, y)$ are sufficiently smooth functions. This equation is equivalent to the statement that a vector field with components $a_{1} (x, y)$ , $a_{2} (x, y)$ , $f (x, y)$ is tangent to the surface $z = u (x, y)$ , which has a normal vector with components $(u_{x}, u_{y}, - 1)$ . Accordingly, for this equation one can write the system of ODE, called equations of characteristics. It has the following form:

\frac{d x}{d t} = a_{1} (x, y), \frac{d y}{d t} = a_{2} (x, y), \frac{d u (x, y)}{d t} = f (x, y) .

This ODE system reduces the solution of the partial differential equation of the first order to the solution of the ODE system of the first order.

To get the equations of characteristics for the eikonal equation one has to perform two steps. At the first stage, the equation should be converted to the form (2), and after that the ODE system may be written down. For the two-dimensional case the conversion of the eikonal equation to (2) is performed by replacing

p_{1} = \frac{\partial u}{\partial x}, p_{2} = \frac{\partial u}{\partial y} .

In this case, the equation itself is converted to form:

| p |^{2} = p_{1}^{2} + p_{2}^{2} = n^{2} (x, y) .

A number of changes should be made.

	$\frac{\partial}{\partial x} (p_{1}^{2} + p_{2}^{2}) = 2 p_{1} \frac{\partial p_{1}}{\partial x} + 2 p_{2} \frac{\partial p_{2}}{\partial x} = 2 n \frac{\partial n}{\partial x},$
	$\frac{\partial}{\partial y} (p_{1}^{2} + p_{2}^{2}) = 2 p_{1} \frac{\partial p_{1}}{\partial y} + 2 p_{2} \frac{\partial p_{2}}{\partial y} = 2 n \frac{\partial n}{\partial y} .$

After that the following system of equations is obtained:

\begin{matrix} p_{1} \frac{\partial p_{1}}{\partial x} + p_{2} \frac{\partial p_{2}}{\partial x} = n \frac{\partial n}{\partial x}, p_{1} \frac{\partial p_{1}}{\partial y} + p_{2} \frac{\partial p_{2}}{\partial y} = n \frac{\partial n}{\partial y}, \end{matrix} ⟹ \begin{matrix} (p, \frac{\partial p}{\partial x}) = n \frac{\partial n}{\partial x}, (p, \frac{\partial p}{\partial y}) = n \frac{\partial n}{\partial y} . \end{matrix}

Since

\frac{\partial p_{1}}{\partial y} = \frac{\partial^{2} u (x, y)}{\partial y \partial x} = \frac{\partial^{2} u (x, y)}{\partial x \partial y} = \frac{\partial p_{2}}{\partial x},

then

\frac{\partial p_{1}}{\partial y} = \frac{\partial p_{2}}{\partial x} .

Using this equality our expressions may be converted in

		$\frac{\partial p}{\partial x} = (\frac{\partial p_{1}}{\partial x}, \frac{\partial p_{2}}{\partial x}) = (\frac{\partial p_{1}}{\partial x}, \frac{\partial p_{1}}{\partial y}) = \frac{\partial p_{1}}{\partial x} = \nabla p_{1},$
		$\frac{\partial p}{\partial y} = (\frac{\partial p_{1}}{\partial y}, \frac{\partial p_{2}}{\partial y}) = (\frac{\partial p_{2}}{\partial x}, \frac{\partial p_{2}}{\partial y}) = \frac{\partial p_{2}}{\partial x} = \nabla p_{2} .$

As a result:

Thus, the goal is achieved — the equation (1) is transformed in two equations of the form (2).

\begin{matrix} p_{1} \frac{\partial p_{1}}{\partial x} + p_{2} \frac{\partial p_{1}}{\partial y} = n \frac{\partial n}{\partial x}, p_{1} \frac{\partial p_{2}}{\partial x} + p_{2} \frac{\partial p_{2}}{\partial y} = n \frac{\partial n}{\partial y}, \end{matrix} ⟹ ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ \begin{matrix} \frac{p_{1}}{n^{2}} \frac{\partial p_{1}}{\partial x} + \frac{p_{2}}{n^{2}} \frac{\partial p_{1}}{\partial y} = \frac{1}{n} \frac{\partial n}{\partial x}, \frac{p_{1}}{n^{2}} \frac{\partial p_{2}}{\partial x} + \frac{p_{2}}{n^{2}} \frac{\partial p_{2}}{\partial y} = \frac{1}{n} \frac{\partial n}{\partial y} . \end{matrix}

The characteristics for each equation may be written down:

$\frac{p_{1}}{n^{2}} \frac{\partial p_{1}}{\partial x} + \frac{p_{2}}{n^{2}} \frac{\partial p_{1}}{\partial y} = \frac{1}{n} \frac{\partial n}{\partial x}$	$\frac{p_{1}}{n^{2}} \frac{\partial p_{2}}{\partial x} + \frac{p_{2}}{n^{2}} \frac{\partial p_{2}}{\partial y} = \frac{1}{n} \frac{\partial n}{\partial y}$
$\frac{d x}{d t} = \frac{p_{1}}{n^{2}}$	$\frac{d x}{d t} = \frac{p_{1}}{n^{2}}$
$\frac{d y}{d t} = \frac{p_{2}}{n^{2}}$	$\frac{d y}{d t} = \frac{p_{2}}{n^{2}}$
$\frac{d p_{1}}{d t} = \frac{1}{n} \frac{\partial n}{\partial x}$	$\frac{d x}{d t} = \frac{1}{n} \frac{\partial n}{\partial y}$

So the ODE system of four equations with four functions: $x (t)$ , $y (t)$ , $p_{1} (t)$ , $p_{2} (t)$ , is derived:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ \begin{matrix} \frac{d x}{d t} = \frac{p_{1}}{n^{2}}, \frac{d y}{d t} = \frac{p_{2}}{n^{2}}, \frac{d p_{1}}{d t} = \frac{1}{n} \frac{\partial n}{\partial x}, \frac{d p_{2}}{d t} = \frac{1}{n} \frac{\partial n}{\partial y} . \end{matrix}

The initial conditions:

		${x (t) \|}_{t = 0} = x_{0},$
		${y (t) \|}_{t = 0} = y_{0},$
		${p_{1} (t) \|}_{t = 0} = c_{1} n (x_{0}, y_{0}),$
		${p_{2} (t) \|}_{t = 0} = c_{2} n (x_{0}, y_{0}) .$

Constants $c_{1}$ and $c_{2}$ are bonded by following relation $c_{1}^{2} + c_{2}^{2} = 1$ . These constants may be presented as $c_{1} = cos (α)$ and $c_{2} = sin (α)$ . The initial conditions give a mathematical description of the source of the rays. For example, to model a point source, we need to fix the initial coordinates $x_{0}$ , $y_{0}$ and change the angle $α$ , which will set the angle of the beam exit from the source-point. To simulate the radiating surface, on the contrary, it is necessary to fix the angle $α$ and change the coordinates $x_{0}$ and $y_{0}$ .

Let us to find the relation between the parameter $t$ and the function $u (x, y)$ . Since:

\frac{d u}{d t} = \frac{\partial u}{\partial x} \frac{d x}{d t} + \frac{\partial u}{\partial y} \frac{d y}{d t},

and

(\frac{\partial u}{\partial x}, \frac{\partial u}{\partial y}) = \nabla u = p,

when

\frac{d u}{d t} = \nabla u \frac{d x}{d t} = (p, \frac{d x}{d t}) = p_{1} \frac{d x}{d t} + p_{2} \frac{d y}{d t} = \frac{p_{1} p_{1}}{n^{2}} + \frac{p_{2} p_{2}}{n^{2}} = \frac{| p |^{2}}{n^{2}} .

Due to the fact that $| p |^{2} = n^{2} (x, y)$ , we obtain:

\frac{d u}{d t} = \frac{| p |^{2}}{n^{2}} = \frac{n^{2}}{n^{2}} = 1 \Rightarrow \frac{d u}{d t} = 1.

The solution of the equation $u_{t} = 1$ is the function $u (x, y) = t + c o n s t$ which implies that the parameter $t$ has a physical meaning of the signal propagation time from the point $(x_{0}, y_{0})$ to the point $(x, y)$

In polar coordinates, the eikonal equation has the following form::

{(\frac{\partial u (r, φ)}{\partial r})}^{2} + \frac{1}{r^{2}} {(\frac{\partial u (r, φ)}{\partial φ})}^{2} = n^{2} (r),

and the corresponding system of ODEs will have the form:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ \begin{matrix} \frac{d r}{d t} = p_{r}, \frac{d φ}{d t} = \frac{p_{φ}}{r}, \frac{d p_{r}}{d t} = n \frac{\partial n}{\partial r} + \frac{p_{φ}^{2}}{r}, \frac{d p_{φ}}{d t} = - \frac{p_{φ} p_{r}}{r} . \end{matrix}

The initial conditions:

		${r (t) \|}_{t = 0} = r_{0},$
		${φ (t) \|}_{t = 0} = φ_{0},$
		${p_{r} (t) \|}_{t = 0} = c_{1} n (r_{0}),$
		${p_{φ} (t) ∣ ∣}_{t = 0} = c_{2} n (r_{0}) .$

Iii Numerical simulation of Luneburg and Maxwell lenses

Let’s consider the examples of lenses kulyabov:2017:sfm:geometrization_maxwell; kulyabov:2018:sfm:lens-calculations.

iii.1 Luneburg lens

Luneburg lens luneburg:1964; morgan:1958:luneberg_lens; lock:2008:luneburg_ray; lock:2008:luneburg_wave is a spherical lens of radius $R$ with the center at point $(X_{0}, Y_{0})$ (consider the projection on the plane Oxy) with a refractive index of the following form

n (x, y) = ⎧ ⎪ ⎨ ⎪ ⎩ \begin{matrix} n_{0} \sqrt{2 - (\frac{r}{R})^{2}}, r ⩽ R, n_{0}, r > R, \end{matrix}

where $r (x, y) = \sqrt{(x - X_{0})^{2} + (y - Y_{0})^{2}}$ is the distance from the center of the lens to an arbitrary point in the $(x, y)$ plane. The formula implies that the coefficient $n$ continuously varies from $n_{0} \sqrt{2}$ to $n_{0}$ starting from the center of the lens and ending with its boundary. The refractive index of the medium outside the lens is constant and is equal to $n_{0}$ . Usually $n_{0}$ is equal to $1$ .

To solve the eikonal equation by the method of characteristics it is necessary to find partial derivatives of the function $n (x, y)$ . For the case of Luneburg lens the partial derivatives are:

\frac{\partial n (x, y)}{\partial x} = - \frac{n_{0}^{2} (x - X_{0})}{R^{2} n (x, y)}, \frac{\partial n (x, y)}{\partial y} = - \frac{n_{0}^{2} (y - Y_{0})}{R^{2} n (x, y)}, r ⩽ R .

Outside the lens region derivatives are equal to $0$ .

iii.2 Maxwell fish eye lens

Maxwell fish eye lens maxwell:1854:fish-eye is also a spherical lens of radius $R$ with the center at point $(X_{0}, Y_{0})$ (consider the projection on the plane Oxy) with a refractive index of the following form:

n (x, y) = ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ \begin{matrix} \frac{n_{0}}{1 + (\frac{r}{R})^{2}}, r ⩽ R, n_{0}, r > R . \end{matrix}

To solve the eikonal equation by the method of characteristics it is necessary to find partial derivatives of the function $n (x, y)$ . For the case of Maxwell lens partial derivatives have the form:

\frac{\partial n (x, y)}{\partial x} = - \frac{2 n^{2} (x, y) (x - X_{0})}{n_{0} R^{2}}, \frac{\partial n (x, y)}{\partial y} = - \frac{2 n^{2} (x, y) (y - Y_{0})}{n_{0} R^{2}}, r ⩽ R .

iii.3 Description of the numerical modeling

Julia programming language joshi:book:learning-julia is used to simulate the trajectories of rays through the Maxwell and Luneburg lenses. We use classical Runge–Kutta methods with constant step to solve the ODE system.

We carry on numerical modeling for lenses with a radius $R = 1$ , the refractive index of the external medium $n_{0} = 1$ , the center of the lens was placed in the point $(X_{0}, Y_{0}) = (2, 0)$ , the boundary region was set as the rectangle $x_{min} = 0$ , $x_{max} = 5$ , $y_{min} = - 1.5$ and $y_{max} = 1.5$ . The point source was placed on the lens boundary at $(x_{0}, y_{0}) = (0, 0)$ . $50$ values of the $α$ parameter have been taken from the interval $[- π / 2 + π / 100, π / 2 - π / 100]$ , which allowed to simulate rays trajectories from a point source within an angle slightly smaller than $180^{\circ}$ . The $t$ parameter was changed within the $[0, 5]$ interval.

Each $α$ parameter value sets new initial conditions for the ODE system. The process of numerical simulation consists in multiple solution of this system for different initial conditions. The numerical solution of the ODE system for a particular initial condition gives us a set of points $(x_{I}, y_{I})$ , $I = 1, \dots, N$ approximating the trajectories of a particular beam. After performing calculations for all the selected initial conditions, we obtain a set of rays. To visualize the rays, it is enough to depict each of the obtained numerical solutions. The result of the simulation can be seen in the Fig. 2 and Fig. 3 (the trajectories of the rays) and Fig. 2 and Fig. 4 (the wavefronts).

Figure 1: The trajectories of the rays in case of Maxwell’s lens for a point source and $n_{0} = 1$

Figure 3: The trajectories of the rays in case of Luneburg lens for a point source and $n_{0} = 1$ .

Figure 4: The wavefronts in case of Luneburg lens for a point source and $n_{0} = 1$ .

To visualize the wave fronts with the resulting numerical data it is necessary to carry out additional recalculations. From each numerical solution, we must select points $(x_{I}, y_{I})$ that correspond to a specific point in time $t_{I}$ .

The use of a numerical method with a fixed step gives an advantage, since each numerical solution will be obtained for the same uniform grid $t_{0} < t_{1} < \dots < t_{i} < \dots < t_{n}$ .

Iv Conclusion

The paper presents the description of the numerical solution of the eikonal equation for the case of Luneburg and Maxwell lenses. The results are visualized as trajectories of rays passing through lenses and as fronts of electromagnetic waves.

Acknowledgements.

The publication has been prepared with the support of the ‘‘RUDN University Program 5-100’’ and funded by Russian Foundation for Basic Research (RFBR) according to the research project No 19-01-00645.

		${x (t) \|}_{t = 0} = x_{0},$
		${y (t) \|}_{t = 0} = y_{0},$
		${p_{1} (t) \|}_{t = 0} = c_{1} n (x_{0}, y_{0}),$
		${p_{2} (t) \|}_{t = 0} = c_{2} n (x_{0}, y_{0}) .$

		${r (t) \|}_{t = 0} = r_{0},$
		${φ (t) \|}_{t = 0} = φ_{0},$
		${p_{r} (t) \|}_{t = 0} = c_{1} n (r_{0}),$
		${p_{φ} (t) ∣ ∣}_{t = 0} = c_{2} n (r_{0}) .$

Numerical analysis of eikonal equation

Authors

A geometric optics ansatz-based plane wave method for two dimensional Helmholtz equations with variable wave numbers

Towards Accuracy and Scalability: Combining Isogeometric Analysis with Deflation to Obtain Scalable Convergence for the Helmholtz Equation

Relative Schwarz Method

Uniform observability of the one-dimensional wave equation for non-cylindrical domains. Application to the control's support optimization

Solitary water wave interactions for the Forced Korteweg-de Vries equation

Numerical study of the stabilization of 1D locally coupled wave equations

A Flawed Dataset for Symbolic Equation Verification

I Introduction

Ii Application of the characteristics method to the eikonal equation solution

ii.1 The eikonal equation

ii.2 Characteristics of the eikonal equation

Iii Numerical simulation of Luneburg and Maxwell lenses

iii.1 Luneburg lens

iii.2 Maxwell fish eye lens

iii.3 Description of the numerical modeling

Iv Conclusion

Acknowledgements.

References

Numerical analysis of eikonal equation

Authors

Related Research

A geometric optics ansatz-based plane wave method for two dimensional Helmholtz equations with variable wave numbers

Towards Accuracy and Scalability: Combining Isogeometric Analysis with Deflation to Obtain Scalable Convergence for the Helmholtz Equation

Relative Schwarz Method

Uniform observability of the one-dimensional wave equation for non-cylindrical domains. Application to the control's support optimization

Solitary water wave interactions for the Forced Korteweg-de Vries equation

Numerical study of the stabilization of 1D locally coupled wave equations

A Flawed Dataset for Symbolic Equation Verification

This week in AI

I Introduction

Ii Application of the characteristics method to the eikonal equation solution

ii.1 The eikonal equation

ii.2 Characteristics of the eikonal equation

Iii Numerical simulation of Luneburg and Maxwell lenses

iii.1 Luneburg lens

iii.2 Maxwell fish eye lens

iii.3 Description of the numerical modeling

Iv Conclusion

Acknowledgements.

References