Numerical solution of the problem of recovering parameters for the limit model of an uninsulated vacuum plane diode

The model of a vacuum diode under the influence of a strong external magnetic field is considered. The uninsulated variant, when a part of electrons emitted from the cathode reaches the anode, is investigated. The model is described by a singular boundary value problem for a system of ordinary differential equations. The sensitivity of the problem solution to the change of input parameters is investigated. A coordinate descent method to restore parameters of the boundary value problem is implemented numerically.


Introduction
The modern statement of the problem of magnetic isolation was formulated by physicists in the late 80s of the last century. The magnetic insulation effect is that under the influence of a strong external magnetic field the electrons emitted from the cathode do not reach the anode and turn back to the cathode. This effect creates an electronic layer outside of which the electromagnetic field is equal to zero. Here two main modes are possible: the first one is when electrons reach the anode ("uninsulated diode"), and the second one is when electrons turn back to the cathode ("insulated diode"). In addition, there is an intermediate mode where some electrons can reach the anode while others cannot.
The greatest interest is in modeling the complete problem, which is a general relativistic Vlasov-Maxwell equations [Abdallah et al., 1998]. This general model is still waiting for its researchers.
In [Abdallah et al., 1998;Sinitsyn, 2001;Semenov et al., 2010] the limit model of magnetic insulation described by the boundary value problem for the system of two nonlinear second order ordinary differential equations is considered. A specific feature of this system is its singularity. In [Abdallah et al., 1998;Sinitsyn, 2001;Semenov et al., 2010;Kosov et al., 2012] a number of analytical results were obtained and numerical experiments were carried out to solve boundary value problems with different parameters. Nevertheless, the problem of magnetic insulation simulation turned out to be rather complicated.
In this paper, a switched algorithm is used for numerical solution of the singular boundary value problem. A series of numerical experiments was carried out to study the sensitivity of the problem. To solve the inverse problem of parameter recovery, the coordinate descent method was proposed and implemented. This method proved to be quite effective for the problem under consideration.

Mathematical Model
The magnetic insulation effect is that under the influence of a strong external magnetic field the electrons emitted from the cathode do not reach the anode and turn back to the cathode. This effect creates an electronic layer outside of which the electromagnetic field is equal to zero. Here two main modes are possible: the first one is when electrons reach the anode ("uninsulated diode"), and the second one is when electrons turn back to the cathode ("insulated diode"). In addition, there is an intermediate mode where some electrons can reach the anode while others cannot. However, the main system, called the limit problem, describing the uninsulated diode mode is given by the following equations: where x ∈ [0, 1] is a independent variable, which determines a relative distance from the cathode (x = 0) to the anode (x = 1); ϕ(x) is the function which describes the change in the electric field potential when moving from the cathode to the anode; the function a(x) denotes the magnetic potential; j x describes the electrical current density via diode. The equation (1) describes the electric and magnetic fields inside the diode, and its solution has to meet the initial and boundary conditions characterizing the natural physics of the process describing the mode of uninsulated diode. So, on the segment x = [0, 1] we have: the initial conditions: the boundary conditions: β is a parameter in the initial conditions (2) which characterizes the magnetic field near the cathode.
The study of the system (1) -(3) is connected with its numerical solutions. The problem simulation has a basic character, as the analysis of the behavior of integral curves can and should be correlated with the known physical effects of the behaviour of vacuum diodes. To overcome these difficulties, the authors have done their work to identify the simplest admissible numerical methods and the development of a special method for finding the parameters of the model by given boundary conditions.
Changing the variables we get the following boundary value problem y 1 (0) = 0, y 2 (0) = 0, y 3 (0) = 0, y 4 (0) = β, Let us define the first statement of the inverse boundary value problem for the limit model (5) -(7), describing the mode of uninsulated vacuum diode: Statement 1: For the given parameters ϕ L and a L it is necessary to restore j x and β.
Note that for the system (5) -(7) we can also define a number of other inverse problems: Statement 2: For the given parameters a L and β it is necessary to restore j x and ϕ L . Statement 3: For the given parameters ϕ L and β it is necessary to restore j x and a L . Statement 4: For the given parameters a L and j x it is necessary to restore β and ϕ L . Statement 5: For the given parameters ϕ L and j x it is necessary to restore β and a L .
In the authors' opinion, simulation of all statements of the inverse problems for the system (5) -(7) will allow to consider all aspects of the uninsulated diode model, which can be correlated with its physical properties. Within the framework of the current work, we will consider only the first statement, as it corresponds to the principal physical evaluation of diode operation parameters in the steady-state (limit) mode relative to the initial given potentials of magnetic and electric fields.

