Advanced Methods for Solving Nonlinear Eigenvalue Problems of Generalized Phase Optimization

Abstract

In the process of solving the problems of generalized phase optimization the necessity to apply an eigenvalue approach often appears. The practical statement of the optimization problems consists of using the amplitude characteristics of functions that are sought. The usual way of optimization is deriving the Euler equation of the functional, which is used as criterion of optimization. As a rule, such equation is an integral one. It is worth pointing out that the integral equations of the generalized phase optimization are nonlinear ones. The characteristic property of such equations is non-uniqueness of solutions and their branching or bifurcation. The determination of branching solutions leads to the investigation of the corresponding homogeneous equations and the respective eigenvalue problem. This problem is nonlinear because of specificity of the statement of the optimization problem. The study of the above problem allows us to determine a set of points, in which the respective eigenvalues are equal to unity that determines the branching points of solutions. The data of calculations testify to the ability of the approach proposed to determine the solutions of nonlinear equations numerically with not large computations.

1. Introduction

The nonlinear eigenvalue approach is used in this chapter for the study of the properties to solutions of the generalized phase problem related to the synthesis of radiation systems through the incomplete data. Such incompleteness is considered here in the example of an indeterminate phase characteristic of function, which characterizes the radiation of the plane antenna arrays.

The problems with an indeterminate phase of the wave field arise in various applications and are widely described in the literature. The most well-known of these is the so-called phase problem (see, for example, [1, 2, 3, 4]). It consists in restoring the phase distribution (argument) of the Fourier transform of a finite function by its amplitude (module) given (measured) along the entire real axis. This problem belongs to the classical problems of recovery (identification) and requires the conditions of existence of a unique solution.

In this chapter, another class of inverse problems is considered, and it can be termed optimization (design) problems. In sense of the Fourier transform, this can be, for example, the problem of finding such a finite complex function, the modulus of its Fourier transform satisfies a certain requirement (e.g., is close to a given positive function). As a rule, such requirements are formulated in the variational form, as the minimization of certain functionals. Obviously, such a formulation does not require a uniform solution. On the contrary, the existence of many solutions is often desirable because the above allows many degrees of freedom to determine an appropriate solution. The characteristic applications of such phase optimization problems include the theory of power transmission lines, field converters, antennas and resonators. The first works dealing with nonlinear inverse problems of such type appeared in the second part of the last century (see, for instance [5, 6, 7, 8, 9, 10]).

In mathematical terms, problems of this type are reduced to the nonlinear integral equations of Hammerstein type [11, 12, 13]. They contain a linear kernel and a nonlinear multiplier that depends on a complex unknown function as an integrand. As a rule, the argument (phase) of this function appears there separately from the module. Similar equations are found in the literature in the context of the mentioned phase problem [14, 15]. They have different solutions, and the study of their structure and process of branching or bifurcation is an interesting mathematical problem [16].

Due to their nonlinearity, the problems under consideration require the development and application of special analytical and numerical methods for their solving. Along with the iterative methods that simulate the physical processes of field formation, the various modifications of Newton’s method could be the most promising in this direction [17]. One such modification, which uses solving the nonlinear eigenvalue problems and searching for the zero curves of respective determinants, is proposed in this Chapter. It allows simultaneously with the finding of the branch of solutions to detect the presence of branching points on it and to determine them approximately, provided by this the initial approximations for more accurate calculation.

The nonlinear eigenvalue problems arise in pure and applied mathematics, as well as in the different areas of science that investigate the nonlinear phenomena [18, 19]. A variety of analytical-numerical methods have been elaborated till now for solving the nonlinear problems in acoustics, electrodynamics, fluid dynamics and other areas of applied science [20, 21]. The methods, developed until that time, were focused mainly on solving one-dimensional problems. The difficulties of analytical and computational nature appear if to apply them to a multidimensional problem. The method of implicit function is one of effective tools that been applied for solving the two- and three-dimensional nonlinear eigenvalue problem in the last two decades [22, 23, 24]. The extension of this method, which leads to solving the Cauchy problem (21) and (22), we apply in Section 3 to solve the nonlinear two-dimensional eigenvalue problem.

2. The operators of direct electrodynamics problem

In the physical relation, the radiation system represents the plane array with the rectangular or hexagonal placement of radiators. Firstly, we consider the array with the rectangular ordering of separate elements (Figure 1a).

Figure 1.

The placement of elements for the considered arrays. a) rectangular. b) hexagonal.

Consider a plane rectangular array consisting of $N_2\cdot M_2=\left(2N+1\right)\cdot \left(2M+1\right)$ identical elements (radiators), which are located in the $\mathit{xOy}$ plane of the Cartesian coordinate system equidistantly for each of the coordinates. Since the radiators are identical, it is possible to formulate the synthesis problem not for the whole three-dimensional vector directivity pattern (DP), but only for some complex scalar function $f\left(x_1,x_2\right)$ that is termed as the array multiplier. This function for a rectangular equidistant array has the form [25]:

$$f\left(x_1,x_2\right)=A\mathbf{I}\equiv \sum _{n=-N}^N\sum _{m=-M}^M{I}_{\mathit{nm}}{e}^{i\left(c_1{\mathit{nx}}_1+c_2{\mathit{mx}}_2\right)}$$E1

