Numerical analysis of eikonal equation

TL;DR

This paper explores the numerical analysis of the eikonal equation, derived from wave optics, and applies geometric methods to analyze optical devices like Maxwell and Luneburg lenses.

Contribution

It demonstrates the transformation of the eikonal equation into an ODE system for optical device analysis using geometric methods.

Findings

Successful transformation of the eikonal equation into an ODE system

Application of geometric methods to Maxwell and Luneburg lenses

Enhanced understanding of wave and geometric optics connection

Abstract

The Maxwell equations have a fairly simple form. However, finding solutions of Maxwell's equations is an extremely difficult task. Therefore, various simplifying approaches are often used in optics. One such simplifying approach is to use the approximation of geometric optics. The approximation of geometric optics is constructed with the assumption that the wavelengths are small (short-wavelength approximation). The basis of geometric optics is the eikonal equation. The eikonal equation can be obtained from the wave equation (Helmholtz equation). Thus, the eikonal equation relates the wave and geometric optics. In fact, the eikonal equation is a quasi-classical approximation (the Wentzel-Kramers-Brillouin method) of wave optics. This paper shows the application of geometric methods of electrodynamics to the calculation of optical devices, such as Maxwell and Luneburg lenses. The eikonal equation, which was transformed to the ODE system by the method of characteristics, is considered. The resulting system is written for the case of Maxwell and Luneburg lenses.

Equations

$$\{\begin{array}{rl}& |\nabla u(\mathbf{r}){|}^2=n^2(\mathbf{r}),\mathbf{r}\in {\mathbb{R}}^3,\\ & u(\mathbf{r})=\phi (\mathbf{r}),\mathbf{x}\in \Gamma \subset {\mathbb{R}}^3.\end{array}$$

$$\{\begin{array}{rl}& {\left(\frac{\partial u(x,y)}{\partial x}\right)}^2+{\left(\frac{\partial u(x,y)}{\partial y}\right)}^2=n^2(x,y),(x,y)\in {\mathbb{R}}^2,\\ & u(x,y)=\phi (x,y),(x,y)\in \Gamma \subset {\mathbb{R}}^2.\end{array}$$

$$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),$$

$$\frac{\mathrm{d}x}{\mathrm{d}t}=a_1(x,y),\frac{\mathrm{d}y}{\mathrm{d}t}=a_2(x,y),\frac{\mathrm{d}u(x,y)}{\mathrm{d}t}=f(x,y).$$

$$p_1=\frac{\partial u}{\partial x},p_2=\frac{\partial u}{\partial y}.$$

$$|\mathbf{p}{|}^2=p_1^2+p_2^2=n^2(x,y).$$

$$\frac{\partial }{\partial x}(p_1^2+p_2^2)=2p_1\frac{\partial p_1}{\partial x}+2p_2\frac{\partial p_2}{\partial x}=2n\frac{\partial n}{\partial x},$$

$$\frac{\partial }{\partial y}(p_1^2+p_2^2)=2p_1\frac{\partial p_1}{\partial y}+2p_2\frac{\partial p_2}{\partial y}=2n\frac{\partial n}{\partial y}.$$

$$\begin{array}{rl}& 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{array}⟹\begin{array}{rl}& \left(\mathbf{p},\frac{\partial \mathbf{p}}{\partial x}\right)=n\frac{\partial n}{\partial x},\\ & \left(\mathbf{p},\frac{\partial \mathbf{p}}{\partial y}\right)=n\frac{\partial n}{\partial y}.\end{array}$$

$$\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},$$

$$\frac{\partial p_1}{\partial y}=\frac{\partial p_2}{\partial x}.$$

$$\frac{\partial \mathbf{p}}{\partial x}=\left(\frac{\partial p_1}{\partial x},\frac{\partial p_2}{\partial x}\right)=\left(\frac{\partial p_1}{\partial x},\frac{\partial p_1}{\partial y}\right)=\frac{\partial p_1}{\partial \mathbf{x}}=\nabla p_1,$$

$$\frac{\partial \mathbf{p}}{\partial y}=\left(\frac{\partial p_1}{\partial y},\frac{\partial p_2}{\partial y}\right)=\left(\frac{\partial p_2}{\partial x},\frac{\partial p_2}{\partial y}\right)=\frac{\partial p_2}{\partial \mathbf{x}}=\nabla p_2.$$

$$\begin{array}{rl}& \left(\mathbf{p},\frac{\partial \mathbf{p}}{\partial x}\right)=n\frac{\partial n}{\partial x},\\ & \left(\mathbf{p},\frac{\partial \mathbf{p}}{\partial y}\right)=n\frac{\partial n}{\partial y},\end{array}⟹\begin{array}{rl}& \left(\mathbf{p},\nabla p_1\right)=n\frac{\partial n}{\partial x},\\ & \left(\mathbf{p},\nabla p_2\right)=n\frac{\partial n}{\partial y}.\end{array}$$