The Sensitivity of the Numerical solution of a Quasi-singular Problem
Note that at least some complete simulation of the limit case as a system of ordinary differential equations (5) - (7) is not available at this moment, except for the initial work of the authors [Varin, 2013]. The situation has developed, in particular, because this problem is described by an autonomous system of ODE of the second order, which has a discontinuity in the starting point.
For given j x and β the simulation problem is determined by the initial problem, which requires two subproblems to be solved: 1. to choose of the numerical integration method; 2. to determine the solvability of the solution with respect to the singularity of initial conditions.
As an integration algorithm it is proposed to use the LSODA switched algorithm developed by Petzold [Petzold, 1980] and included in scipy.integrate package. This algorithm automatically chooses between the soft method of predictor-corrector Adams method and the rigid Gear method [Gear, 1971;Shampine, 1979] with the inverse differential formula (Backward Differentiation Formula (BDF) method) directly during the solution. At the initial stage of the solution, which is not usually rigid, LSODA uses more efficient Adams methods. If the presence of stiffness is detected, then the automatic transition to the method of Gear solution for stiff systems is carried out. Given the stiffness of the (5) model, this property of the algorithm is mandatory. It also be noted that in the calculations by the Adams method, the error made at any step does not tend to exponential growth.
To overcome the problem of singularity at the numerical solution of the boundary value problem it is proposed to choose the values y 1 (0) and y 3 (0) as non-zero, but rather small, having connected the parameter of smallness with the value of the integration step. A detailed study of the system equations (5) shows that supposing ϕ (2) (0) ≡ a (2) (0) = const (the value of the constant is arbitrary) preserves the integrity of the system. It follows from this that all explicit numerical methods of ODE solution are not applicable at the starting point. To solve the problem within the framework of the numerical solution it is proposed to take these parameters small enough and pass to the quasi-singular problem concerning parameters y 1 (0) = ε(0) 1 and y 3 (0) = ε(0) 2 , in which parameters ε(0) 1 and ε(0) 2 are set with a sufficient order of smallness (relative to the integration step) about 0.
For numerical analysis of sensitivity of the method of integration of the quasi-singular initial problem concerning initial conditions y 1 (0) = ε(0) 1 and y 3 (0) = ε(0) 2 we will consider the problem at the following parameters: Figs. 1 and 2 show the results of the numerical analysis of the initial problem solution sensitivity for (5) -(7) at the decrease of the initial parameters y 1 (0) and y 3 (0) at the system integration (5) -(7) with the interval division into 50,000 points. The results of modeling showed the sensitivity of the numerical solution of the quasi-singular initial problem from the initial condition of y 1 (0) (within the 4th order) and weak sensitivity from parameter y 3 (0). Thus for both dependences presented on figures have the tendency to convergence at approaching to zero of initial condition y 1 (0).
Tables 1 and 2, Figs. 3 and 4 show the results of the numerical experiment on the sensitivity of the solution of the initial problem for (5) -(7) with the reduction of the initial parameters y 1 (0) and y 3 (0), as well as the reduction of the integration step.
As can be seen from the analysis of the data presented in Tables 1 and 2, Figs. 3 and 4, the results of our numerical solution for the selected integration method are insensitive to the value of the numerical integration step, which indicates the stability of the method in relation to the choice of the size of the integration grid. In general, we can conclude that the chosen method of integration from the position of the solution of the quasi-singular problem allows taking into account the convergence of the solution when approaching zero initial parameters y 1 (0) and y 3 (0) and the invalidity of the method to the integration step.    of the quasi-singular initial problem of the ODE system was carried out at (5) -(7) at iterative change of initial conditions of j x , β. As a result, we obtained the dependences that determine the domain of permissible values of boundary conditions in the subspace (ϕ L , a L ), presented in Fig. 4. The results of the numerical solution of the quasisingular initial problem allow us to formulate the following conclusions.
1. Restored characteristics of the boundary value problem in the first statement determine the smooth functional in the space of (ϕ L , a L ) 2. Parameter β = y 4 (0) defines the point movement for some "parabolic" curve (ϕ L , a L ) k = G k (β)| jx=j 0 x +k·∆j=const 3. When increasing j x = j 0 x + k · ∆j "minimum" the characteristics of G k (β) with respect to the axis of abcissus ( a L = 0) shifts to the right. 4. All "parabolic" curves G k are restricted by asymptotes, i.e. the solution area at the right edge is limited. So, for example, the solution point ϕ L = 1 and