where $\mathbf{I}=\left\{I_{\mathit{nm}},-N\le n\le N,-M\le m\le M\right\}$ is a set of excitations (currents) in the array’s elements, $x_1=sin\theta cos\phi /sin{\alpha }_1$, $x_2=sin\theta sin\phi /sin{\alpha }_2$ are the generalized angular coordinates, $c_1={\mathit{kd}}_{1}sin{\alpha }_1$, $c_2={\mathit{kd}}_{2}sin{\alpha }_2$, $k=2\pi /\lambda$ is wave number, $d_1$ and $d_2$ are the distances between radiators along the $\mathit{Ox}$ axis and $\mathit{Oy}$ axis respectively, ${\alpha }_1$ and ${\alpha }_2$ are the angular coordinates, within which the desired power DP $P\left(x_1,x_2\right)$ is not equal to zero ($P\left(x_1,x_2\right)\equiv 0$ outside these angles). The function $f\left(x_1,x_2\right)$ possesses $2\pi /c_1-$periodicity with respect to $x_1$ and $2\pi /c_2-$ periodicity with respect to $x_2$. Let us denote the region of change of coordinates $x_1$ and $x_2$ on one period as $\mathrm{\Omega }=\left\{\left(x_1,x_2\right):|x_1|\le \pi /c_1,|x_2|\le \pi /c_2\right\}$. Below, the function $f\left(x_1,x_2\right)$ is termed as the DP of array.

A similar formula can be derived for the array with the hexagonal placement of separate elements (Figure 1b)

$$f\left(x_1,x_2\right)=A\mathbf{I}\equiv \sum _{m=-M_2}^{M_2}\sum _{n=-N_1\left(m\right)}^{N_1\left(m\right)}{I}_{\mathit{nm}}{e}^{i\left(c_1{\mathit{nx}}_1+c_2{\mathit{mx}}_2\right)}$$E2

where $M=2M_2+1$ is quantity of the linear subarrays, then $N=2N_1\left(m\right)+1$ is the number of elements in the $m-$th subarray.

Eqs. (1) and (2) for DP $f\left(x_1,x_2\right)$ represent the result of a linear operator $A$, which acts on a complex-valued space $H_I=C^{N_2\times M_2}$ (rectangular case) or $H_I=C^{N_0\times M}$ (hexagonal case) to the space of complex functions of two variables defined in the domain $\mathrm{\Omega }$. The value $N_0$ determines the number of elements in the central linear subarray in the hexagonal case.

Assume that the desired power DP $P\left(s_1,s_2\right)$ is not equal to zero in some regions $\overline{G}\subseteq \mathrm{\Omega }$, and it is equal to zero outside. The optimization problem is formulated as the minimization problem of the functional

$${\sigma }_{\alpha }\left(\mathbf{I}\right)={‖P-{\left|\mathit{AI}\right|}^2‖}_f^2+\alpha {‖\mathbf{I}‖}_I^2$$E3

where ${‖\cdot ‖}_f$ and ${‖\cdot ‖}_I$ determine the norms in the space of DPs and space of currents respectively, which are defined by the inner products

$${‖f‖}_f^2={\left(f_1,f_2\right)}_f=\underset{\mathrm{\Omega }}{\iint }{f}_1\left(x_1,x_2\right){\overline{f}}_2\left(x_1,x_2\right){\mathit{dx}}_1{\mathit{dx}}_2$$E4

$${‖\mathbf{I}‖}_I^2={\left({\mathbf{I}}_1,{\mathbf{I}}_2\right)}_I=\frac{4{\pi }^2}{c_1{c}_2}\sum _{n=-N}^N\sum _{m=-M}^M{I}_{1\mathit{nm}}{\overline{I}}_{2\mathit{nm}}$$E5

here the values ${\overline{f}}_2\left(x_1,x_2\right)$ and ${\overline{I}}_{2\mathit{nm}}$ are conjugated to $f_2\left(x_1,x_2\right)$ and $I_{2\mathit{nm}}$.

The nonlinear integral equation for the complex vector $\mathbf{I}$of currents in space $H_{\mathbf{I}}$, which is derived using the necessary condition of the minimum of functional (3), has the form [26].

$$\alpha \mathbf{I}+2A^{\ast }\left(|A\mathbf{I}|A\mathbf{I}\right)-2A^{\ast }\left(\mathit{PA}\mathbf{I}\right)=0$$E6

Here $A^{\ast }$ is the operator adjoint to $A$, its form is defined by equality ${\left(A\mathbf{I},f\right)}_f={\left(\mathbf{I},A^{\ast }f\right)}_I$. Using the inner products (4), (5) and Eq. (1) we obtain.

$$\begin{array}{l}{\left(A^{\ast }f\right)}_{\mathit{nm}}=\frac{c_1{c}_2}{4{\pi }^2}\underset{\mathrm{\Omega }}{\iint }f\left(x_1,x_2\right)e^{-i\left(c_1{\mathit{nx}}_1+c_2{\mathit{mx}}_2\right)}{\mathit{dx}}_1{\mathit{dx}}_2,n=-N,-N+1\dots ,\\ N-1,N,m=-M,-M+1,\dots ,M-1,M.\end{array}$$E7

If to act by operator $A$ on both the parts of (6), we get a nonlinear integral equation of Hammerstein type for the function $f$

$$\mathit{\alpha f}+2{\mathit{AA}}^{\ast }\left(|f|f\right)-2{\mathit{AA}}^{\ast }\left(\mathit{Pf}\right)=0$$E8

The kernel of the ${\mathit{AA}}^{\ast }$ operator for the rectangular array is defined as

$$K\left(c_1,c_2,x_1,x_1\text{'},x_2,x_2\text{'}\right)=K_1\left(c_1,x_1,x_1\text{'}\right)K_2\left(c_2,x_2,x_2\text{'}\right),$$E9

where