$$\begin{array}{rl}& 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{array}⟹\{\begin{array}{rl}& \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{array}$$

$$\{\begin{array}{rl}& \frac{\mathrm{d}x}{\mathrm{d}t}=\frac{p_1}{n^2},\\ & \frac{\mathrm{d}y}{\mathrm{d}t}=\frac{p_2}{n^2},\\ & \frac{\mathrm{d}{p}_1}{\mathrm{d}t}=\frac{1}{n}\frac{\partial n}{\partial x},\\ & \frac{\mathrm{d}{p}_2}{\mathrm{d}t}=\frac{1}{n}\frac{\partial n}{\partial y}.\end{array}$$

$${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).$$

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

$$\left(\frac{\partial u}{\partial x},\frac{\partial u}{\partial y}\right)=\nabla u=\mathbf{p},$$

$$\frac{\mathrm{d}u}{\mathrm{d}t}=\nabla u\frac{\mathrm{d}\mathbf{x}}{\mathrm{d}t}=\left(\mathbf{p},\frac{\mathrm{d}\mathbf{x}}{\mathrm{d}t}\right)=p_1\frac{\mathrm{d}x}{\mathrm{d}t}+p_2\frac{\mathrm{d}y}{\mathrm{d}t}=\frac{p_1{p}_1}{n^2}+\frac{p_2{p}_2}{n^2}=\frac{|\mathbf{p}{|}^2}{n^2}.$$

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

$${\left(\frac{\partial u(r,\phi )}{\partial r}\right)}^2+\frac{1}{r^2}{\left(\frac{\partial u(r,\phi )}{\partial \phi }\right)}^2=n^2(r),$$

$$\{\begin{array}{rl}& \frac{\mathrm{d}r}{\mathrm{d}t}=p_r,\\ & \frac{\mathrm{d}\phi }{\mathrm{d}t}=\frac{p_{\phi }}{r},\\ & \frac{\mathrm{d}{p}_r}{\mathrm{d}t}=n\frac{\partial n}{\partial r}+\frac{p_{\phi }^2}{r},\\ & \frac{\mathrm{d}{p}_{\phi }}{\mathrm{d}t}=-\frac{p_{\phi }{p}_r}{r}.\end{array}$$

$${r(t)|}_{t=0}=r_0,$$

$${\phi (t)|}_{t=0}={\phi }_0,$$

$${p_r(t)|}_{t=0}=c_{1}n(r_0),$$

$${p_{\phi }(t)|}_{t=0}=c_{2}n(r_0).$$

$$n(x,y)=\left\{\begin{aligned} &n_{0}\sqrt{2-\Big{(}\dfrac{r}{R}\Big{)}^{2}},\quad r\leqslant R,\\ &n_{0},\quad r>R,\end{aligned}\right.$$

$$\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.$$

$$n(x,y)=\left\{\begin{aligned} &\dfrac{n_{0}}{1+\Big{(}\dfrac{r}{R}\Big{)}^{2}},\quad r\leqslant R,\\ &n_{0},\quad r>R.\end{aligned}\right.$$

$$\frac{\partial n(x,y)}{\partial x}=-\frac{2n^2(x,y)(x-X_0)}{n_0{R}^2},\frac{\partial n(x,y)}{\partial y}=-\frac{2n^2(x,y)(y-Y_0)}{n_0{R}^2},r⩽R.$$

Full text

eikonal equation, Luneburg lens, Maxwell lens, characteristics method, Julia

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:

[TABLE]

where $\mathbf{r}=(x,y,z)^{T}$ is radius-vector, $\varphi(\mathbf{r})$ is the boundary condition, $n(\mathbf{r})$ is the refractive index of the medium. The function $u(\mathbf{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 $Oxy$ plane. In this case, the eikonal equation is reduced to the following two-dimensional form:

[TABLE]

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:

[TABLE]

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:

[TABLE]

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

[TABLE]

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

[TABLE]

A number of changes should be made.

[TABLE]

After that the following system of equations is obtained:

[TABLE]

Since

[TABLE]

then

[TABLE]

Using this equality our expressions may be converted in

[TABLE]

As a result:

[TABLE]

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

[TABLE]

The characteristics for each equation may be written down:

[TABLE]

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

[TABLE]

The initial conditions:

[TABLE]

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(\alpha)$ and $c_{2}=\sin(\alpha)$. 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 $\alpha$, 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 $\alpha$ 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:

[TABLE]

and

[TABLE]

when

[TABLE]

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

[TABLE]

The solution of the equation $u_{t}=1$ is the function $u(x,y)=t+\mathrm{const}$ 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::

[TABLE]

and the corresponding system of ODEs will have the form:

[TABLE]

The initial conditions:

[TABLE]

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

[TABLE]

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:

[TABLE]

Outside the lens region derivatives are equal to [math].

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:

[TABLE]

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:

[TABLE]

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 $\alpha$ parameter have been taken from the interval $[-\pi/2+\pi/100,\pi/2-\pi/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 $\alpha$ 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,\ldots,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).

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}<\ldots<t_{i}<\ldots<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.