Solution of the Inverse Problem by the Method of Coordinate Descent
Consider Statement 1 of the inverse boundary value problem Taking into account the chosen method of numerical integration and the properties of the method's sensitivity when solving the initial quasi-singular problem (5) -(7), let us define the algorithm of the coordinate descent in the domain (ϕ L , a L ) 1. Setting the boundary conditions s L (ϕ L , a L ) for solution of the quasi-singular problem (5) - (7) and numerical error of the solution ε L regarding the Euclidean norm. 2. Selection of arbitrary initial conditions j x  x , ∆β [0] ). 6. Assigning the current iteration of initial parameters and starting the external coordinate descent cycle 7. Checking the cycle continuation criterion for the kth epoch of the coordinate descent by the criterion When the condition is fulfilled, the transition to the internal cycle (item 8) (b) If the condition is not valid, the search is completed (item 12).
8. Internal cycle: a component local coordinate motion relative to the given coordinate step at the current epoch k: 9. Solution of the quasi-singular problem relative to the given coordinate step (5) - (7) for determination of the next coordinate point s j+1 ).

Checking the condition of the internal cycle continuation
11. (a) If the internal cycle condition is not met for all coordinate directions of motion, the coordinate step for the next iteration is reduced and it's reaching the external cycle level (item 7). (b) if the internal loop condition is met, move to the next iteration step (item 8). 12. Output of the result j x [k] , β [k] .
The above algorithm, taking into account the limited area of the solution in the domain (ϕ L , a L ) should also contain the stop conditions in case of its cycling, for example, when choosing such boundary conditions that cannot be achieved due to the space of the solution of the boundary functional (see Fig. 4). Therefore, to stop the algorithm, we can offer a maximum limit on the number of iterations of the external loop. So, for example, if the initial step is ∆z [0] = (∆j [0] x , ∆β [0] ) = (1.0, 1.0), then already at the epoch with the index k = 50 we will have a step of approximation ∆z [50] = (1, 776e−15, 1, 776e−15), which can be considered a good approximation.