$$K_1\left(x_1,x_1^{\text{'}},c_1\right)=\frac{c_1}{\pi }\frac{sin\left(N_2{c}_1\left(x_1-x_1^{\text{'}}\right)/2\right)}{sin\left(c_1\left(x_1-x_1^{\text{'}}\right)/2\right)}$$E10

$$K_2\left(x_2,x_2^{\text{'}},c_2\right)=\frac{c_2}{\pi }\frac{sin\left(M_2{c}_2\left(x_2-x_2^{\text{'}}\right))/2\right)}{sin\left(c_2\left(x_2-x_2^{\text{'}}\right))/2\right)}$$E11

The kernel of the ${\mathit{AA}}^{\ast }$ operator for the hexagonal array is more complicated because we can not to present it in the form of two multipliers

$$\begin{array}{l}K\left(c_1,c_2,x_1,x_1^{\text{'}},x_2,x_{12}^{\text{'}}\right)=\frac{sin\left[c_1\left(N_1\left(0\right)-1/2\right)\left(x_1-x_1^{\text{'}}\right)\right]}{sin\left(1/2c_1\left(x_1-x_1^{\text{'}}\right)\right)}+\\ +2\sum _{m=1}^{M_2}cos{\mathit{mc}}_2\left(x_2-x_2^{\text{'}}\right)\{\begin{array}{c}\frac{sin\left[c_1\left(N_1\left(m\right)-1/2\right)\left(x_1-x_1^{\text{'}}\right)\right]}{sin\left(1/2c_1\left(x_1-x_1^{\text{'}}\right)\right)},N_1\left(m\right)\text{is}\mathrm{odd},\hfill \\ 2\sum _{n=1}^{N_1\left(m\right)}cos\left[c_1\left(n-1/2\right)\left(x_1-x_1^{\text{'}}\right)\right],N_1\left(m\right)\text{is even.}\hfill \end{array}\end{array}$$E12

The kernels (9) and (12) of the integral Eq. (8) are real and degenerate. Since Eqs. (6) and (8) are nonlinear ones, both may have a non-unique solution. The number of solutions and their properties is studied according to the method proposed in [16, 27]. In the practical applications, the solution of Eqs. (6) and (8) is performed by the method of successive approximations. The convergence of the method depends on the parameter $\alpha$, desired DP $P\left(x_1,x_2\right)$, as well as the parameters $c_1$and $c_2$contained in the kernels (9) and (12).

3. Search for the bifurcation curves

We should use the linear integral equation to define the bifurcation curves according to [16]. Based on this equation, we pass to the respective eigenvalue problems, solutions of which allow us to find the characteristic values of parameters $c_1$ and $c_2$ in the kernel of equation, at which the bifurcation appears.

3.1 Description of procedure

The linear equation

$$\mathit{\alpha f}=2{\mathit{AA}}^{\ast }\left(\mathit{Pf}\right)$$E13

is used to study the properties of Eq. (8).

In contrast to a similar equation for the amplitude DP synthesis problem [25], Eq. (8) does not have a trivial nonzero initial solution $f_0$ for all parameters $c_1$ and $c_2$; the trivial solution $f_0$ is zero for it, so in contrast to the problem of synthesis by amplitude DP, we are not talking about the branching of solutions, but about their bifurcation.

The problem of finding bifurcation curves is reduced to the corresponding eigenvalue problem. The equation for eigenfunctions and corresponding eigenvalues, which refers to (13), is

$$g\left(x_1,x_2\right)=2{\mathit{\lambda \alpha }}^{-1}\underset{\mathrm{\Omega }}{\iint }g\left(x_1^{\text{'}},x_2^{\text{'}}\right)K_1\left(c_1,x_1,x_1^{\text{'}}\right)K_2\left(c_2,x_2,x_2^{\text{'}}\right){\mathit{dx}}_1^{\text{'}}{\mathit{dx}}_2^{\text{'}}$$E14

As stated by the branching theory of solutions of the nonlinear equations [16], the bifurcation points can be those values of $c_1$and $c_2$ at which Eq. (14) has nonzero solutions.

Using the properties of the degeneracy of the kernel ${\mathit{AA}}^{\ast }$, we reduce Eq. (14) to the equivalent system of the linear algebraic equations (SLAE). The coefficients of matrix of this equation depend on the parameters $c_1$ and $c_2$ analytically. To this end, the equations for eigenfunctions corresponding to (13) are written as

$$g\left(x_1,x_2\right)=\sum _{n=-N}^N\sum _{m=-M}^M{x}_{\mathit{nm}}{e}^{i\left(c_1{\mathit{nx}}_1+c_2{\mathit{mx}}_2\right)}$$E15

where

$$x_{\mathit{nm}}=\frac{c_1{c}_2}{4{\pi }^2}\underset{\mathrm{\Omega }}{\iint }P\left(x_1^{\text{'}},x_2^{\text{'}}\right)g\left(x_1^{\text{'}},x_2^{\text{'}}\right)e^{-i\left(c_1{\mathit{nx}}_1^{\text{'}}+c_2{\mathit{mx}}_2^{\text{'}}\right)}{\mathit{dx}}_1^{\text{'}}{\mathit{dx}}_2^{\text{'}}$$E16

Multiplying both the parts of (15) on $P\left(x_1^{\text{'}},x_2^{\text{'}}\right)e^{-i\left(c_1{\mathit{kx}}_1^{\text{'}}+c_2{\mathit{lx}}_2^{\text{'}}\right)}$ at $k=-N,-N+1,\dots ,N-1,N,l=-M,-M+1,\dots ,M-1,M$ and integrating over the domain $\mathrm{\Omega }$, we obtain a system of linear algebraic equations to determine the quantities $x_{\mathit{nm}}$

$$\begin{array}{l}{x}_{\mathit{kl}}=\sum _{n=-N}^N\sum _{m=-M}^M{a}_{\mathit{nm}}^{\left(\mathit{kl}\right)}\left(c_1,c_2\right)x_{\mathit{nm}},k=-N,-N+1,\dots ,N-1,N,\\ l=-M,-M+1,\dots ,M-1,M,\end{array}$$E17

where

$$a_{\mathit{nm}}^{\left(\mathit{kl}\right)}=\frac{c_1{c}_2}{4{\pi }^2}\underset{\mathrm{\Omega }}{\iint }P\left(x_1,x_2\right)e^{-i[\left(c_1\left(n-k\right)x_1+c_2\left(m-l\right)x_2\right)}{\mathit{dx}}_1{\mathit{dx}}_2$$E18

and matrix of the coefficients $a_{\mathit{nm}}^{\left(\mathit{kl}\right)}$ is self-adjoint and Hermitian.

Thus, we obtained a two-parameter nonlinear spectral problem corresponding to a homogeneous SLAE (17). This problem can be given as

$$\left(E_M-A_M\left(c_1,c_2\right)\right)\mathbf{x}=0$$E19

where $A_M$ is the matrix of coefficients $a_{\mathit{nm}}^{\left(\mathit{kl}\right)}$, $E_M$ is a unit matrix of dimension $N_2\cdot M_2$.

For the system (19), the equality

$$\mathrm{\Psi }\left(c_1,c_2\right)=det\left[E_M-A_M\left(c_1,c_2\right)\right]=0$$E20

must be met to have a non-zero solution.

One can easy to make sure that the function $\mathrm{\Psi }\left(c_1,c_2\right)$ is real. Moreover, since $A_M\left(c_1,c_2\right)$ is the Hermitian matrix, then $E_M-A_M\left(c_1,c_2\right)$ is Hermitian too. The determinant of the Hermitian matrix is a real number [28]. Thus, $\mathrm{\Psi }\left(c_1,c_2\right)$ is a real function of real arguments $c_1$and $c_2$.

Consequently, the problem to find the eigenvalues of Eq. (14) or to determine the solution of the equivalent SLAE (19) is reduced to finding zeros of function $\mathrm{\Psi }\left(c_1,c_2\right)$.

If to consider the equation $\mathrm{\Psi }\left(c_1,c_2\right)=0$ as a problem of determining an implicit function $c_2=c_2\left(c_1\right)$ in the vicinity of some point $c_1$, we get Cauchy problem [29].

$$\frac{{\mathit{dc}}_2}{{\mathit{dc}}_1}=\frac{{\mathrm{\Psi }}_{c_1}\text{'}\left(c_1,c_2\right)}{{\mathrm{\Psi }}_{c_2}\text{'}\left(c_1,c_2\right)}$$E21

$$c_2\left(c_1^{\left(0\right)}\right)=c_2^{\left(0\right)}$$E22

To retrieve the initial conditions (22) we pass to an auxiliary one-dimensional nonlinear spectral problem if to substitute $c_2$ by $c_2=\gamma c_1$ in Eq. (20) with some real parameter $\gamma$. As a result, we get the one-dimensional eigenvalue problem

$$\left(E_M-A_M\left(c_1,\gamma c_1\right)\right)\tilde{\mathbf{x}}\equiv \left(E_M-{\tilde{A}}_M\left(c_1\right)\right)\tilde{\mathbf{x}}=0$$E23

Eq. (20), which corresponds to Eq. (23), is

$$\mathrm{\Psi }\left(c_1,\gamma c_1\right)=det\left[E_M-{\tilde{A}}_M\left(c_1\right)\right]=0$$E24

Let $c_1^{\left(0\right)}$ be the solution of the Eq. (24), then $\left(c_1^{\left(0\right)},c_2^{\left(0\right)}\right)=\left(c_1^{\left(0\right)},\gamma c_1^{\left(0\right)}\right)$ is the point that corresponds to eigenvalue ${\lambda }_0\approx 1$ of Eq. (15). By solving Eqs. (21) and (22) in a small vicinity of point $\left(c_1^{\left(0\right)},c_2^{\left(0\right)}\right)$, we find the spectral curve of the matrix-function $A_M\left(c_1,c_2\right)$, which is the curve $c_2\left(c_1\right)$ defining a set of the bifurcation points.

The eigenfunctions of Eq. (14) are defined as the eigenvectors of matrix $A_M\left(c_1,c_2\right)$ using the resulting solution of the Cauchy problem with the sought solutions $\mathrm{\Psi }\left(c_1,c_2\right)$. In this procedure, a four-dimensional matrix $A_M\left(c_1,c_2\right)$ is reduced to a two-dimensional one by the relevant renouncement of its elements.

3.2 Defining the area of nonzero solutions

Due to the peculiarity of the problem statement according to desired power DP $P\left(x_1,x_2\right)$, Eq. (8) has zero solution at arbitrary values of the parameters $c_1,c_2,\alpha$. From an engineering point of view, this is a significant drawback, but for some desired DPs $P\left(x_1,x_2\right)$ it is possible to fix an area of parameters $c_1,c_2,\alpha$ at which a nonzero solution exists. At the small $c_1$ and $c_2$, the kernel (9) is given approximately in the form

$$K\left(c_1,c_2,x_1,x_1^{\text{'}},x_2,x_2^{\text{'}}\right)\approx \frac{M_2{N}_2{c}_1{c}_2}{{\pi }^2}$$E25

Assuming that $f\left(x_1,x_2\right)$ is constant, the integral Eq. (8) can be rewritten as (usually for small $c_1$ and $c_{2}f\left(x_1,x_2\right)\approx \text{const}$).

$$\frac{{\pi }^2\alpha }{2M_2{N}_2{c}_1{c}_2}=\underset{-1}{\overset{1}{\int }}\underset{-1}{\overset{1}{\int }}P\left(x_1,x_2\right){\mathit{dx}}_1{\mathit{dx}}_2-4{\left|f(x_1,x_2)\right|}^2$$E26

The area of integration $\mathrm{\Omega }$ in Eq. (8) is reduced in the last formula to the area $\left[-1,1\right]\times \left[-1,1\right]$ because of definition of both the arguments $x_1,x_2$ and parameters $c_1,c_2$.

Taking into account that ${\left|f\left(x_1,x_2\right)\right|}^2$ is positive, we get the following relationship between the function $P\left(x_1,x_2\right)$ and the parameters $c_1$, $c_2$, and $\alpha$:

$$\underset{-1}{\overset{1}{\int }}\underset{-1}{\overset{1}{\int }}P\left(x_1,x_2\right){\mathit{dx}}_1{\mathit{dx}}_2-\frac{{\pi }^2\alpha }{2M_2{N}_2{c}_1{c}_2}>0$$E27

Finally, considering the case $P\left(x_1,x_2\right)\equiv 1$, we obtain:

$$c_1{c}_2>\frac{{\pi }^2\alpha }{8M_2{N}_2}$$E28

In fact; inequality (28) determines the area of parameters $c_1,c_2,\alpha$, where nonzero solutions exist. In Figure 2, the dependence curves $c_2=c_2\left(c_1\right)$ for three different values $M_2$ and $N_2$ are shown. The results are given for array with the number of elements $M_2=N_2=3$ (curve 1), $M_2=N_2=5$ (curve 2) and $M_2=N_2=11$ (curve 3). The area of values $c_1$ and $c_2$, where the existence of zero solutions is possible, according to the estimate (28) is located below and to the left of the presented curves. As can be seen, the area of zero values decreases significantly with increasing $N_2$and $M_2$. The obtained results testify that the zero solutions of Eq. (8) for a given constant power DP can exist either at a small value $c_1{c}_2$ corresponding to low frequencies (at a given size of array), or at the values of $c_1$ that significantly exceeding $c_2$ and vice versa. The last case corresponds to arrays with a large difference in distance of elements along the coordinate axes. Such arrays are usually rarely used in practice.

Figure 2.

The curves $c_2=c_2\left(c_1\right)$ at the different $M_2$ and $N_2$_.

3.3 Determination of bifurcation lines

3.3.1 The case of rectangular array

The finding of bifurcation lines of the nonlinear Eq. (8) was performed for the array containing $N_2\ast M_2=11\ast 11=121$ radiators for the desired power DP $P\left(x_1,x_2\right)=1$ at ${\mathrm{\Lambda }}_c=\left\{\left(c_1,c_2\right),0<c_1,c_2\le 2\right\}$ for the different values of the parameter $\alpha$ in (3).

The search for bifurcation lines can be performed directly by investigating the properties of the determinant (20) as a function of the parameters $c_1$ and $c_2$. In addition, the function (20) depends on the parameter $\alpha$; so the set of its eigenvalues also depends on this parameter, i.e. the set of spectral curves that separate the areas of zero and nonzero solutions.

The behavior of the corresponding curves when changing the parameter $\alpha$ is shown in Figure 3. The behavior of the determinant (24) depending on the parameters $c_1$and $c_2$ at $\alpha =0.5$ is given in Figure 3a; and in Figure 3b–d, the intersection of this function with a plane $\mathrm{\Psi }\left(c_1,\gamma c_1\right)=0$ is illustrated at the different $\alpha$. This results in a set of curves that correspond to a set of spectral lines separating the area of zero and nonzero solutions. At a fixed size of array, the area where zero solutions can exist expands if the parameter $\alpha$ increases, this area is located below the left of the first curve.

Figure 3.

The spectral curves of Eq. (19) at the different $\alpha$.

The curves marked by number 1 correspond to the solutions with constant (zero or even) phase DP; curves with number 2 correspond to the solutions with phase DP that is even with respect to the ${\mathit{Ox}}_1$ axis, and odd with respect to the ${\mathit{Ox}}_2$ axis, and curves numbered by 3 correspond to the solutions with a phase DP odd with respect to two coordinate axes. The proposed procedure is quite approximate, it does not allow to separate the curves that correspond to different types of solutions and thus identify the areas where there is a nonzero solution for the synthesized power DP with the specified phase property.

The method of implicit function proposed in [23] and developed for plane array in [30] is devoid of this drawback.

At the first step of this method, a series of one-dimensional eigenvalue problems is solved, by this the different values of parameter $\gamma$ are prescribed by the relation $c_2=\gamma c_1$and a one-dimensional problem is solved with respect to $c_1$. In Figures 4 and 5, the first four eigenvaluesof the problem at $\gamma =1.0$ and $\gamma =0.2$ are shown. The values $\left(c_1^{\left(i\right)},c_2^{\left(i\right)}=c_1^{\left(i\right)}\right)$, $i=1,2,3,4$, at which ${\lambda }_i=1$, are the bifurcation points in the plane $c_1{\mathit{Oc}}_2$. By this, the set of points $\left(c_1,c_2\right)$at which the eigenvalue ${\lambda }^{\left(i\right)}=1$ is determined approximately from the graphical data.

Figure 4.

The first eigenvalues at the ray $c_1=c_2$, $\alpha =0.5$.

Figure 5.

The first eigenvalues at the ray $c_2=0.2c_1$, $\alpha =0.5$.

The next step is to refine the values $\left(c_1^{\left(i\right)},c_2^{\left(i\right)}\right)$ by solving the transcendental Eq. (20), and the point $\left(c_1^{\left(i\right)},c_2^{\left(i\right)}\right)$, which is considered as the initial approximation.

In the final step, the bifurcation curve in the plane $\left(c_1,c_2\right)$ is determined by solving Eqs. (21) and (22), after specification of the values $\left(c_1^{\left(i\right)},c_2^{\left(i\right)}\right)$. In Figure 6, the bifurcation curves $\left(c_1^{\left(1\right)},c_2^{\left(1\right)}\right)$–$\left(c_1^{\left(4\right)},c_2^{\left(4\right)}\right)$ that correspond to the first four eigenvalues are shown. The curve with number 1 corresponds to the solution with the zero (even) phase of the created DP. This curve corresponds to that is marked by 1 in Figure 3b.

Figure 6.

The bifurcation curves corresponding to set $\left(c_1^{\left(i\right)},c_2^{\left(i\right)}\right),i=1,\dots ,4$, $N_2=M_2=11$.

There are no nonzero solutions with such a phase property for the values $c_1$ and $c_2$ above and to the right of this curve. Curves 2 and 2′ correspond to solutions in which the phase DP is symmetric about one axis and asymmetric about the other axis (obviously, for a plane array, there are two such curves and they are antisymmetrical). Curve 2 corresponds that is marked by 2 in Figure 3b. The curve with number 3 corresponds to a solution with a phase DP antisymmetrical (odd) with respect to both the axes. The location of the areas of zero and non-zero solutions is the same as in Figure 3b. It should be noted that the problem of refining the roots of Eq. (20) is the most time-consuming in computational relation because refining the roots of this equation requires a series of computational experiments with different values $\left(c_{1,0}^{\left(i\right)},c_{2,0}^{\left(i\right)}\right)$ of initial parameters close to approximate values.

3.3.2 The case of hexagonal array

Firstly, we consider the procedure of determination of bifurcation curves by finding zero lines of determinant (24). The results, similar to those are presented in Figure 3a and b for the rectangular array, are shown in Figures 7 and 8. One can see that the behavior of function $\mathrm{\Psi }\left(c_1,c_2\right)$ is more complex than in the case of rectangular array. The obtained graphs testify that the solutions with other different behavior of phase $arg\left(f\left(x_1,x_2\right)\right)$ of the DP appear additionally. One such solution is marked by number 4. Other solutions appear when parameters $c_1$ and $c_2$ increase at the fixed $\alpha$.

Figure 7.

The surface of determinant (24) values for the hexagonal array, $\alpha =0.5$.

Figure 8.

Zero lines of determinant (24) for the hexagonal array, $\alpha =0.5$.

Search of the bifurcation curves is carried out similarly to the case of rectangular array. The numerical results are presented for the array with $N_{\mathit{tot}}=61$ elements for the desired power DP $N_0\left(x_1,x_2\right)=1$ at ${\mathrm{\Lambda }}_c=\left\{\left(c_1,c_2\right),0<c_1,c_2\le 2.0\right\}$ for the different values of $\alpha$ in the functional (3). At the first step, the one-dimensional eigenvalue problems were solved at the different values of parameter $\gamma$. In Figure 9, the first four eigenvalues are shown at $\gamma =1.0$, and in Figure 10, they are shown at $\gamma =0.2$. Similar to the case of rectangular array, the points, in which ${\lambda }_i=1$ are moved to right and the distance between them increases at $\gamma =0.2$. The values $\left(c_1^{\left(i\right)},c_2^{\left(i\right)}=\gamma c_1^{\left(i\right)}\right)$, where $i=1,2,3,4$, are the bifurcation points in the plane $\left(c_1,c_2\right)$. The points $\left(c_1^{\left(i\right)},c_2^{\left(i\right)}\right)$, for which the eigenvalues ${\lambda }^{\left(i\right)}=1$ are determined approximately in this step.

Figure 9.

The first eigenvalues for the ray $c_1=c_2$ at $\alpha =0.5$.

Figure 10.

The first eigenvalues for the ray $c_2=0.2c_1$, $\alpha =0.5$.

The specification of values $\left(c_1^{\left(i\right)},c_2^{\left(i\right)}\right)$ by solving Eq. (20) is carried out in the next step, and the points $\left(c_1^{\left(i\right)},c_2^{\left(i\right)}\right)$ of the graph data from the Figures 9 and 10 are used as initial approximations. The usual numerical half-division method is used for this goal.

The bifurcation points $\left(c_1^{\left(i\right)},c_2^{\left(i\right)}\right)$, $i=1,2,3,4$ for the first four eigenvalues in the rays, $c_2=\gamma c_1$ are shown in Figure 11. The respective curves of bifurcations, which are obtained by solving Eqs. (21) and (22), are shown in Figure 12. As in the case of a rectangular array, to obtain the necessary data, we should carry out precise computations.

Figure 11.

The bifurcation points at the rays $c_2=\gamma c_1$.

Figure 12.

The bifurcation curves in the plane $\left(c_1,c_2\right)$.

4. The engineering applications

The results presented in this Section demonstrate how the knowledge about the point of bifurcation obtained as the solutions of the nonlinear eigenvalue problems allows us to understand better the process of bifurcation and how to get the solutions, which are the most optimal in sense of the used criterion of optimization.

4.1 The method of successive approximations

The properties of solutions to Eq. (8) obtained by using the method of successive approximation are related directly with the properties of phase characteristic of the eigenfunctions, which are determined at solving the eigenvalue problem. Prescribing the initial approximation $f_0$ for the iterative process for solving Eq. (8) with the specified property of the phase $arg{f}_0$, we could receive the solution of Eq. (8) with the same phase property in the wide range of characteristic parameters $c_1$ and $c_2$. This is important for the engineering design of arrays having the fixed phase characteristics of radiation in the defined range of frequencies.

The method of successive approximations

$$f_{n+1}-\beta f_n+\left(1-\beta \right)B\left(f_n\right)=0,n=0,1,2,\dots$$E29

is used for solving Eq. (8) with a set of specific physical parameters of array. In the last formula,

$$B\left(f\right)=\frac{2}{\alpha }\left[{\mathit{AA}}^{\ast }\left[\left(P\cdot f\right)-\left({\left|f\right|}^2\cdot f\right)\right]\right]$$E30

Parameter $\beta \in \left[0,1\right]$ in (29) is used to accelerate the convergence of iterative process. To substantiate the condition of convergence of the iterative process (29), we apply Theorem 2.6.2 [26] (p. 133), which states that the operator ${\mathit{AA}}^{\ast }$ be contraction one. This requirement is met when the inequality

$$\alpha >2\Vert {\mathit{AA}}^{\ast }\left[f\left(N_0-{\left|f\right|}^2\right)\right]\Vert$$E31

met. The results of numerical calculations show that condition (31) is overestimated and for some values of the problem parameters the iterative process (29) converges for values $\alpha$that do not satisfy the estimate (31).

4.2 The case of rectangular array

In Figure 13, the dependence of the convergence of the iterative process (29) on the value of parameter $\alpha$ for a desired power DP $P\left(x_1,x_2\right)=1$ at the fixed values $\beta =0.1$, $c_1=c_2=2.0$, the number of radiators $M_2\cdot N_2=11\cdot 11=121$ is shown. The required accuracy $\epsilon ={10}^{-3}$.

Figure 13.

The character of convergence of iterative process (29) at the different $\alpha$, $M_2=N_2=11$.

The results of solving the optimizing problem for this desired power DP at $\alpha =0.5$ are shown in Figure 14. The approximation quality to a desired DP $P$ significantly depends on the parameters $c_1$ and $c_2$ at both fixed $N_2$ and $M_2$. The mean-square deviation (MSD) (the first term in (3)) is equal to 0.0847 for $c_1=c_2=1.0$, and it is equal to 0.0075 for $c_1=c_2=3.14$.

Figure 14.

The created power RP $\mid f\mid$ (a) and optimal distribution of currents $\mid I_{\mathit{nm}}\mid$ (b) at $c_1=c_2=1.0$.

The synthesized DP $\mid f\mid$for larger $c_1,c_2$ has not only a more optimal mean-square approximation, but it is also closer to the shape of the desired DP $P$. The optimal amplitudes $\mid I_{\mathit{nm}}\mid$ of currents in the array’s elements are close to constant at such parameters $c_1$ and $c_2$.

When solving the optimizing problem for desired power DP of a more complex form, the quality of the approximation significantly depends on both the parameter $\alpha$ and the type of initial approximation for the phase of a given DP. The results are shown for the desired power DP

$$P\left(x_1,x_2\right)=\mid sin\left(\pi x_1\right)\cdot sin\left(\pi x_2\right)\mid ,-1\le x_1\le 1,-1\le x_2\le 1$$E32

at $c_1=c_2=3.14$ and $\alpha =0.2$ in Figure 15. Despite the fact that the shape of desired DP $P\left(x_1,x_2\right)$ is more complex that in the previous example, decrease of $\alpha$ (from 0.5 to 0.2) and simultaneous increase of $c_1$ and $c_2$ (from 1.0 to 3.14) allows us to get the amplitude $\mid f\left(x_1,x_2\right)\mid$ of created DP, which is very close to the $P\left(x_1,x_2\right)$. The optimal distribution of currents’ amplitudes $\mid I_{\mathit{nm}}\mid$ (Figure 15b) approaches the shape of the created DP.

Figure 15.

The created power RP ${\left|f\right|}^2$ at $c_1=c_2=3.14$.

We have used an additional optimization parameter $\beta$ in Eq. (29), which, as shown by the results of numerical calculations, accelerates the convergence of iterative process significantly. In Figure 16, the results of the study of the influence of this parameter on the rate of convergence at a fixed value of the parameter $\alpha =0.5$ are shown. The results are given for the desired power DP $P\left(x_1,x_2\right)=1$ at $c{}_1=c_2=2.0$. In order to achieve the accuracy ${10}^{-3}$ of calculations, one needs 157 iterations at $\beta =0.01$. If parameter $\beta$ increases to a certain value, the number of iterations decreases significantly, so at $\beta =0.05$, $\beta =0.10$, $\beta =0.15$, $\beta =0.20$, and $\beta =0.25$ one requires 62 iterations, 37 iterations, 28 iterations, 20 iterations, and 16 iterations, respectively. At the subsequent increase, the number of required iterations begins to increase and already at $\beta =0.30$ the iterative process begins to diverge. Numerical calculations show that the limit valueof $\beta$, at which the iterative process (29) begins to diverge, significantly depends on the value of the parameter $\alpha$. So, if this parameter decreases, the threshold value of $\beta$ increases. The dependence of the convergence on the array’s parameters ($c_1,c_2,M_2,N_2,d_1,d_2$) is not so significant.

Figure 16.

The convergence of iterative process (29) for the different $\beta$.

The values of the functional (3) for the created DP with different phases are shown in Figure 17. The solid curve corresponds to the phase DP even with respect to two axes, the dotted curve corresponds to phase DP even with respect to one axis and odd with respect to the other, the dashed curve corresponds to the phase DP odd with respect to both the axes.

Figure 17.

The values of functional (3) versus the phase of created DP.

One can see that the values of functional at the fixed $c$ (frequency) significantly depend on the phase of the created DP. The value $\sigma =0.2$ is achieved for the “even-even” solution at $c=0.79$, for the “even-odd” solution at $c=0.71$ and for the “odd-odd” solution at $c=0.625$. That is, within the used criterion, the latter type of solution is 21% better than the first one. From this fact, it follows that at a fixed distance between the radiators for the desired DP $P\left(x_1,x_2\right)=1$, the number of array’s elements can be reduced by 21% with the same value of MSD. A similar situation is observed for the characteristics of DP at $\sigma =0.1$, i.e. “odd-odd” solution is better on 19.4% than “even-even”.

4.3 The case of hexagonal array

The results of solution of the optimization problem for two given power DPs $P_1\left(x_1,x_2\right)\equiv 1$ and

$$P_2\left(x_1,x_2\right)=\left\{\begin{array}{c}2\sqrt{x_1^2+x_2^2}\sqrt{1-x_1^2-x_2^2},x_1^2+x_2^2\le 1,\\ 0,x_1^2+x_2^2>1,\end{array}\right.$$E33

in the form of body of rotation are shown in Figures 18 and 19 at $\alpha =0.5$.

Figure 18.

The amplitude of created DP ${\left|f\right|}^2$ for $P_1$ at $c_1=2.0,c_2=2.236$.

Figure 19.

The amplitude ${\left|f\right|}^2$ of created DP for $P_2$ at $c_1=2.0,c_2=2.236$.

As previously, the optimization problem consists of solving Eq. (8) by the method of successive approximation (29). The MSD (the value of the first term in (3)) for the first desired DP is equal to $0.3774$, and it is equal to $0.2218$ for the second desired DP.

Similar to the case of rectangular array, the approximation quality to the desired DP $P$ depends on both the parameters $c_1$, $c_2$, and $\alpha$. The characteristic of MSD of DPs for $\alpha$ at the different $c_1$ on the ray $c_2=1.118c_1$ is shown in Figures 20 and 21. The chosen relation between $c_1$ and $c_2$ provides the regularity of the array’s geometry, and as the numerical computations have shown, gives the ability to get the close characteristics of radiation in the planes $x_1$and $x_2$.

Figure 20.

The MSD versus the parameter $\alpha$ for the DP $P_1$.

Figure 21.

The MSD versus parameter $\alpha$ for DP $P_2$.

The largest MSD for the $P_1$ is achieved at $\alpha =1.0$ for $c_1=0.5$, and it is equal to $1.96443$, it diminishes almost linearly if parameter $\alpha$ decreases. The largest MSD for DP $P_2$ is equal to $1.43685$. One should note that the value of MSD diminishes if $\alpha$ decreases, but the norm ${‖\mathbf{I}‖}_{H_{\mathbf{I}}}$ of current growth that is unacceptable from the engineering point of view.

The approximation quality to the desired DP depends also on the type of initial data, which are prescribed for the iterative process (29). The dependence of the values of functional (3) on the parity of phase of the initial approximation $f_0$ for the DP $P_1$ is shown in Figure 22. The results are shown for four types of initial approximation: odd with respect to both the axes (dotted curve), even with respect to the ${\mathit{Ox}}_1$ axis and odd with respect to the ${\mathit{Ox}}_2$ axis (dashed curve), odd with respect to the ${\mathit{Ox}}_1$ axis and even with respect to the ${\mathit{Ox}}_2$ axis (dash-dot curve), and even with respect to both the axis (solid curve). The initial approximation $f_0$, corresponding to the even phase with respect to both the axis, is optimal for this DP, moreover the values of ${\sigma }_{\alpha }$ for the small values of parameters $c_1$ and $c_2$ differ significantly, but starting from $c_1=0.8$ this difference does not exceed 10%.

Figure 22.

The values of ${\sigma }_{\alpha }$ versus the initial approximation of initial approximation for the iterative process (29).

The dependence of convergence of the iterative process (29) on the parameter $\alpha$ at the fixed $\beta =0.1$ is shown in Figure 23. As in the case of a rectangular array, the iterative process converges most slowly at $\alpha =1.5$, one needs 67 iterations to achieve the accuracy that is equal to ${10}^{-3}$. The number of iterations decreases if $\alpha$ diminishes. For example, one needs 30 iterations to achieve the same accuracy at $\alpha =0.2$.

Figure 23.

The convergence of iterative process (29) versus $\alpha$.

The dependence of convergence on the parameter $\beta$ is studied too. The necessary number of iterations that needs to achieve the accuracy ${10}^{-3}$at $\beta =0.01$, $\beta =0.05$, $\beta =0.10$, $\beta =0.15$, $\beta =0.20$, and $\beta =0.30$ (curves 1–6 respectively) is shown in Figure 24. It is substantiated that the iterative process converges most slowly at $\beta =0.01$, it is necessary 152 iterations to achieve the prescribed accuracy. The most optimal among the considered $\beta$ is $\beta =0.20$, one needs 20 iterations only to achieve the above accuracy. The iterative process becomes slow at the subsequent growth of $\beta$. For example, one needs 37 iterations to achieve this accuracy; the iterative process (29) becomes convergent at $\beta >0.40$. This testifies that in the process of computations one should to limit by non-large values of $\beta$ ($\beta \le 0.20$) that guarantees the convergence and considerably grows its speed on the contrast with small $\beta$ ($\beta \le 0.01$).

Figure 24.

The convergence of iterative process (29) versus number of iteration, $\alpha =0.5$.

More information about the problem under investigation one can find in [30, 31, 32, 33, 34].

5. Conclusions

The problem of finding the solutions to the nonlinear integral equations and their properties is reduced to nonlinear two-dimensional eigenvalue problems that lead to the subsequent application of an implicit function method for solving the Cauchy problem for the linear differential equation. The area of non-zero solutions to the above equations is determined by involving the solving transcendental equation, which is got by equating to zero of determinant related to the eigenvalue problem. The results of solving the nonlinear eigenvalue problems are applied subsequently for specification of the bifurcation points and obtaining the bifurcation curves. The approach does not depend on the form of operator determining the radiation properties of physical system (plane rectangular and hexagonal arrays). The obtained results are the constructive basis on which a series of practical engineering problems of optimization was solved numerically.

Conflict of interest

The author declares no conflict of interest.