CC BY
28
Discrete & Continuous Models 2019 27 (1) 42—48
. „ . http://journals.rudn.ru/miph
& Applied Computational Science
Research article
UDC 512, 519.6
DOI: 10.22363/2658-4670-2019-27-1-42-48
Software for the numerical solution of first-order partial
differential equations
Yaroslav Yu. Kuziv
Department of Applied Probability and Informatics Peoples' Friendship University of Russia (RUDN University) 6, Miklukho-Maklaya str., Moscow, 117198, Russian Federation
(received: June 29, 2019; accepted: October 21, 2019)
Partial differential equations of the first order, arising in applied problems of optics and optoelectronics, often contain coefficients that are not defined by a single analytical expression in the entire considered domain. For example, the eikonal equation contains the refractive index, which is described by various expressions depending on the optical properties of the media that fill the domain under consideration. This type of equations cannot be analysed by standard tools built into modern computer algebra systems, including Maple.
The paper deals with the adaptation of the classical Cauchy method of integrating partial differential equations of the first order to the case when the coefficients of the equation are given by various analytical expressions in the subdomains G1,..., Gk, into which the considered domain is divided. In this case, it is assumed that these subdomains are specified by inequalities. This integration method is implemented as a Python program using the SymPy library. The characteristics are calculated numerically using the Runge-Kutta method, but taking into account the change in the expressions for the coefficients of the equation when passing from one subdomain to another. The main functions of the program are described, including those that can be used to illustrate the Cauchy method. The verification was carried out by comparison with the results obtained in the Maple computer algebra system.
Key words and phrases: eikonal, partial differential equation, SymPy
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 G\,... ,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].
1. 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:
1) the intersection of any pair of subdomains is empty, that is, Gi n Gj = 0
(i = j),
2) the closure of the union of these domains gives the closure of G, i.e.,
G = uGi.
We call such a partition elementary, if for any i = 1,... ,k we can specify an elementary expression gi such that
Gi = {x e 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 f 1,..., fk exist such that
fk = f Vi =l,...,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
dz dz
p = q =
dx dy
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:
dz
—dp
—dq
dx dy
Fp = Fq = p ■ Fp + q ■ Fq = Fx + p ■ Fz = Fy + q ■ Fz
(2)
where
F = dF, F = dF, F = dF, F = dF, F = dF.
x dx' y dy ' z dz ' p dp' q Sq '
According to the well-known Cauchy theorem, a characteristic that intersects 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 = /(x,y^ p = dx(x,y)'
df ( ) q = dy (x,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.
2. 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 = l,...,k).
Thus, to set the left-hand side of the equation (1) means to set the inequalities 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:
1) we find the coordinates of N points on the curve C;
2) 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.
3. 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 = r, z = br, p = ar, 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.
Figure 1. The solution of Eq. (3) using Figure 2. Discrepancy between numerical our software and analytical solutions of Eq. (3)
0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
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
0.4
0.2
0.0
0.2
0.4
x
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
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.
The paper presents the original software for the numerical solution of firstorder 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 Vi = l,..., 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 8x8 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 piecewise-specified refractive index. It is this application that we see as the natural scope of the software presented in this article.
The publication has been prepared with the support of the "RUDN University Program 5-100".
1. M. Born, E. Wolf, Principles of optics. Electromagnetic theory of propagation, interference and diffraction of light, 6th Edition, Elsevier, 1980.
2. B. Enquist, O. Runborg, Computational high frequency wave propagation, 6th Edition, Cambridge University Press, 2003.
x2 + y2 < 1, x2 + y2 > 1.
Conclusion
Acknowledgments
References
3. A. A. Egorov, L. A. Sevastianov, Structure of modes of a smoothly irregular integrated-optical four-layer three-dimensional waveguide, Quantum Electronics 39 (6) (2009) 566. doi:10.1070/QE2009v039n06ABEH013966.
4. A. A. Egorov, K. P. Lovetsky, L. A. Sevastianov, A. L. Sevastianov, Simulation of guided modes (eigenmodes) and synthesis of a thin-film generalised waveguide Luneburg lens in the zero-order vector approximation, Quantum Electronics 40 (9) (2010) 830-836. doi:10.1070/QE2010v040n09ABEH014332.
5. A. A. Egorov, L. A. Sevastianov, A. L. Sevastianov, Method of adiabatic modes in research of smoothly irregular integrated optical waveguides: zero approximation, Quantum Electronics 44 (2) (2014) 167-173. doi:10.1070/QE2014v044n02ABEH015303.
6. A. D. Polyanin, V. E. Nazaikinskii, Handbook of linear partial differential equations for engineers and scientists, 2nd Edition, CRC Press, Boca Raton, London, 2016.
7. E. Goursat, Cours d'analyse mathématique, 3rd Edition, Vol. 2, Gauthier-Villars, Paris, 1918.
8. PDEplot for Maple (2019).
URL http://www.maplesoft.com
9. Python library for symbolic mathematics SymPy (2019). URL http://www.sympy.org
10. M. D. Malykh, On integration of the first order differential equations in finite terms, IOP Conf. Series: Journal of Physics: Conf. Series 788 (2017) 012026. doi:10.1088/1742-6596/788/1/012026.
11. IarKuz at Github (2019).
URL https://github.com/IarKuz/PostGradeCode/blob/MF_Solver_ PDE/MF_Solver_PDE.ipynb
12. E. Hairer, G. Wanner, S. P. N0rsett, Solving ordinary differential equations, 3rd Edition, Vol. 1, Springer, New York, 2008.
13. W. H. Press, S. A. Teukolsky, W. T. Vetterling, B. P. Flannery, Numerical recipes in C: the art of scientific computing, 2nd Edition, Cambridge University Press, 1992.
14. J. Lock, Scattering of an electromagnetic plane wave by a Luneburg lens. II. Wave theory, Journal of the Optical Society of America A: Optics Image Science and Vision 25 (2008) 2980-2990. doi:10.1364/J0SAA.25.002980.
15. S. Cornbleet, Geometrical optics reviewed: A new light on an old subject, Proceeding of the IEEE 71 (1983). doi:10.1109/PR0C.1983.12620.
For citation:
Ya. Yu. Kuziv, Software for the numerical solution of first-order partial differential equations, Discrete and Continuous Models and Applied Computational
Science 27 (1) (2019) 42-48. DOI: 10.22363/2658-4670-2019-27-1-42-48.
Information about the authors:
Yaroslav Yu. Kuziv (Russian Federation) — postgraduate student of
Department of Applied Probability and Informatics of Peoples' Friendship
University of Russia (RUDN University) (e-mail: yaroslav.kuziw@yandex.
ru, phone: +7(495)9550927)
