Software for the numerical solution of first-order partial differential equations

Full Text

Introduction Partial differential equations (PDE) of the first order arise in many applied problems. The specificity of the problems in optics and optoelectronics is that the coefficients of these equations are not defined by a single analytical expression in the entire considered domain. Thus, for example, the eikonal equation is a nonlinear first-order PDE, one of whose coefficients is the refractive index, which is described by various expressions depending on the optical properties of the media filling the considered domain [1]. The same applies to the differential equations used to calculate the wave fronts and © Kuziv Ya. Yu., 2019 This work is licensed under a Creative Commons Attribution 4.0 International License http://creativecommons.org/licenses/by/4.0/ amplitudes in wave optics [2] or in the framework of the adiabatic waveguide method (AWM) proposed by A. A. Egorov and L. A. Sevastyanov [3-5]. The standard method for integrating nonlinear PDE is the characteristic method [6, 7], proposed in the 19th century by Cauchy. This method is implemented in many computer algebra systems (CAS), including Maple [8]. Although the formulation of the Cauchy method as such does not imply this, in all these implementations, the coefficients of the PDE are assumed to be given by unique analytical expressions throughout the entire considered domain. Therefore, to solve optical problems, it is necessary to develop new software capable of integrating first-order PDEs in the case when the coefficients of the equation are given by different analytical expressions in the subdomains G1, . . . , Gk , into which the considered domain is divided. This paper provides a mathematical description of this class of differential equations, offers the adaptation of the Cauchy integration method to this class of equations and its implementation as a Python program using the library SymPy [9]. Piecewise elementary partial differential equations We restrict our consideration to some domain G in the space Rn. We say that this domain is divided into subdomains G1, . . . , Gk if two conditions are met: the intersection of any pair of subdomains is empty, that is, Gi ∩ Gj = ∅ (i /= j), the closure of the union of these domains gives the closure of G, i.e., G = ∪Gi. We call such a partition elementary, if for any i = 1, . . . , k we can specify an elementary expression gi such that Gi = {x ∈ G : gi(x) < 0} . The function f with the domain of definition G will be called given piecewise elementary if an elementary partition of the domain G into several subdomains is specified, and in each of these domains an elementary expression is given for f , i.e., an elementary partition G1, . . . , Gk of the domain G and such elementary functions f1, . . . , fk exist such that f |Gi = f ∀i = 1, . . . , k. Now let x, y, z, p, q be five independent variables whose set of real values will be interpreted as a point in R5. Let f be a piecewise elementary function defined in some domain G of this five-dimensional space. Then the expression F (x, y, z, p, q) = 0, (1) where ∂z ∂z p = , q = , ∂x ∂y is a partial differential equation of the 1st order. We will refer to this kind of PDE as piecewise elementary. The algorithm for integrating a first-order PDE is not related in any way with the form for specifying the differential equation and consists in finding the characteristics [6, 7]. Recall that characteristics are curves in the five- dimensional space xyzpq, which are integral curves for the characteristic system of ordinary differential equations (ODE) for characteristics, which is composed based on the given function F using arithmetic operations and differentiation: where dx dy = Fp Fq dz = = p · Fp + q · Fq -dp = Fx + p · Fz - dq , (2) Fy + q · Fz Fx = ∂F , Fy = ∂x ∂F , Fz = ∂y ∂F , Fp = ∂z ∂F ∂F , Fq = . ∂p ∂q According to the well-known Cauchy theorem, a characteristic that inter- sects the plot of the solution of the equation (1) cannot leave this surface. The 19th century authors integrated the ODE system (2) and using well- known integrals reconstructed the solution of PDE (1). Of course, the class of ODEs that are integrable in a symbolic form is very small [10], therefore, in modern CAS this system is solved numerically, thus combining numerical and analytical methods for integrating the PDE (1). This technique can be described as follows. Let us know the curve C through which the desired solution of (1) should pass, and the value of p and q on the curve, or, equivalently, the curve in the space xyzpq, through which the manifold z = f (x, y), p = ∂f (x, y), q = ∂x ∂f (x, y) ∂y generated by the solution z = f (x, y) of the considered PDE should pass. Then the desired solution of PDE is woven from the characteristics released from this curve. Replace the curve C with a broken line and from each of its vertices let out an integral curve, solving the initial problem for the characteristic system of ODE numerically, for example, according to the Runge-Kutta method. In this case, we will use the coordinates x, y, z of the vertex and the corresponding values of p and q as the initial data. In the general case, as a result, we get a network in the five-dimensional space xyzpq, the projection of which into the space xyz gives the skeleton of the surface, which approximates the graph of the exact solution. Therefore, the construction of this projection can complete an approximate solution of the PDE with the data on the curve C. Degeneration happens only when the curve C itself is a characteristic. In the case when the PDE is piecewise elementary, it is necessary to make only a clarification in this scheme: when solving an ODE using the finite difference method, it is necessary at each step to find out to which domain the found point (x, y, z, p, q) of the integral curve belongs and to calculate the next point (xˆ, yˆ, zˆ, pˆ, qˆ) using the corresponding expression Fi instead of F . Unfortunately, the standard functions built in computer algebra systems, e.g., in CAS Maple, do not allow this refinement of their algorithm, and are essentially useless in solving optical problems. Therefore, we wrote our implementation of this algorithm in the Python programming language using the SymPy [9] library. The code is laid out github.com [11] and is available under a free license. We called this package MF_Solver_PDE. MF_Solver_PDE Our software is focused on finding a solution to a piecewise elementary PDE (1) passing through the curve C in the space xyzpq. The right side of the equation (1) is considered as a piecewise elementary function of five variables (x, y, z, p, q). This means that the space xyzpq is divided into several domains G1, . . . , Gk , which are specified by inequalities g1 < 0, . . . , gk < 0. For each of these domains, the symbolic expressions F1, . . . Fk are given such that F |Gi = Fi (i = 1, . . . , k). Thus, to set the left-hand side of the equation (1) means to set the inequal- ities describing the domains G1, . . . , Gk and the expressions F1, . . . , Fk . In the SymPy package, like in any computer algebra system, we can work with such data types. The input data for our algorithm for the solution of PDE are: the number k; the inequalities describing domains G1, . . . , Gk ; the expressions F1, . . . , Fk ; the parametric representation for the curve C in R5, through which the desired solution of PDE must pass. The auxiliary parameters are: N , the number of points that will be taken on the curve C with an equal parameter step; h, the ODE discretisation step. To find the PDE solution: we find the coordinates of N points on the curve C; for each of them we form the initial problem of finding the characteristic that goes out of this point, and solve it using the Runge-Kutta method. We used the explicit fourth-order method with 4 stages [12, 13]. As a result of performing the described steps, a set of points in R5 is obtained, their projection into xyz space giving a set of points lying on the plot of the desired solution z = z(x, y). By connecting with straight line segments the adjacent points lying on the same characteristic and the points found at the same step on two characteristics released from neighbouring points on the curve C, we get a two-dimensional skeleton of the surface. It can be turned into a surface plot using standard graphical tools available in SymPy. Examples and verification As the first example, consider the PDE az = pq. (3) We will search for its solution passing through a straight line given para- metrically in xyzpq space C : x = 1, y = τ, z = bτ, p = aτ, q = b. This solution can be written explicitly as z = (ax + b - a)y. The same solution obtained using our software is shown in Figure 1. The discrepancy between the numerical and analytical solutions is presented in Figure 2. It is clearly seen that the difference between numerical and analytical solutions is less than 10-7 and grows as the error is accumulated with each step. 1e-7 10 1.0 6 U 8 U 0.8 0.6 4 0.4 2 0.2 0 0.0 6 6 5 5 4 4 Y 3 2 1 X 0 1.0 1.1 1.2 1.3 1.4 1.5 1.6 Y 3 2 1 0 1.0 1.1 1.2 1.3 1.4 1.5 1.6 X Figure 1. The solution of Eq. (3) using our software Figure 2. Discrepancy between numerical and analytical solutions of Eq. (3) 0.4 0.2 y 0.0 -0.2 -0.4 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 x Figure 3. The surface describing the partial solution of the eikonal equation Figure 4. The surface describing the partial solution of the eikonal equation in the Maple program As the second example, let us consider a problem from geometrical optics. The eikonal equation has the form p2 + q2 = n2(x, y). Here n is the refractive index. Let for definiteness  2 2 n2 = 2 - x - y , x2 + y2 < 1, 1, x2 + y2 > 1. In the present example, it is assumed that we have two domains: a circular Lüneburg lens and the medium in which it is placed (e.g., air) [14, 15]. The solution obtained using our software is demonstrated in Figure 3. To check the correctness of the program work we compare the graphical solution, obtained using our program, with the Maple result. Certainly, in Maple the equations were solved only within the circle, and the boundary conditions have been chosen artificially. The Maple solution is shown in Figure 4. As seen from the figures, the solutions visually coincide. Conclusion The paper presents the original software for the numerical solution of first- order partial differential equations (PDE), the characteristic feature of which is that the coefficients of equations in different domains G1, . . . , Gk are given by different analytical expressions F1, . . . Fk such that F |Gi = Fi ∀i = 1, . . . , k. The software was tested on benchmark problems taken from geometrical optics. A comparison with the results obtained by Maple was carried out. The considered examples are illustrative. The software is potentially suitable for solving any first-order PDEs, the left-hand parts of which are piecewise described by very long symbolic expressions. Very complex equations of this type arise, for example, in the framework of the adiabatic waveguide method (AWM) proposed by A. A. Egorov and L. A. Sevastianov [3-5]. In this theory, the equation that plays the same role as the eikonal equation in geometrical optics is obtained by equating the determinant of a 8×8 matrix to zero. Even composing the symbolic expression for the right-hand side turns out to be a challenge for computer algebra systems. However, as soon as this expression is found, the proposed algorithm allows constructing the propagation of electromagnetic waves of the AWM in media with piecewisespecified refractive index. It is this application that we see as the natural scope of the software presented in this article.

About the authors

Yaroslav Yu Kuziv

Peoples’ Friendship University of Russia (RUDN University)

Author for correspondence.
Email: yaroslav.kuziw@yandex.ru

postgraduate student of Department of Applied Probability and Informatics

6, Miklukho-Maklaya str., Moscow, 117198, Russian Federation

References