{'fraction_non_alphanumeric': 0.052481994638762394, 'fraction_numerical': 0.021800213157639765, 'mean_word_length': 4.144376142216315, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Levin method for logarithmically singular oscillatory integrals

We propose a new stable Levin method to compute oscillatory integrals with logarithmic singularities and without stationary points. To avoid the singularity, we apply the technique of singularity separation and transform the singular ODE into two non-singular ODEs, which can be solved efficiently by the collocation method. Applying the equivalency of the new Levin method for the singular oscillatory integrals and the Filon method when the oscillator is linear, we consider the convergence of the new Levin method. This new method shares the proposition that less error for higher oscillation. Several numerical experiments are presented to validate the efficiency of the proposed method.


Introduction
We consider in this paper highly oscillatory integrals of Fourier type of the form where a is a positive real number, f and g are suitably smooth, g (x) = 0, x ∈ [0, a] and w is a real parameter whose absolute value can be extremely large. If the integral is over another bounded domain [a, b] with finitely many logarithmically singular points, it can be written in sums of integrals of the form I Oscillatory integrals with some logarithmic singularities occurs frequently in the numerical process of solving many problems of science and engineering such as electromagnetic and acoustic scattering. Since the antiderivatives of the integrands are unknown in most of cases, they have to be computed numerically. However, the high oscillation and the weak singularity of the integrands make the classic numerical integral methods such as Gauss quadrature hard to derive an acceptable approximation within a limited cost. The computation of integrals of this type I [0,a] log,w [f, g] is regarded as a challenging issue which requires special focus.
Many effective methods have been proposed for the oscillatory integrals in order to overcome the difficulty caused by the high oscillation such as Filon-type methods [10,11,27], Levin methods [15,20], the generalized quadrature rule [3], numerical steepest descent methods [9]. We refer the interested reader to [8] for a review of these methods. There is a series of papers that develop quadratures for non-oscillatory integrals with singularities by using graded meshes [12], or by Euler-Maclaurin summation formula [1,14,24,25]. An interesting hybrid Gauss-Trapezoidal quadrature rule was introduced in Alpert [2] for the integrand with algebraic or logarithmic singularity, and for improper integrals with oscillatory weight e iγx , where the quadrature nodes and weights are computed by solving a nonlinear system. It is high time-consuming especially for highly oscillatory integral since the minimum sampling was taken to be two points per period. The hybrid Gauss-trapezoidal rule is quite accurate for integrand without highly oscillation, but fails to computation of the highly oscillatory integral when the frequency is much bigger than the number of the nodes (see Table 1).   We next review the development of quadratures regarding both the properties of oscillation and logarithmic singularity. The asymptotic behavior of Fourier integrals involving logarithmic singularities was obtained by repeated integration by parts in [6,19]. This method is unstable. To overcome the difficulty of the singularity completely, a new Filon method was proposed in [23,30,4,13] for the case when the modified moments including the singularity can be obtained numerically for the linear oscillator g. However, it might be impossible for the general cases of the oscillator g since the modified moments are hard to be obtained.
Recently, by partition of the integration interval based on the singularity and the oscillation of the integrand, the composite Filon methods were considered in [5] and [18] based on the Filon-Clenshaw-Curtis quadrature and moment-free Filon method [27], respectively. The composite moment-free Filon-type methods developed in [5,18] are efficient for computing oscillatory integrals with weakly singular integrand and stationary points. Furthermore, the new methods in [18], unlike the existing Filon-type methods, do not have to compute the inverse of the oscillator, and have a polynomial order or exponential order of convergence. However, the main disadvantage of these composite methods is that subintervals near the singular point in the specified mesh have very small lengths and thus may cause serious round-off error problems.
A special Gauss-type quadrature, based on the numerical steepest method, has been proposed for the highly oscillatory integrals with algebraic singularities [7,31,32] for linear oscillators, not applied to general oscillators. There is still much work to compute the oscillatory integrals with logarithmic singularities in efficiency and accuracy.
Different from the existing methods, the purpose of this paper is to design an efficient quadrature rule based on the classic Levin method for the computation of integrals of the form I log,w [f, g]. The method developed in this paper requires no graded meshes and the computation of modified moments and can be easily extended to the case with a complicated oscillator g which inherit from the merits of the classic Levin method.
The Levin method, proposed in [15], is very efficient for computation of integrals of the form is not singular and g (x) = 0 for x ∈ [a, b], where to evaluate the integral is transformed into a certain ODE problem. In particular, one of the solutions of the ODE system is non-oscillatory and can be solved by a collocation technique. Compared with the other methods for the oscillatory integrals, the Levin method can be applied to a more general oscillator without explicit computation of the moments. It was also found that the Levin method is equivalent to the Filon method when applied to oscillatory integrals with a linear oscillator [28]. In addition, the Levin method can be implemented stably through the Chebyshev collocation method with TSVD [17] or by the GMRES method [22,21]. However, it cannot be applied directly to oscillatory integrals of the type I [0,a] log,w [f, g] with logarithmic singularity. For example, we compute by the classic Levin method a simple logarithmically singular integral where Si and Ci denote the sine and cosine integral functions, respectively. Due to the singularity of the integrand at x = 0, we use the modified Chebyshev-Gauss-Radau points t j = 1 + cos 2πj 2n−1 /2, j = 0, 1, . . . , n − 1 as the collocation points in the classic Levin method. The relative errors, shown in Table 2, reveal that the classic Levin method loses the spectral accuracy when w is small and it fails at all when the frequency is large enough. Our main idea of the new Levin method is the separation of the singularity from the solution of the singular ODE to avoid the influence from the singular forcing function. By using the technique of singularity separation, the solution of the singular ODE is transformed into solutions of two nonsingular ODEs based on the principle of superposition. The linear systems obtained in solving these two ODEs by the collocation method share the same matrix and thus they can be solved efficiently with a little increased cost compared with the classic Levin method. The new Levin method for oscillatory integrals with logarithmic singularities keeps nearly all of the merits of the classic Levin method for the non-singular oscillatory integrals: (1) It does not require the computation of moments and is applicable for the nonlinear oscillator case.
(2) It converges nearly superalgebraically with respect to the number of collocation points when f is smooth, g is linear and the ODEs are solved by the Chebyshev collocation method. Numerical experiments validate that it is also true for the general oscillator g without stationary points.
(3) It shares the important property: the higher of the frequency the more accuracy of the method. The asymptotic order is O(w −2 log(1 + |w|)) with respect to the frequency.
This paper is organized as follows. In section 2, we develop the new Levin method for integrals of the type I

A new Levin method for logarithmically singular oscillatory integrals
In this section, we propose a new Levin method for logarithmically singular oscillatory integrals with two algorithms for linear and nonlinear oscillators, respectively. We further assume that the oscillator satisfies g(0) = 0 and g (x) > 0, x ∈ [0, a], i.e. g is an increasing function starting from the origin. If g(0) = 0, then g(x) is replaced by g(x) − g(0) and if g (x) < 0, then g(x) is replaced by −g(x) and w by −w. For the notation simplification, log,w [f, g] is shorten as I log [f, g] in this section. The spirit of the Levin method for is based upon the fact that if f were of the form Thus it transforms the evaluation of the integral into an ODE problem, Lp(x) = f (x). It has been proven in [15,16] that the ODE possesses at least a non-oscillatory solution which inspires the Levin method, or the Levin collocation method. That is to find a polynomial p n (x) with degree ≤ n − 1 such that Lp n (c j ) = f (c j ), j = 0, 1, . . . , n − 1, (2.2) and to compute I The convergence rate of the Levin method has been studied extensively [16,28,21]. In addition, the Filon method and the Levin method were proved to be identical when the oscillator is a linear function in [28]. Furthermore, two numerically stable algorithms, the Levin-Chebyshev collocation method using TSVD and the GMRES-Levin collocation method are presented in [17] and [21], respectively.
To compute I log [f, g] according to the spirit of the Levin method, we hope to find a function p such that p(x)e iwg(x) = f (x) log(x)e iwg(x) , or (2.4) Then I log [f, g] = p(a)e iwg(a) − p(0)e iwg(0) . However, the direct application of the Levin-Chebyshev collocation method to the computation of I log [f, g] shall cause the great error. It is because that the logarithmically singular function log(x) leads to a singular Levin ODE which can not be solved efficiently by the collocation method based on polynomials. The main purpose of this paper is to address this difficulty and to propose efficient new Levin methods for singular integrals of this kind. The basic idea in solving the singular ODE of the Levin method is to separate the singularity and to transform the singular ODE to non-singular ODEs. Inspired by the equality (x(log x − 1)) = log x, the solution p should also possess the logarithmic singularity. With this property, we introduce a new form for the function p p(x) = q(x) log(x) + h(x). (2.5) By substituting (2.5) into (2.4), we obtain that To obtain a particular solution p, we consider naturally a particular case of (2.6), Once we get the particular solutions q and h of (2.7) and (2.8), respectively, we have a particular solution for (2.4). We next consider the system of ODEs (2.7) and (2.8). For this purpose, we recall a result in [16] about the existence of a 'non-oscillatory' solution of the ODE (2.7) which is understood as a solution whose many derivatives are bounded uniformly for w.
a] and the (2n + 1) derivatives of g are bounded uniformly in w, for |w| > 1. Then there exists a constant c independent of w and a solution q(x) of (2.7) satisfying where D is the differential operator and · denotes the maximum norm on [0, a].
It is clear from Lemma 2.1 that there exists a particular non-oscillatory solution satisfying the ODE (2.7) which can be approximated well by polynomials based on the collocation method no matter how large is the absolute value of w. However, if the solution obtained from (2.7) does not vanish at x = 0, equation (2.8) is still a strongly singular ODE. To loose this singularity of (2.8), we restrict q(0) = 0. It is proved in the following lemma that q(x) x possesses good regularity when q is smooth enough with the restriction q(0) = 0. Note that when x = 0, the value of q(x) x is defined by taking the limitation as x → 0.
Proof: When x = 0, it is deduced directly from the general Leibniz rule for the high derivative of a product of two factors that We next verify by induction on k that when x = 0, the k-th derivative ofq equals q (k+1) (0) k+1 for 0 ≤ k < n. When k = 0, it is obvious according to the definition ofq. Assuming that the case k − 1 is established for 1 ≤ k < n, i.e.q (k−1) (0) = q (k) (0) k , we then consider the k-derivative ofq at x = 0. According to the definition of the k-derivative, it is known that where the L'Hospital's rule was used. Again by the use of the L'Hospital's rule, it is easily validated that lim x→0q for k < n and the desired result follows by setting k = n − 1. 2 We thus obtain an initial problem , and q(0) = 0, (2.10) whose exact solution is represented by the highly oscillatory integral, i.e.
Instead of computing the integrals directly, we find a non-oscillatory and well-behaved particular solution q 1 satisfying the ODE (2.7) without the initial condition according to Lemma 2.1 and the expression for q, the solution of the initial problem (2.10), is then given by We turn to equation (2.8).
Note that when f is sufficiently smooth, equation (2.8) is non-singular. However, the solution q obtained from (2.10) is still highly oscillatory which would influence the solution of (2.8). To deal with the high oscillation separately, the solution q is rewritten in the form where q 2 is determined by It is known from Lemma 2.2 and its proof that q 2 is non-oscillatory and has good regularity if q 1 is non-oscillatory and behaves well. Especially, when q 1 is a polynomial, so is q 2 .
With the substitution of q in (2.8), we get that (2.14) By the linear superposition, the solution h(x) can be split into two parts h( where h j , j = 1, 2 satisfy, respectively, the equations According to Lemma 2.2 and the property of q 1 , the ODE (2.15) is a non-singular equation which can be solved efficiently by the classic Chebyshev collocation method with the TSVD. We next focus on equation (2.16) to present an analytic particular solution for (2.16) when g(x) = x. Since possessing an oscillatory forcing function, the equation (2.16) is very difficult to be solved numerically for general oscillator g due to the oscillatory properties from iwg and the oscillatory forcing function. To get around the obstacle, we turn to the special functions for help when g is linear and consider the case of general oscillators later. To this end, we introduce the complementary incomplete Gamma function and its property. The complementary incomplete Gamma function, denoted as Γ(α, z), is defined by and it has a known series expansion when α = 0, where γ is the Euler's constant, equaling approximately 0.57721566490153286060651 and Log(z) := log(|z|) + i arg(z) and arg(z) denotes the principle argument of z for a complex z.. We are now ready to derive a particular solution of (2.16) for the case with g(x) = x. It is well-known that the solution of (2.16) satisfying h 2 (0) = 0 has a closed form, Substituting the Taylor's series expansion of e iwt and using the expansion (2.17) of Γ(0, z), it is obtained that Thus a particular solution for (2.16) when g(x) = x is given by Combining the solutions of (2.7), (2.15) and (2.16), a particular solution p is derived successfully which reads The logarithmically singular and oscillatory integral with linear oscillator follows directly . We summarize the first algorithm of the new Levin method for logarithmically oscillatory integrals with g(x) = x as follows.
We next consider the second algorithm of the new Levin quadrature for the case with a general oscillator g(x). Since it is hard to derive h 2 from (2.16) numerically or analytically for the general oscillator, the Levin method for the linear case can not be applied directly. To overcome this difficulty, we split the integral I log [f, g] into two parts, By the Hospital's rule, there exists the limit Besides, the product function f (x) log x g(x) has good regularity if f and g is suitably smooth according to Lemma 2.2 since log x g(x) = − log g(x) x . It reveals that the integral I 1 (f ) is readily computed efficiently by the classic Levin method.
For the second integral I 2 (f ), we follow the same idea of singularity separation and pursue a particular solution p with the form p(x) = q(x) log(g(x)) + h(x) such that Similarly, it is obtained two ODEs for q(x) and h(x), The solution q can be represented by a sum of a non-oscillatory particular solution and a multiple of the general solution. To facilitate the solution of equation (2.25), the solution q(x) is formed as follows where q 1 is a function satisfying the ODE and and then a solution of h is the sum of h 1 and h 2 . With the help of the complementary incomplete Gamma function, a particular solution of ODE (2.29) is given explicitly by and h 2 (0) = 0. Since ODEs (2.26) and (2.28) are non-singular and possess at least one non-oscillatory solution, they can be solved efficiently by the Chebyshev-collocation methods. Once q 1 , h 1 and h 2 are obtained, the integral I 2 (f ) is readily computed according to the spirit of Levin idea. We conclude the algorithm of the new Levin method for integrals with a general oscillator g. Let

Algorithm 2.4 (Levin algorithm for a general oscillator)
Given a function f ∈ C 1 [0, a], g(x) ∈ C 1 [0, a], g (x) > 0 and g(0) = 0 and a positive integer n: where e j denotes a unit column vector of size n × 1 whose j-th element is 1 while the others 0.
In a word, the new Levin method adopts the separation of singularity and oscillation to get around the singular difficulty by transforming the singular ODE into three non-singular ODEs with two of whom can be solved efficiently by the Chebyshev-collocation methods and the other one is solved analytically. It makes the Levin idea applicable for oscillatory integrals with logarithmic singularities.

Error analysis
In this section, we present error analysis for the two algorithms of the new Levin method proposed in section 2.
We first reveal the relationship between the new Levin method and the Filon method when computing oscillatory integrals with a linear oscillator g(x) = x. To this end, denote the numerical solutions of q 1 , q 2 in (2.12) and h 1 in (2.15) obtained through Algorithm 2.3 based on arbitrary pointŝ x = {0 ≤x 0 <x 1 < . . . <x n−1 ≤ a} byq 1 ,q 2 andĥ 1 . It is obvious that Letĥ 2 denote the approximation of h 2 in (2.18) with q 1 (0) replaced byq 1 (0). The algorithm of the new Levin method for a 0 f (x) log(x)e iwx dx can be expressed as In the following, we present the relation between Q With the definition (2.18) ofĥ 2 , we have that Using the relation betweenq andq 2 , it is derived directly from (3.34) and (3.35) that Combining (3.32) and (3.36), it is obtained that where the second equality is assured by integration by parts.
Using the Fundamental Theorem of Algebra again, it has been proven in [28] that Since L(q 1 (0)e −iw· ) = 0, it is obtained from (3.38) and the equalityq(x) =q 1 (x) −q 1 (0)e −iwx that Lq =f . It has been proved in [28] that the Levin method is equivalent to the Filon method when calculating non-singular oscillatory integrals with the linear oscillator. Theorem 3.1 tells that this property is also kept for the new Levin method when dealing with the singular oscillatory integrals with the linear oscillator. Note that Theorem 3.1 is only true for the case when g(x) is linear which is needed in (3.34) and (3.38) in the proof. With this equivalence, the error analysis of the new Levin algorithms is readily obtained.
We next present the error analysis for Algorithm 2.3 with arbitrary pointsx when g(x) = x based on Theorem 3.1. Let E n (f ) := I For the purpose of bounding E n (f ), we recall two basic lemmas in the numerical analysis for the computation of oscillatory integrals. Lemma 3.2 (van der Corput-type lemma [26, p.332,334] [30]) Suppose that f ∈ C 1 [0, a] and g ∈ C 2 [0, a] satisfying |g (x)| ≥ 1, x ∈ [0, a] and g (x) is monotonic, then for all w > 0, there exists a constant C independent of w such that We are ready to analyze the absolute error E n (f ). Let C denote a generic constant independent of n and w whose value may be changed in each appearance.
Using the expression for interpolation errors, it is clear that where ξ 1 , ξ 2 ∈ [0, a] depending on the value of x. According to Lemma 3.2, there exists a constant C independent of n and w such that (3.43) Sincef (n) ≡ 0, the desired inequality (3.59) follows directly. When x 0 = 0 and x n−1 = a, there have ψ(0) = 0 and ψ(a) = 0. We can derive by integration by parts that Again using Lemma 3.2, we obtain whereψ(x) := ψ(x)/x. By Taylor's expansions, there exists ξ 3 ∈ (0, x) and ξ 4 ∈ (0, x) for a given x ∈ [0, a] such that It is derived by a direct computation that Combining the discussion above, we get that By Rolle's theorem, there exist z j ∈ (y j , y j+1 ) such that ψ (z j ) = 0, j = 0, 1, . . . , n − 3, and then we derive from the interpolation errors that where ξ 5 ∈ [0, a]. The bound for E n (f ) follows naturally. 2 It is obvious that when both end points are included in collocation points in the new Levin method, the asymptotic order is about O(w −2 (1 + log |w|)) according to Theorem 3.3.
We next give the error analysis for Algorithm 2.3 in which the collocation points are selected to be the modified Chebyshev-Gauss-Lobatto points. To this end, we recall the errors of f −f and its derivative wheref is the interpolant on [−1, 1] based on the Chebyshev points. Let · T be the Chebyshev-weighted 1-norm defined by , k > 4, (3.47) (ii) If f is analytic with |f (z)| ≤ M in the region bounded by the ellipse with foci ±1 and major and minor semiaxis lengths summing to ρ > 1, then there exists a constant C independent of n and ρ such that at−t 2 dt = V k < ∞ for some k ≥ 1 , then the numerical integral computed by Algorithm 2.3 satisfies for n ≥ k + 2 and k > 4, where C is a constant independent of n and w.
(ii) If f ((· + 1)a/2) is analytic with |f ((z + 1)a/2)| ≤ M in the region bounded by the ellipse with foci ±1 and major and minor semiaxis lengths summing to ρ > 1, then the numerical integral computed by Algorithm 2.3 satisfies for each n ≥ 1, where C is a constant independent of n, ρ and w.
Proof: Let F (t) := f ((t + 1)a/2), t ∈ [−1, 1] and the interpolant of F based on Chebyshev-Gauss-Lobatto points x j , j = 0, 1, . . . , n − 1 is denoted byF . It is easily obtained that The error bounds follows directly by combining the results of Lemma 3.4 and the inequality (3.44) in the proof of Theorem 3.3. 2 We easily conclude from Theorem 3.5 that the new Levin method possesses the quasi-superalgebaric convergence with respect to the number of collocation points for logarithmically singular and oscillatory integrals with a linear oscillator when the Chebshev points are adopted and f is analytic.
In the left of this section, we discuss the approximation error of Algorithm 2.4 for singular and oscillatory integrals with general oscillators. For this purpose, we first present an equivalent algorithm for Algorithm 2.4. Let Q . Proof: Applying a change of variables, y = g(x), to the second part of (2.22), we get that where (x) = x.

2
Theorem 3.7 tells that the new Levin algorithm for a general oscillator also possesses the same asymptotic order as that for linear oscillator. Numerical experiments in the later section will show that it also has the quasi-superalgebraic convergence for the general oscillator with respect to the number of collocation points when f is smooth.

Numerical experiments
In this section, we present four numerical experiments to verify the efficiency of the proposed new Levin method for singular oscillatory integrals. To achieve this goal, we also compare the computational performance of the proposed methods with that of the quadrature rules proposed in [4,18]. The numerical results presented below were all obtained by using Matlab 2017b on a laptop that has a Intel(R) Core(TM) i7-6500U CPU with 8GB of Ram memory. T m (x) log x 2 e iwx dx whose is the base for the Filon method in computation of the oscillatory integrals with logarithmic singularities. Their exact values can be computed by the recurrence relations [4].
These moments can be rewritten as 2 x dx which can be computed by the new Levin method, i.e. by Algorithm 2.3. Since T m is a polynomial of degree no more than m, the moment can be computed exactly by the new Levin method with m + 1 collocation nodes in theory. We present in Table 3 the absolute errors of the new Levin method in computing the integrals for m = [2 3 4 5 6] and w = [10 10 2 10 3 10 4 ] with n = m + 1 collocation nodes. It is found that the absolute errors of the moments computed by the new Levin method attaches the machine precision for different settings of m and w. It validates the prediction numerically and also shows the potential efficiency of the new method in computing the oscillatory integrals with logarithmically singularities.
where E 1 (z) := ∞ 1 e −zt tdt is an exponential integral. Since the classic Levin method for oscillatory integrals without singularities is well developed, we use the classic Levin method with 32 points to form the reference value for the second oscillatory integral with a nonlinear oscillator.
Numerical results of absolute errors are shown in Table 4 and Figure 1 for different values of w and n. For the dependence on n, it is shown that the absolute errors of the new Levin method decay drastically, as n increases slowly for fixed w = 10 2 and 10 5 no matter the oscillator is linear or nonlinear. It is consistent with the theory that the error possesses the quasi-superalgebraic convergence. As shown in Figure 1, the absolute errors scaled by w 2 1+log |w| are bounded for fixed n which validates numerically that the asymptotic order on the frequency is O(w −2 (1 + log |w|)) for both linear and nonlinear cases which matches well with the theoretical results in Theorems 3.3 and 3.5.  We next compare the performance of the new Levin method with that of the existing methods, the Filon-Clenshaw-Curtis method (FCC) in [4] and the composite moment-free Filon-type quadrature (CMFP) with an polynomial order of convergence in [18]. The reason of the choice of CMFP instead of the composite moment-free Filon-type quadrature (CMFE) with an exponential order of convergence is that the CMFE is not stable due to the use of the interpolation of large order in each segment as shown in our numerical experiments which are not presented here. We then recall the quadrature formulas of [4,18]. The FCC for the integral where p n is a polynomial of degree n − 1 which interpolates f at Chebyshev points cos jπ n−1 , j = 0, 1, . . . , n − 1. To introduce the CMFP method, we first simply review the (composite) moment-free Filon-type method [27] and the Gauss-Legendre quadrature rule. The moment-free Filon method approximates the integral b a f (x)e iwg(x) by where p n is a polynomial of degree n − 1 which interpolates (f /g ) • g −1 at g(t j ), j = 0, 1, . . . , m and t j , j = 0, 1, . . . , m are a set of distinguish points on [a, b]. The composite moment-free Filon-type rules used in CMFP reads The Gauss-Legendre quadrature rule for integral b a f (x)dx is given by where w j and t j are the standard weights and points of the Gauss-Legendre rule on the domain [−1, 1]. Suppose for a nonnegative integer r, the function g ∈ C r+1 [0, 1] has a single stationary point at zero and satisfies g (j) (0) = 0 for j = 1, . . . , r and g (r+1) (x) = 0 for x ∈ [0, 1]. Let σ r := g (r+1) ∞ /(r +1)!, w r := max{kσ r , k}, and λ r := w −1/(r+1) r . The CMFP method for integral 1 0 f (x)e ikg(x) dx is established by where x 0 = 0, x j = (j/s) p , p = (2m+1)/(1+µ), µ is the index of singularity of f , j = 1, 2, . . . , s, ϕ(x) = f (λ r x)e iwg(λrx) , y j = w (j−n)/n/(r+1) r , j = 0, 1, . . . , n, N j = q m/(m−1) j , q j = max{|g (y j−1 )|, |g (y j )|}x j−1 /g(x j−1 ), j = 1, 2, . . . , n and ν is the index of singularity of (f /g ) • g −1 . When f has only the logarithmic singularity, the value of µ is set to be 0.
which is also considered in [4]. The reference value of the integral is obtained by Mathematica with 50 digits. We present in Table 5 and 6 the relative errors and the CPU time for different values of w and of n by using different methods, respectively. To illustrate the dependence of each method on n, we introduce some notations. Setting f 1 (x) = 2 cos(4x) x 2 +x+1 , f 2 (x) = f 1 (−x) and n 1 = 2 n/2−3 , there exists 1 −1 cos(4x) For a given n in Table 5   x 2 +x+1 log x 2 e iwx dx (The settings for each method are as the same as those in Table 5 The settings for the CMFP is chosen according to those used in [18]. The results in Table 5 show that the accuracy of the proposed Levin method is comparable with that of the FCC and is better than the CMFP. For the CPU time, it is shown in Table 6 that the proposed Levin method outperforms the other two methods. Hence, the new Levin method is more efficient in computing oscillatory integrals with a linear oscillator. which is considered in [18]. The reference value of the integral is obtained by Mathematica with 50 digits. Since the moments are unknown, the FCC is not applicable in this example and then we compare only with the CMFP. Numerical results of the relative errors and the CPU time are shown in Table 7 and 8 for different values of w and of n by using different methods, respectively. Setting f (x) = 1 and g(x) = 2x + sin πx 2 /3, the proposed Levin method and the CMFP are implemented for a given n in Table 7 through Q [0,1],L log,w,n [f, g] and Q CM F P w,n 1 ,n 1 ,4,4 [f log(·), g], respectively, where n 1 = 2 n/2−1 . It is shown clearly that the new proposed method is more accurate than the CMFP and cost less computation time. Therefore, the new method is also more efficient in dealing with oscillatory integrals with a nonlinear oscillator.
Besides, extra numerical results show that when the number of n 1 increases up to 2 14 , the error of CMFP in Examples 4.3 and 4.4 starts to increase for w = 10 2 which is due to the round-off error of tiny meshes while the new method is free of the problem caused by tiny meshes.

Conclusions
We have constructed a numerically stable Levin method for computing highly oscillatory integrals with logarithmically singularities, which does not require knowledge of the derivatives of f . This method retains the most vital computational property: higher frequency requires less work. The proposed method possesses the asymptotic order with respect to w of O(w −2 (1 + log |w|)) and the quasi-superalgebraic convergence when f is analytic. As shown in the algorithms, it only needs the comparable computation cost of the classic Levin method.   In the future, we hope to generalize these results for computing oscillatory integrals with other singularities and stationary points.