CALCULATION OF SPECIAL FUNCTIONS ARISING IN THE PROBLEM OF DIFFRACTION BY A DIELECTRIC BALL Текст научной статьи по специальности «Математика»

Research article

UDC 519.872:519.217

PACS 07.05.Tp, 02.60.Pn, 02.70.Bf

DOI: 10.22363/2658-4670-2021-29-2-146-157

Calculation of special functions arising in the problem of diffraction by a dielectric ball

Ksaverii Yu. Malyshev

Skobeltsyn Institute of Nuclear Physics Lomonosov Moscow State University GSP-1, Leninskie Gory, Moscow, 119991, Russian Federation

(received: April 29, 2021; accepted: May 25, 2021)

To apply the incomplete Galerkin method to the problem of the scattering of electromagnetic waves by lenses, it is necessary to study the differential equations for the field amplitudes. These equations belong to the class of linear ordinary differential equations with Fuchsian singularities and, in the case of the Lüneburg lens, are integrated in special functions of mathematical physics, namely, the Whittaker and Heun functions.

The Maple computer algebra system has tools for working with Whittaker and Heun functions, but in some cases this system gives very large values for these functions, and their plots contain various kinds of artifacts. Therefore, the results of calculations in the Maple'11 and Maple'2019 systems of special functions related to the problem of scattering by a Lüneburg lens need additional verification. For this purpose, an algorithm for finding solutions to linear ordinary differential equations with Fuchsian singular points by the method of Frobenius series was implemented, designed as a software package Fucsh for Sage. The problem of scattering by a Lüneburg lens is used as a test case. The calculation results are compared with similar results obtained in different versions of CAS Maple.

Fuchs for Sage allows computing solutions to other linear differential equations that cannot be expressed in terms of known special functions.

Key words and phrases: linear differential equations, Whittaker functions, Heun functions

1. Introduction

The problem of diffraction of a plane electromagnetic wave by a ball with an arbitrary radially symmetric filling allows the construction of an analytical solution by the incomplete Galerkin method [1]. Let us make use of a spherical coordinate system and assume that the dielectric constant of the ball e depends only on r, and the permeability ^ is constant. Let us expand the Borgnis

© Malyshev K. Yu., 2021

This work is licensed under a Creative Commons Attribution 4.0 International License http://creativecommons.org/licenses/by/4-0/

potentials u, v for electric- and magnetic-type fields in terms of the system of spherical functions:

oo

(r)P^ (0) sin 4>, v = ^vn (r)P^ (0) sin 0.

n=1 n=l

Functions un(r), vn(r) satisfy linear differential equations of the second order: for TE-polarized wave

d 1 du,r

+

