The mathematical modeling of epidemics has a long history; it began with the Kermack-McKendrick model [8] , introduced in 1927. In this seminal paper, the whole population is divided into Susceptible, Infectious, and Recovered sub-populations. Then, some ordinary differential equations are formulated specifying the time evolution of the functions representing these sub-populations. Wavelet analysis is now frequently used to extract information from epidemiological and other time series. Grenfell et al. [7] introduced wavelet analysis for characterizing non-stationary epidemiological time series.
Cazelles et al. [1] use the Morlet wavelets for applications in epidemiology. Lavrova et al. [10] modeled the disease dynamics caused by Mycobacterium tuberculosis in Russia using a sum of two logistic functions (3) (bi-logistic model).
SARS-CoV-2 initially emerged in China, at the end of 2019; after Chinese scientists identified the sequence of the new virus [17] , this information was shared with the international community. Since then, a lot of articles were written and published, describing from different points of view, the new SARS-CoV-2 coronavirus and the COVID-19 disease, caused by the virus. We will point out only some of them. Fokas et al. FDK used a generalization of the logistic function for forecasting the number of individuals reported to be infected with SARS-CoV-2 in different countries. Krantz et al. [9] proposed a two-phase procedure (combining discrete graphs and Meyer wavelets) for constructing true epidemic growth. A method similar to that one from Lavrova et al. [10] was used by E. Vanucci and L. Vanucci [16] for predicting the end date of Covid-19 disease in Italy.
The outline of the present paper is as follows. In Sec. 2, we discuss the basic properties of Riccati's equation, logistic equation, and logistic curve. For this purpose, we use Eulerian numbers. Sec. 3 and Sec. 4 are devoted to logistic wavelets. In Sec. 5 we model the cumulative number of persons reported to be infected by SARS-CoV-2 in the United States as a sum of several logistic functions.
We use the following convention for the Fourier transform:
where f ∈ L 1 (R) ∩ L 2 (R). 
The logistic equation is defined as
where t is time, u = u(t) is the unknown function, s, u max are constants. The constant u max is called the saturation level. The integral curve u(t) fulfilling condition 0 < u(t) < u max is known as the logistic function. After solving (2) we get the logistic function in the following form
where t 0 is its inflection point, which is related to the initial condition u(0) = u 0 = u max 1 + e st0 , therefore
Equation (2) is a particular case of Riccati's equation with constant coefficients
The constants r = 0, u 1 , u 2 can be generally real or complex numbers. If u(t) is a solution of (4) then it is known a formula for the nth derivative u (n) (t) (n = 2, 3, 4, . . .) of u(t) expressing it as a polynomial of the function u(t) itself:
where n = 2, 3, . . . and n k denotes the Eulerian number (number of permutations of the set {1, 2, . . . , n} having k, (k = 0, 1, 2, . . . , n − 1) permutation ascents, see Graham et al [6] ). The first few Eulerian numbers are given in the Table 1 . Formula (5) was discussed during the Conference ICNAAM 2006 (September 2006) held in Greece and it appeared, with an inductive proof, in paper [11] (see also [12] ). Independently the formula has been considered and proved, with the proof based on generating functions, by Franssens [5] . The polynomial of u, of order (n + 1), appearing on the right-hand side of (2) is known in the literature as a kind of the so-called derivative polynomials. It is easy to see that all (n + 1) roots of the polynomial are simple and lie in the interval [u 1 , u 2 ]. The derivative polynomials have been recently intensively studied.
Formula (5) applied to the particular case of the logistic equation (2) is as follows:
The polynomial of the variable u and of order (n + 1) on the right hand side of (6) is uniform in the sense of the following. Remark 1. If u 0 is a root of the polynomial on the right hand side of (6), i.e.,
then dividing both sides of (7) by u n+1 max we get
Thus u 0 is a root of the derivative polynomial on the right hand side of (6) if u 0 /u max is the root of the polynomial
Let us write down, using formula (6) and the notation of (8), the first few derivatives of the logistic function, which fulfills equation (2) . By Remark 1 we can assume, without loss of the generality, that u max = 1 and s = 1.
We obtain successively:
All roots of the polynomials P k (u) for k = 3, 4, 5, 6 can be calculated explicitly, so the polynomials can be factored and we get
Therefore the minimal positive root of the polynomial
Thus by using Remark 1 we see for example that if at a minimal time t 1 , u ′′′ (t 1 ) = 0 (t 1 is simultanously a maximum of u ′′ (t)) then the value of the logistic function at this point is
Similar conclusions can be drawn for the smallest zero of the u (4) (t) (polynomial P 5 (u)) or u (5) (t) (polynomial P 6 (u)) using constants (9).
Let a wavelet ψ 2 (x) (see Figure 1 ) be the second derivative of the logistic function u(x) = 1 1+e −x . Since u ′ (x) = −u(u − 1), then by (5) or directly we get
and by (11) it follows that the wavelet has the following exact form
Changing the variable u = 1 1+e −x , u ′ (x) = u(1 − u) in the following three integrals we calculate
which proves that ψ 2 (x) ∈ L 1 (R) ∩ L 2 (R). In fact ψ 2 (x) ∈ S(R) (the space of rapidly decreasing functions on R). We will discuss this in the next section. By ψ + 2 (x) we denote the positive part of ψ 2 (x), i.e.,
The Fourier transform of ψ 2 (x) is as follows:
It is well known (see [2] ) that a wavelet ψ(x) ∈ L 1 (R) ∩ L 2 (R) should satisfy the following admissibility condition
We will show that for ψ 2 (x) the condition (15) is satisfied and even the integral can be expressed in a closed form in terms of the Riemann zeta function. Namely, using (14) and the following formula from Dwight's Tables [3] (item no 860.519):
we have
We generate a doubly-indexed family of wavelets from ψ 2 by dilating and translating,
where a, b ∈ R, a > 0 and denote by ψ + a,b 2 (x) the positive part of ψ a,b 2 (x).
Similarly as in the previous section we define a wavelet ψ n (x) to be the nth (n = 3, 4, ...) derivative of the logistic function u(x) = 1 1+e −x . Figure 2 shows graph of the wavelet ψ 3 (x). Thus (5) gives 
By definition, the function ψ n (x) is an even function for odd n and an odd function when n is even. The numerator of the expression (18) is a polynomial of degree n of the variable e −x , while the denominator of degree n + 1. Therefore for any polynomial p(x) we have lim x→−∞ p(x)ψ n (x) = 0. Since ψ n (x) has the symmetry property then also lim x→∞ p(x)ψ n (x) = 0. The last conclusion can also be drawn from multiplying the numerator and the denominator of (18) by e (n+1)x . From this and from the fact that ψ k+1 (x) = ψ ′ k (x) for any integer k ≥ 2 it follows that ψ n (x) ∈ S(R), (n = 2, 3, . . .).
By (14) we haveψ
Now using (19) and once again formula (16) we can calculate the integral of the admissibility condition (15) as follows:
As usually we generate a doubly-indexed family of wavelets from ψ n by dilating and translating,
where a, b ∈ R, a > 0, n = 2, 3, . . ..
Denote by y * n total cumulative number of individuals reported to be infected up to nth day in a country or a region and by y n the 7-day central moving arithmetic average for the sequence y * n , i.e.,
We will look, in the sequence (y n ), for points corresponding to the zeros of the second or the third derivative of the logistic function. This is equivalent to detect the points, where the sequence of second differences, ∆ 2 y n = y n+1 − 2y n + y n−1 , takes a value close to zero or a maximum respectively. We will find these points either directly by observing the sequence of second differences (∆ 2 y n ) or detect them by using the wavelet ψ 2 (x) and its positive part ψ + 2 (x). From the considerations in Sec. 2 and from (10) it follows that parameter b should be determined as that point where the sequence (∆ 2 y n ) changes sign. Parameter a should be chosen in such a way that the distance between the zero and the maximum of (∆ 2 y n ) was approximately 1.319a. Thus, we obtain two parameters defining the first logistic function (first wave) approximating the time series (y n ). It remains to determine the third parameter of the first wave, i.e., its saturation level y max . Assuming that (y n ) initially follows a logistic function y n ≈ y(n) = y max
and since by definition it holds
then by (13) we get successively n ∆ 2 y n ψ + a,b
Using (21) we can estimate y max as follows
Parameters a and b can also be estimated by maximizing locally the integral on the left-hand side of (21). Thus we find in the sequence ∆ 2 y n the best pattern corresponding to the positive part of the wavelet ψ a,b 2 . To avoid the situation that the next wave, immediately following the previous one, could distort our findings we use here the positive part ψ + a,b 2 , not the whole wavelet ψ a,b 2 . The saturation level of the first wave can also be estimated as twice the value of the sequence (y n ) at the point where (∆ 2 y n ) changes signs (inflection point) or its maximal value multiplied by 1/0.211 (zero of the third derivative). Having found the values of the parameters a, b, and y max for the first wave, we create a new time series by subtracting the first wave from y n , i.e., z n = y n − y max 1 + exp(− n−b a )
, and with the sequence z n we proceed in the same way as with y n , calculating successive logistic waves. After this we use the nonlinear Generalized Reduced Gradient method to optimize the values of saturation levels (but not a's and b's). All data were collected from the https://www.worldometers.info/coronavirus/ platform. 
Let us use the theory to build a model for the total cumulative number of individuals reported to be infected by SARS-CoV-2 in the USA. We assumed the observation period of 189 days, from March 13 (n = 1, the first day when the number of cases exceeded 2, 000) to September 17, 2020 (n = 189). All calculations were performed in Excel. Using the above-described procedure we have got the approximating function as a sum of the following waves f (x) = .
(23) Figure 3 shows the positive part of the wavelet (scaled) ψ + a,b 2 , a = 25, b = 6.6 fitted to the second differences (∆ 2 y n ). Figure 4 shows the total cumulative number of individuals reported to be infected by SARS-CoV-2 in the USA in the period: March 13 (n = 1) -September 17 (n = 189), (blue points) and the approximating function f (x) (23), (red points).
The |y n − f (n)| y n + f (n) = 0.01195.
The same procedure applied to the total cumulative number of individuals reported to be infected by SARS-CoV-2 in the United Kingdom for the period of 201 days from March 13 (n = 1) to September 29 (n = 201) gives the following approximating function (see Figure 5 ) .
The approximation errors are: MAD = 1, 443.37, SMAPE 1 = 0.05955 (for the shorter period of 182 days, beginning from April 1, SMAPE 1 = 0.007123), SMAPE 2 = 0.003047.
In the paper, we have proved, by using the Eulerian numbers, properties of the logistic function related to zeros of its successive derivatives. We also used logistic wavelets to estimate the parameters of a logistic curve that best fits to a given time series with a potential logistic trend. Then we described, based on the data from the United States and the United Kingdom, that the total reported number of SARS-CoV-2 infections can be modeled, in a natural way, as a sum of several logistic functions. The theory and the procedure can be applied to model the number of infections in any country or a region.
In our further work, we intend to use, in a similar way, the logistic wavelets of higher order (see Sec. 4). Using some appropriate special numbers we are going to define analogous wavelets for the Gompertz function (see some initial calculations [13] , [14] ) or for the fractional logistic functions (some preliminary theorems see [15] ).