Results of the Inverse Problem Numerical Solution
Here are examples of numerical solution trajectories for (ϕ(x) and a(x) under different boundary conditions. For computational experiments let us take test points: The experiment 1: ϕ L = 1, a L = 1. The experiment 2: ϕ L = 8, a L = 3. The experiment 3: ϕ L = 0.3, a L = 0.8.
The specified test points lie in the area of boundary functional solutions (4). The point ϕ L = 0.3, a L = 0.8 approaches the upper asymptotic boundary. Experiment 1. Let us set ϕ L = 1, a L = 1 belonging to the solution area at the right end. We need to define j x , β. The starting point for the algorithm is s  (1; 1), the numerical error of the solution is ε L = 1e − 9 Results of the solution of the boundary value problem by the coordinate descent algorithm are given below.

Number of epochs k = 30
Euclidean norm error ε L = 6.40e − 10 Restored characteristics: Experiment 2. Let us set ϕ L = 8, a L = 3, belonging to the solution area at the right end. We need to define j x , β. The starting point for the algorithm is s  (1; 1), the numerical error of the solution is ε L = 1e − 8 Results of the so-lution of the boundary value problem by the coordinate descent algorithm are given below.

Number of epochs k = 28
Euclidean norm error ε L = 6.58e − 09 Restored characteristics: Graphs of the numerical solution are presented in Fig.  7, and the process of the coordinate descent for the numerical solution is illustrated in Fig. 8: Experiment 3. Let us set ϕ L = 0.3, a L = 0.8, belonging to the solution area at the right end (these parameters characterize the limit solution of the boundary value problem in relation to the domain of admissible values, as they are close to the asymptotic solution). We need to define j x , β. The starting point for the algorithm is s Number of epochs k = 28 Euclidean norm error ε L = 5.52e − 09 Restored characteristics: During the testing of the algorithm it was determined that the proposed algorithm showed good convergence within the boundary functional area (Fig. 4). At the same time, there were problems with convergence at the boundaries of the boundary functional (limiting asymptotes of the G k functional) and the proposed algorithm of coordinate search hovered relatively small values of the step within the proposed criterion of the Euclidean norm. To understand this aspect we considered functional dependencies ϕ L (j x ) and a L (j x ), which are obtained by fixing parameter y 4 (0) = β = const (Fig. 9).
Dependencies analysis ϕ L (j x ) and a L (j x ) has shown that the following three areas can be separated in the solution space of the boundary value problem Figure 11. The internal area of solving the boundary value problem at growth β: a) 2d projection in a space ϕ L , a L ; b) 3d projection in a coordinate space ϕ L , a L , j x Left asymptotic region of the boundary value problem solution relative to the points obtained for the points with coordinates at in which nonlinear harmonic fluctuations in characteristics are observed. The right stable area of the boundary value problem solution relative to the points obtained at that shows smooth, increasing characteristics. Transition area between the left and right regions relative to the received points argmax [min (ϕ L (j x )) , min (a L (j x ))] Among the presented areas of interest are the asymptotic and internal areas of solution. The internal area of the solution characterizes the parameters of the work of "uninsulated diode" in the steady-state mode relative to the initial given potentials of magnetic and electric fields. Within the limits of the specified area it is possible to draw a conclusion on a stable mode of work of a diode and the unique decision of inverse boundary value problem within the limits of the first statement.
The key features of the solution are the following.
The unique solution inside the boundary functional area G = f (ϕ L , a L ) exists. The solution at the boundary of f the boundary functional area G k = f (ϕ L , a L ) is "singular", as there is uncertainty of the solution concerning the boundary functional asymptote. It happens because of the fact that the edge characteristics with different current j x fall on the asymptotes (5 3d graph.). In the author's opinion, the approach of the fast channel parameter optimization technique [Altsybeyev et al., 2018] may also prove to be effective in the problem under consideration.

Conclusion
Thus, numerical solution of the considered inverse problem for the "uninsulated diode" demonstrates an interesting boundary effect associated with the uncertainty of the boundary solution. This effect due to the harmonic instability of boundary parameters ϕ L and a L on the asymptotes of the solution with respect to the area of acceptable solutions of the functional G k can give one and the same solution of the inverse boundary value problem at different boundary parameters of the current density through the diode j x and the potential of the own magnetic field near the cathode β (Fig. 11). In this case the inverse problem is in restoration of set of decisions and research of borders of the given set. If to recollect principles of functioning of a vacuum diode and remarks about its own noises the considered mathematical model in the given area shows that small change of a magnetic field (for example external effects) can considerably change characteristics of a diode that characterizes approximations of characteristics of a diode to magnetic-isolating. Also, the limiting effect is manifested in the coherent fluctuation of the characteristics of the dependences of the numerical solution relative to the first derivative ofφ (x),ȧ(x) (Fig. 12) These similarities require the introduction of special methods of singular bifurcation analysis, which will be discussed in subsequent articles in the framework of the solution of other formulations of the boundary value problem.