dr [i dr and for TM-polarized wave

d 1 dvn ---- +

dr e dr

^ n(n+1)

k2 e —

¡ir2

un 0,

„2.. n(n+1)

k2 ^ —

erz

vn = 0.

(1)

(2)

Not only in the classical case of a ball with constant filling, but also in the case of a Lüneburg lens, when

\ 2-r2, r <1, {1, r>1,

these equations are integrated in special functions of mathematical physics. In the classical case, the cylindrical functions are obtained [2]-[4]. For a Lüneburg lens, (1) was integrated in the paper by Lock [5] in terms of Whittaker functions. In the Maple system [6] both equations (1) and (2) for a Lüneburg lens are integrable: the solution to Eq. (1) is expressed in terms of Whittaker functions and the solution to Eq. (2) in terms of Heun functions.

When working with these functions in the Maple system, we faced the fact that they take on huge values that grow indefinitely with an increase in n, and all sorts of artifacts (gaps and ripples) appear in the plots. This observation raised doubts about the adequacy of the algorithms used in Maple.

However, it is easy to see that both equations (1) and (2) are linear differential equations of the 2nd order that have a Fuchsian singular point [7] at r _ 0. Therefore, their solutions upon the analytic dependence of e on r can be expanded in the Frobenius series [7]. The coefficients of the equations (1) and (2) have finite nonzero singular points only where the functions e and 1/t have them, respectively. In particular, if e is a polynomial, then the Frobenius series for the solution to Eq. (1) converges for all r, and the Frobenius series for the solution of Eq. (2) converges in a circle, on the boundary of which the complex zero with the least modulus lies [7]. For example, in the case of the Lüneburg lens, the Frobenius series for the solution to Eq. (1) converges for all r, and for the solution of Eq. (2), the radius of convergence is \[2. Therefore, to find the field in the lens, it is quite sufficient to calculate the coefficients of the Frobenius series for the solutions of Eqs. (1) and (2).

2. Fuchs for Sage package

In general, the problem of finding Frobenius series can be formulated as follows: for a linear differential equation

x2y" + xpy' + qy = 0, (3)

with meromorphic functions p, q, construct a solution in the form of a series

y = xr(1 + u1 X + u2x2 + ... ), where r is a root of characteristic equation

r(r-1)+p(0)r + q(0) = 0. (4)

Usually, both roots of this equation correspond to solutions, but in some cases, when choosing a root with a smaller real part, a solution in the form of a series does not exist. In this case, there is one solution in the form of a Frobenius series, and the second solution is obtained from the Liouville formulas connecting two solutions of a second-order linear differential equation [7]. The solution to this problem was implemented in CAS Sage [8], as the Fuchs for Sage software package.

Given symbolic expressions p and q, the function fuchs_order(p,q) returns the roots of Eq. (4) as elements of the algebraic number field.

sage: load('fuchs.sage') None

sage: fuchs_order(1,x~2-3~2) [-3, 3]

sage: fuchs_order(1,x~2-8) [-2.828427124746190?, 2.828427124746190?]

Function fuchs_order(p,q,r,N) returns symbolic expression

y = xr (1 + u1 x + u2 x2 + —+ unxn ),

given symbolic expressions p and q, order r, and number N of the sought series terms.

This function supports operation with non-rational roots:

sage: R=fuchs_order(1,x~2-8) sage: fuchs_series(1,x~2-8,R[1],10)

(-(1.41301303805500?e-9)*x~10 + (2.212333918945947?e-7)*x'8

- 0.00002417081750580247?*x~6 + 0.001690534180537310?*x'4

- 0.06530096874093536?*x~2 + 1)*x~2.828427124746190?

Let us give an example of calculating the solution of Eq. (1) for the case of a Lüneburg lens with numerical values of the parameters k, n:

sage: k=4

sage: n=3

sage: p=0

sage: q=k~2*x~2*(2-x~2) -n*(n+1)

sage: fuchs_order(p,q) [-3, 4]

sage: expand(fuchs_series(p,q,4,10))

-3712/19305*x~14 + 1928/3861*x"12 - 448/429*x~10 + 164/99*x"8 - 16/9*x~6 + x"4

The terms of this series do not form a monotonie sequence (see figure 1), therefore, when summing the series, it is important to take a sufficient number of terms.

UnXn + r

6el4

4el4

2el4

-2el4

-4el4

-6el4

■50

t-I«

100

n

150

200

250

300

Figure 1. Solution for the Lüneburg lens: distribution of the values of the terms of the series over n — the summand number of the Erobenius series. Calculations are performed

for the point x = 4

3. Calculation of special functions in the Lüneburg case

Fuchs for Sage allows describing the partial solutions of Eqs. (1) and (2), bounded at zero, in the form of the first N terms of the Frobenius series. In our experiments, a moment was revealed after which a further increase in the number of terms did not lead to a noticeable change in the sum at the reference points. Of course, it would be very convenient to get not greatly overstated theoretical estimates for the number of terms that must be kept in the series in order to preserve the chosen accuracy.

In order to exclude the influence of round-off errors, calculations are carried out in the field of rational numbers, the summation of the series was carried out at rational points of the real axis. For visibility, the rational values obtained without round-off errors are presented in the plots below in logarithmic scale.

Let us compare the results of computer experiments in Sage and Maple for the case of a Lüneburg lens. As noted above, in this case, solutions in the form of a series can be compared with solutions in special functions of mathematical physics that are used in the Maple system.

3.1. TE field

In the case of the Lüneburg lens, the coefficients of Eq. (1) have only two singular points (r = 0 and r = œ) and therefore this series can be expressed in terms of the degenerate hypergeometric function or the Whittaker function [9, §13.14] [10]. The Maple system uses the expression

1 „T1 . , /fc 2rc + 1 , 2\ Urn = — WhittakerM ( -,-, kr2 ) .

Vr \2 4 ' )

This solution to Eq. (1) differs from those obtained using the Frobenius method by the multiplicative constant, which can be found in symbolic form using the Taylor series expansion.

0.4 0.6 0.8 1.0 1.2 1.4

Figure 2. Function ln \un(x)\ in case 3, k = 40, n = 25, plotted in Sage (left) and Maple'2019 (right). In Fuchs for Sage we took 400 terms and the calculations were

performed at 221 points

We compared the results of calculating the solution by the Frobenius method and by means of Maple'2019 for 3 representative cases: 1) k = 4, n = 3,2) k = 10, n = 3, 3) k = 40, n = 25. Despite some difficulties associated with the difference in the automatic selection of proportions of the plots, in all cases the calculation using our package leads to the same results as using Maple, at least with graphic accuracy, see, e.g., figure 2. For a more accurate analysis, we calculated the values at the reference points, see Table 1. It is clearly seen that the values found in Sage and Maple match up to a round-off error. Thus, the numbers to which the Frobenius series converge coincide with the results of the built-in Maple algorithms with high accuracy. Thus, both the algorithm for calculating the Frobenius series and the high orders of magnitude of the sought functions observed in numerical experiments are verified.

Values of function In \un(x)\

Case 1: k = 4, n = 3.

x 1/2 1 3/2 3

N 50 100 100 200

Sage -3.21179025713161 -1.66169222058965 -1.32135176796381 8.42740419841879

Maple -3.211790257 -1.661692220 -1.321351768 8.427404199

Case 2: k = 10, n = 3.

x 1/2 1 3/2 3

N 50 100 100 500

Sage -7.92874739435119 -5.90970510172370 -5.14868973498956 51.3761521918252

Maple -7.928747395 -5.909705102 -5.148689031 51.37615219

Case 3: k = 40, n = 25.

x 1/2 1 3/2 3

N 100 300 400 1000

Sage -26.5511644045404 -27.0297223651032 -23.8706312720468 79.7748111646792

Maple -26.55116440 -27.02972236 -23.87063127 79.77481117

Table 2

Values of function vn(x). Case 1: k = 10, n = 3

x 1/2 3/4 1 7/5

N 100 100 300 700

Sage 0.000368374895008100 -0.00110378333099445 0.00201228022550062 -0.00181338274195031

Maple 0.0003683748946 -0.001103783338 0.002012279981 -0.001813534799

Table 3

Partial sums of Erobenius series for vn(x) at the singular point x = a/2. Case 1: k = 10, n = 3

N 100 300 700 1000

Sage -0.00181372503052628 -0.00181282882557290 -0.00181273282789891 -0.00181272165872556

3.2. TM field

Equation (2) has 4 singular points in the complex plane: r = 0, r = and r = <x>. Therefore. its solutions are expressed in terms of the Heun confluent function [9. §31.12] [11]. In Maple the expression

vn = ekr2/2HeunC(2A;, n + 1/2, -2, k2, -k2 + 3/4, r2/2)rn+1 (5)

is used. To find the multiplicative constant. by analogy with above. we tried to find the first non-vanishing term of the Taylor series of this function in the vicinity of the point r = 0. Unfortunately. the standard function taylor cannot be applied to the vn function for expansion in Taylor series.

Remark. An attempt of using taylor to find the first three partial sums of the Taylor series at the point r = 0 leads to the result 0(r4). An attempt to find the fourth partial sum at the point r = 0 leads to an error Error, (in series/function) unable to compute series. We observed this kind of failure in all test examples without exception. When choosing a point other than zero as the center of the expansion. the function taylor works without errors.

Nevertheless. computer experiments show that the function (5). calculated in Maple. coincides with a high accuracy with the sum of the Frobenius series. see Table 2. This makes us to assume that the desired constant is 1.

X

Figure 3. Plot of function ln lvn(x)l at k = 18. n = 1; on the right - Frobenius series with N = 200 terms. calculated at 41 points. on the left - plot from Maple'2019

Unlike Eq. (1). Eq. (2) has a finite singular point. namely x = \[2. We will

now show that for x > \[2 the Maple system produces incorrect results for representing the function vn (x).

The Frobenius exponents in the vicinity of the Fuchsian singular point x = \f2 of Eq. (2) are equal to 0 and 2. and do not depend on the parameters n and k. In the first case. the application of our function leads to division by zero. which indicates the presence of a logarithmic singularity in the expression of the desired function [7]. [11]. In this case. the solution corresponding to the Frobenius exponent equal to 0 is bounded at the singular point x = \[2. In

the second case, according to Fuchs theorem [7] we obtain a Frobenius series with a nonzero radius of convergence. Therefore, the function vn (x) continues analytically beyond the point x = V2, having a removable singularity at this point for real values of x. The plot of the Frobenius series sum (figure 3) fully confirms that vn has a finite limit at the point x = \[2\ it cannot give more, since it diverges at x > \[2.

A numerical experiment in Fuchs for Sage indicates the convergence of the Frobenius series for the function vn(x) at the point x = \[2 in the field of real numbers, see Table 3. However, this convergence is noticeably slower than at rational points x <\[2.

The situation in Maple looks much less clear. Substitution of the value x = \[2 in the Maple system yields the result Float(infinity)+Float(infinity)*I for all the values of k,n we checked. In this case, the plot of vn(x), obtained by means of Maple'2019, is cut off at

this point. Substitution of values x > \[2 results in complex numbers of astronomical orders of magnitude, for example evalf(eval(V, r = 3/2)) yields -5.630587096*10~12 + (1.389900117*10~15)*I, for the parameter values k=10,n = 3.

Moreover, when calculating the Heun function in Maple, problems appear that are characteristic of the accumulation of round-off errors when summing a power series. In older versions of Maple, these problems were encountered for all considered parameters, and their first appearance took place already at k = 10, see figure 4. In new versions, the same shortcomings in the display of plots appear in a different range of parameters: the highest harmonic n = 1 is displayed incorrectly at k > 18 (see figure 3, left, and 5). At the same time, even a relatively small number of terms of the Frobenius series makes it possible to obtain a smooth curve without artifacts, see figure 3, right.

Figure 4. Heun functions vn for values Figure 5. Calculation of functions

k = 10, n= 1,2,3 in Maple'11 ln K(x)|. Case 3, k = 40, n = 25

0.8

0.85

0.9

0.95

1.05

1.1

1.15

1.2

-1C

3C

2C

1C

c

Figure 6. Plot of function ln \vn(x)\ at k = 40, n = 25, left — plot from Maple'2019, right — Frobenius series with 250 terms and 82 values of x used

In this paper, we present the Fuchs for Sage package , which allows calculating the solutions of second-order equations having the Fuchsian singularity with high accuracy and, in particular, looking for the potentials of the fields scattered by a ball with a wide class of dielectric constant radial dependences. This is confirmed by a comparison with the case of a Lüneburg lens, when the solution is expressed in terms of the known functions of mathematical physics.

At the same time, we made sure in an independent way that the solutions of the equations (1) and (2) in the case of a Lüneburg lens can change by many orders, as it happened in Maple.

The results of our research can be useful for analyzing the electromagnetic field in the vicinity of the focus of the Lüneburg lens [5].

It should be especially noted that the symbolic expression for TM fields in terms of Heun confluent functions found in Maple is dangerous for use, since it gives results that are not completely correct from a theoretical point of view. At the same time, even a small number of terms of the Frobenius series allows good approximations to the values of the sought functions, see figure 6.

The author thanks M. D. Malykh (RUDN University) for setting an interesting problem and constant attention to this work. The calculations presented in the article were performed in systems Maple [6] and Sage [8].

[1] A. G. Sveshnikov and I. E. Mogilevsky, Matematicheskiye zadachi teorii difraktsii [Mathematical problems of the theory of diffraction]. Moscow: MSU, 2010, In Russian.

4. Conclusion

Acknowledgments

References

[2] G. Mie, "Beiträge zur Optik trüber Medien, speziell kolloidaler Metal-lösungen," Annalen der Physik, vol. 25, no. 3, pp. 377-445, 1908. DOI: 10.1002/andp.19083300302.

[3] C. F. Bohren and D. R. Huffman, Absorption and scattering of light by small particles. New York: John Wiley & Sons, 1983.

[4] C. Mätzler, MATLAB Functions for Mie Scattering and Absorption, Institute of Applied Physics, University of Bern, June 2002. Research Report No. 2002-08 http://arrc.ou.edu/~rockee/NRA_2007_ website/Mie-scattering-Matlab.pdf, 2002.

[5] 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, vol. 25, pp. 2980-2990, 2008. DOI: 10.1364/JOSAA.25.002980.

[6] Waterloo Maple (Maplesoft), Symbolic and numeric computing environment Maple, https://www.maplesoft.com/, 2019.

[7] F. G. Tricomi, Differential equations. London: Blackie & Sons ltd., 1961.

[8] The Sage Developers, SageMath, the Sage Mathematics Software System (Version 7.4), https://www.sagemath.org, 2016.

[9] National Institute of Standards and Technology (NIST), United States, Digital Library of Mathematical Functions. Version 1.1.1, https://dlmf. nist.gov, 2021.

[10] A. F. Nikiforov and V. B. Uvarov, Special Functions of Mathematical Physics. A Unified Introduction with Applications. Springer Basel AG, 1988.

[11] S. Y. Slavyanov and W. Lay, Special functions: unified theory based on singularities. Oxford: OUP, 2000.

For citation:

K. Yu. Malyshev, Calculation of special functions arising in the problem of diffraction by a dielectric ball, Discrete and Continuous Models and Applied Computational Science 29 (2) (2021) 146-157. DOI: 10.22363/2658-46702021-29-2-146-157.

Information about the authors:

Malyshev, Ksaverii — engineer, Skobeltsyn Institute of Nuclear Physics, Lomonosov Moscow State University (e-mail: kmalyshev08102@mail.ru, ORCID: https://orcid.org/0000-0001-8823-9136)

УДК 519.872:519.217

PACS 07.05.Tp, 02.60.Pn, 02.70.Bf

DOI: 10.22363/2658-4670-2021-29-2-146-157

О вычислении специальных функций, возникающих при исследовании задачи дифракции на диэлектрическом шаре

К. Ю. Малышев

Научно-исследовательский институт ядерной физики имени Д. В. Скобельцына Московский государственный университет имени М.В. Ломоносова Ленинские горы, д. 1, стр. 2, Москва, ГСП-1, 119991, Россия

Для применения неполного метода Галёркина к задаче о рассеянии электромагнитных волн на линзах необходимо исследовать дифференциальные уравнения для амплитуд полей. Эти уравнения принадлежат к классу линейных обыкновенных дифференциальных уравнений с фуксовыми особенностями и, в случае линзы Люнеберга, интегрируются в специальных функциях математической физики — функциях Уиттекера и Гойна.

В системе компьютерной алгебры Maple имеются инструменты для работы с функциями Уиттекера и Гойна, однако в ряде случаев эта система выдаёт очень большие значения для этих функций, а их графики содержат разного рода артефакты. Поэтому результаты вычислений в системах Maple'11 и Maple'2019 специальных функций, связанных с задачей рассеяния на линзе Люнеберга, нуждаются в дополнительной проверке. С этой целью был реализован алгоритм нахождения решений линейных обыкновенных дифференциальных уравнений с фуксовыми особыми точками методом рядов Фробениуса, оформленный в виде пакета программ Fuchs for Sage. Задача рассеяния на линзе Люнеберга используется в качестве тестового примера. Результаты расчётов сопоставляются с аналогичными результатами работы в CAS Maple разных версий.

Пакет Fucsh for Sage позволяет вычислять решения и других линейных дифференциальных уравнений, решения которых не выражаются через известные специальные функции.

Ключевые слова: линейные дифференциальные уравнения, функции Уиттекера, функции Гойна