Construction of conformal maps based on the locations of singularities for improving the double exponential formula

TL;DR

This paper introduces a new conformal mapping technique based on Schwarz-Christoffel transformations to enhance the double exponential formula's accuracy near singularities, validated through numerical experiments.

Contribution

It develops a generalized transformation formula using explicit Schwarz-Christoffel maps to improve the DE formula's handling of singularities.

Findings

Enhanced accuracy in numerical integration near singularities.

The new transformation generalizes existing DE transformations.

Numerical experiments confirm improved performance.

Abstract

The double exponential formula, or the DE formula, is a high-precision integration formula using a change of variables called a DE transformation; whereas there is a disadvantage that it is sensitive to singularities of an integrand near the real axis. To overcome this disadvantage, Slevinsky and Olver (SIAM J. Sci. Comput., 2015) attempted to improve it by constructing conformal maps based on the locations of singularities. Based on their ideas, we construct a new transformation formula. Our method employs special types of the Schwarz-Christoffel transformations for which we can derive their explicit form. Then, the new transformation formula can be regarded as a generalization of the DE transformations. We confirm its effectiveness by numerical experiments.

Tables

Table 1. Table 1. Examples of DE transformations which are written as ϕ ​ ( t ) = ψ ​ ( π 2 ​ sinh ⁡ ( t ) ) italic-ϕ 𝑡 𝜓 𝜋 2 𝑡 \phi(t)=\psi\!\left(\frac{\pi}{2}\sinh(t)\right) .

interval	integrand	$\varphi \left(t\right)$	$\psi$
$(-1,1)$	$f\left(x\right)$	$\tanh \left(\frac{\pi }{2}\sinh \left(t\right)\right)$	$\tanh (\cdot )$
$(-\infty ,\infty )$	$f\left(x\right)$	$\sinh \left(\frac{\pi }{2}\sinh \left(t\right)\right)$	$\sinh (\cdot )$
$(0,\infty )$	$f\left(x\right)$	$\exp \left(\frac{\pi }{2}\sinh \left(t\right)\right)$	$\exp (\cdot )$
$(0,\infty )$	$f_1\left(x\right)e^{-vx}\left(v>0\right)$	$\log \left(\exp \left(\frac{\pi }{2}\sinh \left(t\right)+1\right)\right)$	$\log \left(\exp (\cdot +1)\right)$

Table 2. Table 2. Parameters of Theorem   1 .

	DE	SO	New	New2
$\gamma$	1	1	1	1
$d$	0.346	$\pi /2$	$\pi /2$	?
${\beta }_2$	0.785	0.0695	0.252	0.252

Table 3. Table 3. Parameters of Theorem   1 .

	DE	SO	New	New2
$\gamma$	1	1	1	1
$d$	0.0976	$\pi /2$	$\pi /2$	?
${\beta }_2$	1.57	$5.77\times {10}^{-6}$	$5.12\times {10}^{-3}$	$5.12\times {10}^{-3}$

Table 4. Table 4. Parameters of Theorem   1 .

	DE	SO	New	New2
$\gamma$	1	1	1	1
$d$	0.155	?	$\pi /2$	?
${\beta }_2$	0.393	0.196	0.0377	0.0377

Table 5. Table 5. Parameters of Theorem   1 .

	DE	New	New2
$\gamma$	1	1	1
$d$	0.0139	$\pi /2$	?
${\beta }_2$	0.157	$1.85\times {10}^{-6}$	$1.85\times {10}^{-6}$

Equations

$${\int }_{-1}^{1}f(x)\mathrm{d}x\approx h\sum _{j=-n}^{n}f(\varphi (jh)){\varphi }^{\prime }(jh),$$

$$\varphi (t)=\tanh \left(\frac{\pi }{2}\sinh (t)\right).$$

$$f(\varphi (t)){\varphi }^{\prime }(t)=\mathrm{O}(\exp (-\beta {\mathrm{e}}^{|t|}))(t\to \pm \infty )$$

$${\mathcal{D}}_d=\{z\in \mathbb{C}\mid |\mathrm{Im}z|<d\}.$$

$$H^{\infty }({\mathcal{D}}_d,\omega )=\{f:{\mathcal{D}}_d\to \mathbb{C}\mid f(z)\text{is analytic in}{\mathcal{D}}_d,\text{and}||f||<\infty \},$$

$$||f||=\underset{z\in {\mathcal{D}}_d}{\sup }\left|f(z)/\omega (z)\right|.$$

$$|f(z)|\le ||f|||\omega (z)|(z\in {\mathcal{D}}_d)$$

$${\mathcal{E}}_{N,h}^{\mathrm{T}}(H^{\infty }({\mathcal{D}}_d,\omega ))=\underset{\begin{array}{c}f\in H^{\infty }({\mathcal{D}}_d,\omega )\\ ||f||\le 1\end{array}}{\sup }\left|{\int }_{-\infty }^{\infty }f(x)\mathrm{d}x-h\sum _{j=-n}^{n}f(jh)\right|.$$

$${\int }_{-d}^d|g(x+\mathrm{i}y)|\mathrm{d}y\to 0(x\to \pm \infty )$$

$$\underset{y\to d-0}{\lim }$$

$${\alpha }_1\exp (-{\beta }_1{\mathrm{e}}^{\gamma |t|})\le |\omega (t)|\le {\alpha }_2\exp (-{\beta }_2{\mathrm{e}}^{\gamma |t|}),$$

$${\mathcal{E}}_{N,h}^{\mathrm{T}}(H^{\infty }({\mathcal{D}}_d,\omega ))\le C_{d,\omega }\exp \left(-\frac{\pi d\gamma N}{\log (\pi d\gamma N/{\beta }_2)}\right),$$

$$h=\frac{\log (2\pi d\gamma n/{\beta }_2)}{\gamma n}.$$

$$\omega (t)=\mathrm{O}(\exp (-\beta \exp (\gamma |t|)))as|t|\to \infty ,t\in \mathbb{R},$$

$$\underset{H}{\max }$$

$$\displaystyle\quad\left\{\begin{aligned} &f(\psi(H(\cdot)))\psi^{\prime}(H(\cdot))H^{\prime}(\cdot)\in H^{\infty}(\mathscr{D}_{d},\omega)\\ &d>0\\ &\omega\ \text{satisfies the conditions of \lx@cref{creftypecap~refnum}{thm:Sugihara_DE}}\\ &d\gamma\leq\pi/2\ (\text{limitation of \lx@cref{creftypecap~refnum}{thm:Sugihara_nexists}}).\end{aligned}\right.$$

$$\frac{\pi d\gamma N}{\log (\pi d\gamma N/{\beta }_2)}\approx \frac{\pi d\gamma N}{\log N}(N\to \infty ),$$

$$\tilde{S}=\{{\tilde{\delta }}_j\pm ({\widetilde{ϵ}}_j+2k\pi )\mathrm{i}\mid j=1,\dots ,m,k\in \mathbb{Z}\},$$

$$s\notin H({\mathcal{D}}_{\pi /2})(s\in \tilde{S}\cup S_{\psi }).$$

$$H_{\mathrm{DE}}(t)=\frac{\pi }{2}\sinh (t).$$

$$H_{\mathrm{DE}}({\mathcal{D}}_{\pi /2})=\mathbb{C}\backslash \{y\mathrm{i}\mid |y|\ge \frac{\pi }{2}\}.$$

$$H_{\mathrm{SO}}(t)=C\sinh (t)+\sum _{k=1}^m{u}_k{t}^{k-1},$$

$$\max$$

$$H_{\mathrm{SO}}\left(x_j+\frac{\pi }{2}\mathrm{i}\right)={\tilde{\delta }}_j+{\widetilde{ϵ}}_j\mathrm{i}(j=1,\dots ,m).$$

$$H_{\mathrm{SC}}(z)=C{\int }_0^z\exp \left(\frac{1}{2}({\theta }_{+}-{\theta }_{-})\zeta \right)\left\{\prod _{j=1}^M{\cosh }^{{\alpha }_j-1}(\zeta -{\tau }_j)\right\}\mathrm{d}\zeta +D$$

$$P\cap (\tilde{S}\cup S_{\psi })=\varnothing .$$

$$\tilde{S}\cup S_{\psi }=\{{\tilde{\delta }}_j\pm ({\widetilde{ϵ}}_j+2k\pi )\mathrm{i}\mid j=1,\dots ,m,k\in \mathbb{Z}\},$$

$${\int }_0^t{\mathrm{e}}^{\Delta \theta \tau }\prod _{j=1}^M{\cosh }^{{\alpha }_j-1}(\tau -{\tau }_j)\mathrm{d}\tau$$

$$=$$

$${\int }_0^t\prod _{j=1}^M{\cosh }^{{\alpha }_j-1}(\tau -{\tau }_j)\mathrm{d}\tau =\frac{1}{2^{\bar{\theta }-1}{\bar{\theta }}^2}\sinh (\bar{\theta }t-\sum _{j-1}^M({\alpha }_j-1){\tau }_j)+\mathrm{O}(1)$$

$$H_{\mathrm{New}}(z)=C{\int }_0^z\frac{{\prod }_{j=1}^m\cosh (z-a_j)}{{\prod }_{j=1}^{m-1}\cosh (z-b_j)}\mathrm{d}z+D,$$

$$H_{\mathrm{New}}(t)$$

Full text

1. Introduction

The double exponential formula, or the DE formula, is a numerical integration formula using a change of variables called a DE transformation and the trapezoidal rule [9]. For example, an integral on the interval $(-1,1)$ is calculated as follows:

[TABLE]

where the DE transformation corresponding to this interval is

[TABLE]

DE transformations are changes of variables which make the transformed integrands decay double exponentially:

[TABLE]

for some $\beta>0$.

The advantages and disadvantages of the DE formula have been formulated by Sugihara [8]. He estimated them using a parameter $d$, the width of the domain around the real axis in which the transformed integrand is analytic. He has shown that the error of the DE formula converges on the order of $\mathrm{O}(\mathrm{e}^{-kN/\log N})$ as $N\to\infty$, where $N$ is the number of the nodes for the trapezoidal rule and $k$ is proportional to $d$. From this formulation, we see that the DE formula makes the error converges rapidly regardless of end-point singularities; whereas it has a disadvantage that is sensitive to singularities of the integrand near the real axis since they make the parameter $d$ small.

To overcome this problem, Slevinsky and Olver [7] proposed to improve the DE formula by modifying the DE transformations. They have revealed relations between the parameter $d$ and singularities, and proposed to make polynomial adjustments to the DE transformations based on the locations of singularities.

In this paper, we construct a new transformation formula based on their idea. First, we list the options of the transformation formulas using the idea of the Schwarz-Christoffel transformation. Then, we choose the optimal one from the perspective of the precision and ease of numerical integration. This transformation formula is not only a modification of the DE transformations, but also can be considered to be a generalization of them. We confirm its effectiveness by numerical experiments.

The rest of the paper is organized as follows. In Section 2, we summarize the Sugihara’s analysis. In Section 3, we describe the idea to improve the DE formula and introduce the method by Slevinsky and Olver. In Section 4, we show the proposed methods. In Section 5, we show numerical experiments. Finally, we conclude this paper in Section 6.

We show proofs and calculations omitted in this paper in the appendix. Programs for the proposed methods are available from [3].

2. Precision of the double exponential formula

On the basis of theorems in [8], we formulate the precision of the DE formula by evaluating the error of the trapezoidal formula in the case where the integrands decay double exponentially.

We define a family of integrands. Let $d$ be a positive number and let $\mathscr{D}_{d}$ denote the strip region of width $2d$ about the real axis:

[TABLE]

Let $\omega$ be a non-vanishing function defined on the region $\mathscr{D}_{d}$, and define the weighted Hardy space $H^{\infty}(\mathscr{D}_{d},\omega)$ by

[TABLE]

where the norm of $f$ is given by

[TABLE]

For the following discussions, we assume that the function $\omega$ decays double exponentially. Then, since

[TABLE]

holds from the definition of the norm eq. 6, the weighted Hardy space $H^{\infty}(\mathscr{D}_{d},\omega)$ represents the family of integrands which decay double exponentially.

Let $N=2n+1$ be the number of the nodes for numerical integration. The error of the $N$-point trapezoidal formula is estimated using an error norm. Let $\mathscr{E}_{N,h}^{\mathrm{T}}(H^{\infty}(\mathscr{D}_{d},\omega))$ denote the error norm in $H^{\infty}(\mathscr{D}_{d},\omega)$:

[TABLE]

The following theorem gives the upper bound of this error norm. Let $B(\mathscr{D}_{d})$ denote the family of functions $g$ which are analytic in $\mathscr{D}_{d}$ and satisfy

[TABLE]

and

[TABLE]

Theorem 1** (Sugihara [8]).**

Suppose that the function $\omega$ satisfies the following three conditions:

1. (1)

$\omega\in B(\mathscr{D}_{d})$; 2. (2)

$\omega$ does not vanish at any point in $\mathscr{D}_{d}$ and takes real values on the real axis; 3. (3)

the decay rate on the real axis of $\omega$ is specified by

[TABLE]

where $\alpha_{1}$, $\alpha_{2}$, $\beta_{1}$, $\beta_{2}$, $\gamma>0$.

Then the upper bound of the error norm is given as

[TABLE]

where $N=2n+1$, $C_{d,\omega}$ is a constant depending on $d$ and $\omega$, and the mesh size $h$ is chosen as

[TABLE]

This theorem shows that the error of the trapezoidal formula converges exponentially according the parameters $d,\gamma$ and $\beta_{2}$. The larger these parameters are, the better the convergence rate becomes. However, it has also been shown that the parameters have a restriction as shown by the following theorem.

Theorem 2** (Sugihara [8]).**

There exists no function $\omega$ that satisfies the following three conditions:

1. (1)

$\omega\in B(\mathscr{D}_{d})$; 2. (2)

$\omega$ does not vanith at any point in $\mathscr{D}_{d}$ and takes real values on the real axis; 3. (3)

the decay rate on the real axis of $\omega$ is specified by

[TABLE]

where $\beta>0$ and $\gamma>\pi/(2d)$.

DE transformations are changes of variables which make the transformed integrand $f(\phi(\cdot))\phi^{\prime}(\cdot)$ be a member of $H^{\infty}(\mathscr{D}_{d},\omega)$. However, they do not necessarily guarantee the optimal convergence rate of Theorem 1. In the following sections, we improve the DE formula by modifying the DE transformations so that the convergence rate will be better.

3. Improvement of the double exponential formula

We consider DE transformations which are written as $\phi(t)=\psi\!\left(\frac{\pi}{2}\sinh(t)\right)$ for some function $\psi$. We show examples of such DE transformations in Table 1, where the fourth was introduced in [5]. In these cases, $\psi$ are periodic with period $2\pi$ in the direction of the imaginary axis. We improve the DE formula by constructing a function $H$ and changing these transformations to $\phi(t)=\psi(H(t))$.

On the basis of theorems in Section 2, Slevinsky and Olver [7] proposed an idea to construct $H$ appropriately. In this section, we describe their idea and method from our point of view.

3.1. Construction of the transformation formula

We construct the transformation formula $H$ so that the transformed integrand $f(\psi(H(\cdot)))\psi^{\prime}(H(\cdot))H^{\prime}(\cdot)$ will be a member of $H^{\infty}(\mathscr{D}_{d},\omega)$. Then, the convergence rate of Theorem 1 is applied. We wish to maximize it under the limit of Theorem 2. Thus, we wish to determine the function $H$ according to the following optimization problem:

[TABLE]

However, it is difficult to solve this optimization problem generally. In order to make the problem more simply, we consider the asymptotic form of eq. 13 as $N\to\infty$. Since it is written asymptotically as

[TABLE]

the value of $d\gamma$ is dominant. Thus, we construct the function $H$ according to the following method:

• •

First, we restrict the options of $H$ so that $d\gamma=\pi/2$ will be satisfied. Here, we assume that $d=\pi/2$ and $\gamma=1$.

• •

Then, we choose $H$ from these options so that the parameter $\beta_{2}$ will be larger.

The condition $d=\pi/2$ is equivalent to the condition that the transformed integrand is analytic in $\mathscr{D}_{\pi/2}$. It is attained by avoiding singularities. We assume that $f$ has a finite number of singularities which are symmetric with respect to the real axis. We write these singularities as $S=\{\delta_{j}\pm\epsilon_{j}\mathrm{i}\mid j=1\dots m\}$. Let $\tilde{S}$ denote the preimage of $S$ by $\psi$. By the periodicity of $\psi$, we write elements of $\tilde{S}$ as

[TABLE]

where $\tilde{\delta}_{1}<\dots<\tilde{\delta}_{m}$ and $0<\tilde{\epsilon}_{j}\leq\pi\,(j=1,\ldots,m)$. Also, if $\psi$ has singularities, we write them as $S_{\psi}$. For example, $\psi=\tanh(\cdot)$ has singularities $S_{\psi}=\{(\pm\pi/2+2k\pi)\mathrm{i}\mid k\in\mathbb{Z}\}$. In order to make the transformed integrand analytic in $\mathscr{D}_{\pi/2}$, we need to make the image $H(\mathscr{D}_{\pi/2})$ avoid the singularities in $\tilde{S}$ and $S_{\psi}$; that is, we construct the function $H$ so that it will satisfy

[TABLE]

Figure 1 summarizes this condition.

The parameters $\gamma$ and $\beta_{2}$ appears in the coefficients of the transformation formula $H$. For example, if it is written as $H(t)=C\sinh(\gamma^{\prime}t)+\mathrm{o}(\mathrm{e}^{\gamma^{\prime}t})$ as $|t|\to\infty$, then we see that $\gamma^{\prime}=\gamma$ and that $C$ is proportional to $\beta_{2}$. Thus, we fix $\gamma^{\prime}$ to 1 and consider how to make $C$ larger under the condition eq. 17.

3.2. DE transformations

We review the DE transformations in the context of Section 3.1. We can write the function $H$ of the DE transformations by

[TABLE]

The image $H_{\mathrm{DE}}(\mathcal{D}_{\pi/2})$ is the entire complex plane with a pair of slits:

[TABLE]

DE transformations have advantages that they can avoid the singularities in $S_{\psi}$ and the parameter $C$ is large. However, they have a significant problem that they cannot avoid the singularities in $\tilde{S}$. Specifically, the parameter $d$ may be quite small if the integrand $f$ has singularities near the real axis.

3.3. Method of Slevinsky and Olver

Slevinsky and Olver [7] improved the DE formulas by making polynomial adjustments to them. They used the following formula as an option of the function $H$:

[TABLE]

where $C>0$ and $u_{1},\dots,u_{m}\in\mathbb{R}$. Then, they chose these parameters by solving an optimization problem. In it, they maximized $C$ under the condition that the image of the boundary $\partial\mathscr{D}_{\pi/2}$ under $H_{\mathrm{SO}}$ could pass through the singularities $\{\tilde{\delta}_{j}\pm\mathrm{i}\tilde{\epsilon}_{j}\}_{j=1}^{m}$. This is formulated as follows:

[TABLE]

Algorithms to solve this is shown in [6].

However, we find some points to be improved in their methods as follows:

• •

They does not consider the singularities $S_{\psi}$. Thus, the parameters $C$ and $d$ may be smaller than the optimal values.

• •

The limiting conditions of the optimization problem eq. 22 does not imply the condition of the singularities eq. 17. There are cases where $H_{\mathrm{SO}}$ is not an injection and $d<\pi/2$ even though the conditions eq. 22 are satisfied.

• •

There are cases where we cannot solve the optimization problem by [6] because of the difficulty of finding the solution numerically. In these cases, we cannot use their methods.

We show experiments which cause the second and third problem in Sections 5.3 and 5.4, respectively.

4. Proposed method

In this section, we propose a new method to construct the function $H$. In it, we use a conformal map from the domain $\mathscr{D}_{\pi/2}$ to a polygon $P$ as an option of $H$. Then, choosing the function $H$ corresponds to choosing the polygon $P$. It enables us to handle the condition of the singularities eq. 17 directly. Also, using this conformal map has an advantage that it is written explicitly using the idea of the Schwarz-Christoffel transformation.

4.1. Schwarz-Christoffel transformation

The conformal mapping $H$ is written explicitly using a modified version of the Schwarz-Christoffel transformation [2]. We assume that the polygon $P$ is symmetric with respect to the real axis and that it has vertices at $\pm\infty$. Let $P$ have $M$ pairs of the vertices other than $\pm\infty$ and let $\{\alpha_{j}\pi\}_{j=1}^{M}$ denote the interior angles of them. We assume that $0<\alpha_{j}\leq 2$ if the corresponding vertex is finite and that $-2\leq\alpha_{j}<0$ if it is infinite. Let $\theta_{+}$ and $\theta_{-}$ denote the divergence angles of $\pm\infty$. We assume that $0\leq\theta_{+},\theta_{-}\leq 1$. Figure 2 shows an example of the polygon $P$. Then, the conformal mapping from the domain $\mathscr{D}_{\pi/2}$ to the polygon $P$ is written as the following theorem.

Theorem 3**.**

For a given polygon $P$, there are some real numbers $\tau_{1}<\tau_{2}\cdots<\tau_{M}$ such that a mapping

[TABLE]

is a conformal mapping from the domain $\mathscr{D}_{\pi/2}$ to the polygon $P$ which satisfies $H_{\mathrm{SC}}(+\infty)=+\infty$ and $H_{\mathrm{SC}}(-\infty)=-\infty$.

Proof.

A conformal mapping from the strip region $\{\zeta\mid 0<\mathrm{Im}[\zeta]<1\}$ to the polygon $P$ has been shown in [2]. We obtain eq. 23 by transforming the domain linearly. ∎

The problem of how to choose $\tau_{1},\ldots,\tau_{m}$ is known as the Schwarz-Christoffel parameter problem, which has been studied in [2, 10].

From Theorem 3, the problem of how to choose the function $H$ is changed into the problem of how to choose the polygon $P$. We discuss it in the following subsection.

4.2. Suitable polygon for numerical integration

We consider how to construct the suitable polygon for numerical integration according to the discussions in Section 3.1. The condition of the singularities eq. 17 is written simply as

[TABLE]

We choose the polygon $P$ so that the corresponding parameters $\gamma$ and $C$ will be larger under this condition. For the following discussions, we rewrite

[TABLE]

where $\tilde{\delta}_{1}<\dots<\tilde{\delta}_{m}$ and $0<\tilde{\epsilon}_{j}\leq\pi\,(j=1,\ldots,m)$.

Relations between the polygon $P$ and the parameter $\gamma$ are obtained by asymptotic expansion of $H_{\mathrm{SC}}$.

Theorem 4**.**

Let $\{\alpha_{j}\}_{j=1}^{M},\theta_{+}$, and $\theta_{-}$ be the parameters introduced in Section 4.1. We write $\bar{\theta}=(\theta_{+}+\theta_{-})/2$ and $\Delta\theta=(\theta_{+}-\theta_{-})/2$. Then,

[TABLE]

holds as $|t|\to\infty,t\in\mathbb{R}$. Specifically, when the divergence angles of the both sides are equal, i.e., $\theta_{+}=\theta_{-}=\bar{\theta}$,

[TABLE]

holds as $|t|\to\infty,t\in\mathbb{R}$.

The proof is given in Appendix A.

From this theorem, we see that $\gamma=\min\{\theta_{+},\theta_{-}\}$. Thus, we fix $\theta_{+}$ and $\theta_{-}$ to 1.

Relations between the polygon $P$ and the parameter $C$ are rather complex. Although it is difficult to formulate them, it is observed experimentally that the larger the area of the polygon $P$ is, the larger the parameter $C$ is. We show the experiments in Appendix B.

For these reasons, we propose to construct a new transformation formula $H_{\mathrm{New}}$ so that the image $H_{\mathrm{New}}(\mathscr{D}_{\pi/2})$ will avoid singularities by $m$ pairs of slits as Figure 3. Let $P_{\mathrm{New}}$ denote this. The function $H_{\mathrm{New}}$ coincides with DE transformations if $\tilde{S}=\emptyset$ and $S_{\psi}=\{(\pm\pi/2+2k\pi)\mathrm{i}\mid k\in\mathbb{Z}\}$. It can be considered to be a generalization of DE transformations.

The corresponding transformation to $P_{\mathrm{New}}$ is written as follows. The polygon $P_{\mathrm{New}}$ has $(2m-1)$ pairs of vertices other than $\pm\infty$ at the turning points of the slits (with angles $2\pi$) and the points at infinity (with angles 0). Let $\{a_{j}\pm\frac{1}{2}\pi\mathrm{i}\}_{j=1}^{m}$ and $\{b_{j}\pm\frac{1}{2}\pi\mathrm{i}\}_{j=1}^{m-1}$ denote the preimages of these vertices under $H_{\mathrm{New}}$, respectively. Then, we can write $(\tau_{1},\ldots,\tau_{2m-1})=(a_{1},b_{1},a_{2},\ldots,b_{m-1},a_{m})$ and $(\alpha_{1},\ldots,\alpha_{2m-1})=(2,0,2,\ldots,0,2)$. The transformation $H_{\mathrm{New}}$ is obtained by substituting these parameters into eq. 23, that is,

[TABLE]

where $C>0$ and $a_{1}<b_{1}<\cdots<b_{m-1}<a_{m}$.

Using Theorem 4, the asymptotic form of $H_{\mathrm{New}}$ is written as

[TABLE]

where $T=a_{1}-b_{1}+a_{2}-\cdots-b_{m-1}+a_{m}$.

The parameters are determined as follows. First, we can choose the parameter $T$ arbitrarily. We determine it so that the parameter $\beta_{2}$ will be the largest. Then, we determine the other parameters $C,a_{1},\ldots,a_{m}$ and $b_{1},\ldots,b_{m-1}$ so that the image $H_{\mathrm{New}}(\mathscr{D}_{\pi/2})$ will match $P_{\mathrm{New}}$. It has known that these parameters are uniquely determined with the value of $T$ fixed [2]. We show how to calculate these parameters in Sections 4.3 and 4.4, respectively.

Incidentally, there is another advantage to construct the transformation $H_{\mathrm{New}}$ in this way. The formula $H_{\mathrm{New}}$ is written rather simply. This is important when we use $H_{\mathrm{New}}$ as a change of variables. We show this in Section 4.4.

4.3. Determination of the parameter T

Here, we consider an integral of the interval $(-1,1)$. The cases of the other intervals are shown in Appendix C.

We assume that the integrand $f$ is smooth on the interval $(-1,1)$ and satisfies

[TABLE]

for some $p,q>-1$. The change of variables is given as

[TABLE]

The decay rate of the transformed integrand is estimated as

[TABLE]

for arbitrary $\varepsilon>0$. Then we see that the parameter $\beta_{2}$ satisfies

[TABLE]

To make the parameter $\beta_{2}$ larger, we we make $\varepsilon$ go to 0 and determine $T$ as

[TABLE]

Then the supremum of the parameter $\beta_{2}$ is estimated as

[TABLE]

4.4. Determination of the other parameters

First, we discuss the case of $m=2$:

[TABLE]

The integrand of eq. 36 is rearranged as

[TABLE]

where $2L_{1}=\mathrm{e}^{-a_{1}+b_{1}-a_{2}}\!\left(\mathrm{e}^{2a_{1}}+\mathrm{e}^{2a_{2}}-\mathrm{e}^{2b_{1}}-\mathrm{e}^{2\left(a_{1}-b_{1}+a_{2}\right)}\right)$. Then, the transformation $H_{\mathrm{New}}$ is written as

[TABLE]

where $T=a_{1}-b_{1}+a_{2}$, $2D_{1}=CL_{1}$, and $D_{0}=D-CL_{1}\tan^{-1}\!\left(\mathrm{e}^{-b_{1}}\right)$. The value of the parameter $T$ has been determined in Section 4.3. We determine the other parameters so that the image of the upper boundary of the strip region $\mathscr{D}_{\pi/2}$ will match the upper slits of Figure 3. For this reason, we consider the image of $z=x+\frac{\pi}{2}\mathrm{i},x\in\mathbb{R}$:

[TABLE]

The real part of eq. 41 is given by

[TABLE]

Thus, we determine $D_{0}$ and $D_{1}$ as

[TABLE]

The imaginary part of eq. 41 is given by

[TABLE]

The function eq. 44 has local minimum points in $(-\infty,b_{1})$ and $(b_{1},\infty)$, which correspond to the parameters $a_{1}$ and $a_{2}$, respectively. The function values at them correspond to $\tilde{\epsilon}_{1}$ and $\tilde{\epsilon}_{2}$. Thus, we determine the parameters $C,a_{1},b_{1}$, and $a_{2}$ by solving

[TABLE]

under the constraints $C>0$ and $a_{1}<b_{1}<a_{2}$. Here, the condition $T=a_{1}-b_{1}+a_{2}$ is automatically satisfied by solving eq. 45. We show it later in the general case.

Then, we extend the discussions to the general case. The following proposition shows that we can deform the integrand similarly to the case of $m=2$.

Proposition 1**.**

We write $T=a_{1}-b_{1}+\cdots-b_{m-1}+a_{m}$. Then,

[TABLE]

holds for some $L_{1},\ldots,L_{m-1}\in\mathbb{R}$.

Proof.

We define $m$th degree polynomials $F_{1}$ and $F_{2}$ as $F_{1}(Z)=\prod_{j=1}^{m}(Z+\mathrm{e}^{2a_{j}})$ and $F_{2}(Z)=(Z+\mathrm{e}^{2T})\prod_{j=1}^{m-1}(Z+\mathrm{e}^{2b_{j}})$. Since the leading terms and constants of $F_{1}$ and $F_{2}$ coincide, We can write

[TABLE]

for some polynomial $G$ of which the degree is $(m-2)$ or less. Also, using Lagrange polynomials, the polynomial $G$ is written as

[TABLE]

Then, we obtain eq. 46 by

[TABLE]

where $L_{j}=\mathrm{e}^{-T}l_{j}/2$. ∎

Using this proposition, we can rewrite $H_{\mathrm{New}}$ similarly as

[TABLE]

for some $D_{0},D_{1},\ldots,D_{m}\in\mathbb{R}$.

The parameters are also determined similarly to the case of $m=2$. The parameter $T$ has been determined in Section 4.3. The parameters $D_{0},D_{1},\ldots,D_{m-1}$ are determined based on the real part of $H_{\mathrm{New}}\!\left(x+\frac{\pi}{2}\mathrm{i}\right)$:

[TABLE]

The other parameters are determined based on the imaginary part of $H_{\mathrm{New}}(x+(\pi/2)\mathrm{i})$, that is, the parameters $C$, $a_{1},\ldots,a_{m}$ and $b_{1},\ldots,b_{m-1}$ are determined by solving

[TABLE]

under the constraints $C>0$ and $a_{1}<b_{1}<\cdots<b_{m-1}<a_{m}$. Here, the condition $T=a_{1}-b_{1}+\cdots-b_{m-1}+a_{m}$ is automatically satisfied by solving eq. 54. The following theorem shows this.

Theorem 5**.**

We assume that a system of equations

[TABLE]

holds for some real numbers $a_{1},\ldots,a_{m}$, $b_{1},\ldots,b_{m-1}$, $D_{0},\ldots,D_{m-1}$, C and T which satisfy $a_{1}<b_{1}<\cdots<b_{m-1}<a_{m}$ and $(C,D_{1},\ldots,D_{m})\neq(0,\ldots,0)$. Then,

[TABLE]

holds.

The proof is given in Appendix D.

In general, it is difficult to solve a non-linear system of equations such as eq. 54 numerically under constraints. Thus, we replace the parameters with $x\in\mathbb{R}^{2m}$ as

[TABLE]

and we solve eq. 54 as a system of equations of $x$. This is a similar method to numerical computations of Schwarz-Christoffel transformation [10]. We solve it using NLsolve, a Julia program for solving non-linear systems of equations [4].

4.5. Approximation of the proposed transformation

The proposed transformation $H_{\mathrm{New}}$ has a disadvantage that calculating the terms of $\tan^{-1}$ takes a lot of time. For this reason, we also propose approximating them by

[TABLE]

where the values and the derivatives at $t=0$, and the limits as $t\to\pm\infty$ of both sides coincide. We construct an approximation formula $H_{\mathrm{New2}}$ by approximating $H_{\mathrm{New}}$ using eq. 58:

[TABLE]

where the parameters are determined by the methods of Subsection Sections 4.3 and 4.4.

5. Numerical experiments

We compare the effectiveness of the transformation formulas $H$ by some numerical experiments. In Sections 5.1 and 5.2, we deal with the same examples as those in [7]. In Section 5.3, we show an example that the transformation $H_{\mathrm{SO}}$ is not an injection. In Section 5.4, we show an example to which we cannot apply the method of Slevinsky and Olver.

5.1. Integral on a finite interval

We consider an integral on a finite interval [7]:

[TABLE]

where $\delta_{1}\pm\epsilon_{1}\mathrm{i}=-0.5\pm\mathrm{i}$ and $\delta_{2}\pm\epsilon_{2}\mathrm{i}=-0.5\pm 0.5\mathrm{i}$. The change of variables for this integral is $\phi(t)=\tanh(H(t))$. There are singularities as follows:

[TABLE]

First, the formulas $H$ are given by

[TABLE]

Figure 4 shows the images $H(\mathscr{D}_{\pi/2})$.

Then, we compare the performances of the formulas $H$ as transformation formulas for integration. Table 2 shows the parameters of Theorem 1, where the parameter $d$ of DE is calculated as

[TABLE]

Figure 5 shows the original and transformed integrands.

Finally, we compare the errors of the numerical integration. Figure 6(a) shows relations between orders $n$ and the errors. Figure 6(b) shows relations between time for the calculation of numerical integration and the errors. Here, we assume that $d_{\mathrm{New2}}=\pi/2$ when we calculate the mesh size of the trapezoidal formula eq. 12.

5.2. Integral on the infinite interval

We consider an integral on the infinite interval [7]:

[TABLE]

where $\delta_{1}\pm\epsilon_{1}\mathrm{i}=-2\pm\mathrm{i}$, $\delta_{2}\pm\epsilon_{2}\mathrm{i}=-1\pm 0.5\mathrm{i}$, $\delta_{3}\pm\epsilon_{3}\mathrm{i}=1\pm 0.25\mathrm{i}$, and $\delta_{4}\pm\epsilon_{4}\mathrm{i}=2\pm\mathrm{i}$. The change of variables for this integral is $\phi(t)=\sinh(H(t))$. There are singularities as follows:

[TABLE]

First, the formulas $H$ are given by

[TABLE]

Figure 7 shows the images $H(\mathscr{D}_{\pi/2})$.

Then, we compare the performances of the formulas $H$ as transformation formulas for integration. Table 3 shows the parameters of Theorem 1, where the parameter $d$ of DE is calculated as

[TABLE]

Figure 8 shows the original and transformed integrands.

Finally, we compare the errors of the numerical integration. Figure 9(a) shows relations between orders $n$ and the errors. Figure 9(b) shows relations between time for the calculation of numerical integration and the errors. Here, we assume that $d_{\mathrm{New2}}=\pi/2$ when we calculate the mesh size of the trapezoidal formula eq. 12.

5.3. Integral on a semi-infinite interval (i)

We consider an integral on a semi-infinite interval:

[TABLE]

where $\delta_{1}\pm\epsilon_{1}\mathrm{i}=0.3\pm 0.2\mathrm{i}$, $\delta_{2}\pm\epsilon_{2}\mathrm{i}=0.5\pm 0.6\mathrm{i}$, $\delta_{3}\pm\epsilon_{3}\mathrm{i}=0.8\pm 0.5\mathrm{i}$, and $\delta_{4}\pm\epsilon_{4}\mathrm{i}=1.2\pm 0.3\mathrm{i}$. The change of variables for this integral is $\phi(t)=\sinh(H(t))$. There are singularities as follows:

[TABLE]

First, the formulas $H$ are given by

[TABLE]

Figure 10 shows the images $H(\mathscr{D}_{\pi/2})$.

Then, we compare the performances of the formulas $H$ as transformation formulas for integration. Table 4 shows the parameters of Theorem 1, where the parameter $d$ of DE is calculated as

[TABLE]

Figure 11 shows the original and transformed integrands.

Finally, we compare the errors of the numerical integration. Figure 12(a) shows relations between orders $n$ and the errors. Figure 12(b) shows relations between time for the calculation of numerical integration and the errors. Here, we assume that $d_{\mathrm{SO}},d_{\mathrm{New2}}=\pi/2$ when we calculate the mesh size of the trapezoidal formula eq. 12.

5.4. Integral on a semi-infinite interval (ii)

We consider an integral on a semi-infinite interval:

[TABLE]

where $\delta_{1}\pm\epsilon_{1}\mathrm{i}=1\pm 0.1\mathrm{i}$, $\delta_{2}\pm\epsilon_{2}\mathrm{i}=2\pm 0.5\mathrm{i}$, $\delta_{3}\pm\epsilon_{3}\mathrm{i}=3\pm 0.3\mathrm{i}$, $\delta_{4}\pm\epsilon_{4}\mathrm{i}=4\pm 0.5\mathrm{i}$, $\delta_{5}\pm\epsilon_{5}\mathrm{i}=5\pm 0.2\mathrm{i}$, $\delta_{6}\pm\epsilon_{6}\mathrm{i}=6\pm 0.5\mathrm{i}$, $\delta_{7}\pm\epsilon_{7}\mathrm{i}=7\pm 0.1\mathrm{i}$, $c_{2}=0.8$, $c_{3}=0.2$, $c_{4}=0.5$, $c_{5}=0.1$, and $c_{6}=0.5$. The change of variables for this integral is $\phi(t)=\log(\exp(H(t))+1)$. There are singularities as follows:

[TABLE]

First, the formulas $H$ are given by

[TABLE]

where the method of Slevinsky and Olver is omitted because we could not solve the optimization problem by their program [6]. Figure 13 shows the images $H(\mathscr{D}_{\pi/2})$.

Then, we compare the performances of the formulas $H$ as transformation formulas for integration. Table 5 shows the parameters of Theorem 1, where the parameter $d$ of DE is calculated as

[TABLE]

Figure 14 shows the original and transformed integrands.

Finally, we compare the errors of the numerical integration. Figure 15(a) shows relations between orders $n$ and the errors. Figure 15(b) shows relations between time for the calculation of numerical integration and the errors. Here, we assume that $d_{\mathrm{New2}}=\pi/2$ when we calculate the mesh size of the trapezoidal formula eq. 12.

6. Conclusion

We improved the DE formula in the case where the integrand has finite singularities. We improved it by constructing new transformation formulas $H_{\mathrm{New}}$ and $H_{\mathrm{New2}}$. The transformation $H_{\mathrm{New}}$ could be considered to be a generalization of the DE transformations and $H_{\mathrm{New2}}$ was an approximation of it. By numerical experiments, we confirmed the effectiveness of them even when the methods of the previous research failed.

In future work, we will need to consider cases where we do not know the locations of singularities or integrands have infinite singularities.

Appendix A Proof of Theorem 4

We define $I,J,I_{0}$ and $J_{0}$ as

[TABLE]

First, we consider asymptotic expansion of $I$ and $J$. We rearrange them as

[TABLE]

and

[TABLE]

Some of the terms in eq. 94 and eq. 98 can be ignored because they are on the order of $\mathrm{O}(1)$ as $|t|\to\infty$. Indeed, the polygon $P$ in Section 4.1 has $(2M+2)$ vertices including $\pm\infty$. Then, we see

[TABLE]

and specifically,

[TABLE]

Since $\theta_{+}$ and $\theta_{-}$ satisfy $0\leq\theta_{+},\,\theta_{-}\leq 1$, we see that

[TABLE]

We can show $I(\alpha_{1},\ldots,\alpha_{M-1},\alpha_{M}-2)$, $J(\alpha_{1},\ldots,\alpha_{j-1},\alpha_{j}-1,\alpha_{j+1},\ldots,\alpha_{M-1},\alpha_{M}-1)$, and $J(\alpha_{1},\ldots,\alpha_{M-1},\alpha_{M}-2)=\mathrm{O}(1)\ (|t|\to\infty)$ in a similar manner.

Furthermore, since the first term of eq. 94 is rearranged as

[TABLE]

as $t\to+\infty$ and

[TABLE]

as $t\to-\infty$, we can write

[TABLE]

as $|t|\to\infty$. Similarly, the first term of eq. 98 is written as

[TABLE]

By solving eq. 94 and eq. 98 with respect to $I$ and using eq. 104 and eq. 105, we obtain

[TABLE]

Appendix B Examination with respect to the parameter C

In this section, we show relations between polygons $P$ and the corresponding parameters $C$ experimentally. We calculate the parameter $C$ using sc-toolbox, a MATLAB program to solve the Schwarz-Christoffel parameter problem [1].

Let $\tilde{\delta}_{1}\pm\tilde{\epsilon}_{1}\mathrm{i}=-2\pm\mathrm{i}$ and $\tilde{\delta}_{2}\pm\tilde{\epsilon}_{2}\mathrm{i}=2\pm\mathrm{i}$ be singularities which the polygon $P$ need to avoid. Then, let $\eta$ be a positive number and we consider the following two kinds of polygons $P$:

1. (a)

A polygon which connects the singularities and 4 vertices $(\pm 1,\pm\eta)$. 2. (b)

A polygon which avoids the singularities by slits with width $2\eta$.

We show these polygons in Figure 16.

Figure 17 shows relations between the parameters $\eta$ and $C$. We see that the larger the area of $P$ is, the larger the parameter $C$ is.

Appendix C Determination of the parameter T for the other intervals

In the paper, we show how to determine the parameter $T$ in the case where the interval is $(-1,1)$. In this section, we show the other cases.

C.1. Infinite interval

We consider an integral of the interval $(-\infty,\infty)$. We assume that the integrand $f$ is smooth on the interval $(-\infty,\infty)$ and satisfies

[TABLE]

for some $r,s<-1$. The change of variables is given as

[TABLE]

The decay rate of the transformed integrand is estimated as

[TABLE]

for arbitrary $\varepsilon>0$. Then we see that the parameter $\beta_{2}$ satisfies

[TABLE]

To make the parameter $\beta_{2}$ larger, we we make $\varepsilon$ go to 0 and determine $T$ as

[TABLE]

Then the supremum of the parameter $\beta_{2}$ is estimated as

[TABLE]

C.2. Semi-infinite interval (i)

We consider an integral of the interval $(0,\infty)$. We assume that the integrand $f$ is smooth on the interval $(0,\infty)$ and satisfies

[TABLE]

for some $r<-1$ and $q>-1$. The change of variables is given as

[TABLE]

The decay rate of the transformed integrand is estimated as

[TABLE]

for arbitrary $\varepsilon>0$. Then we see that the parameter $\beta_{2}$ satisfies

[TABLE]

To make the parameter $\beta_{2}$ larger, we we make $\varepsilon$ go to 0 and determine $T$ as

[TABLE]

Then the supremum of the parameter $\beta_{2}$ is estimated as

[TABLE]

C.3. Semi-infinite interval (ii)

We consider an integral of the interval $(0,\infty)$. We assume that the integrand $f$ is smooth on the interval $(0,\infty)$ and satisfies

[TABLE]

for some $q>-1,v>0$, and arbitral $\varepsilon>0$. The change of variables is given as

[TABLE]

The decay rate of the transformed integrand is estimated as

[TABLE]

for arbitrary $\varepsilon>0$. Then we see that the parameter $\beta_{2}$ satisfies

[TABLE]

To make the parameter $\beta_{2}$ larger, we we make $\varepsilon$ go to 0 and determine $T$ as

[TABLE]

Then the supremum of the parameter $\beta_{2}$ is estimated as

[TABLE]

Appendix D Proof of Theorem 5

For simplicity, we show the case of $m=4$.

We rearrange the given system of equations as follows:

[TABLE]

where $A_{1}=\mathrm{e}^{2a_{1}},A_{2}=\mathrm{e}^{2a_{2}},A_{3}=\mathrm{e}^{2a_{3}},A_{4}=\mathrm{e}^{2a_{4}},B_{1}=\mathrm{e}^{2b_{1}},B_{2}=\mathrm{e}^{2b_{2}},B_{3}=\mathrm{e}^{2b_{3}},T^{\prime}=\mathrm{e}^{2T},C^{\prime}=\frac{1}{16}\mathrm{e}^{-T-b_{1}-b_{2}-b_{3}},D^{\prime}_{1}=\frac{1}{4}\mathrm{e}^{-b_{2}-b_{3}},D^{\prime}_{2}=\frac{1}{4}\mathrm{e}^{-b_{1}-b_{3}}$, and $D^{\prime}_{3}=\frac{1}{4}\mathrm{e}^{-b_{1}-b_{2}}$. The equations eq. 128 can be seen as a linear system of equations with respect to $C^{\prime},D^{\prime}_{1},D^{\prime}_{2}$, and $D^{\prime}_{3}$. Since it has non-trivial solutions, we see that $\det X=0$, where

[TABLE]

Here, we define $X_{0}$ as

[TABLE]

Then, from the properties of the determinant, we see that

[TABLE]

Also, since $\det X_{0}$ is a polynomial of degree 9 $(=(m-1)^{2})$ and is divisible by $(A_{i}-A_{j})$ and $(B_{i}-B_{j})$ for arbitral $i<j$, we can write

[TABLE]

for some real number $x_{0}$.

We show that $x_{0}\neq 0$ as follows. Letting $B_{1}=A_{2}$, $B_{2}=A_{3}$, and $B_{3}=A_{4}$ formally, the matrix $X_{0}$ is an upper triangular matrix of which the diagonal components are $(A_{1}-A_{2})(A_{1}-A_{3})(A_{1}-A_{4}),(A_{2}-A_{3})(A_{2}-A_{4}),(A_{3}-A_{2})(A_{3}-A_{4})$, and $(A_{4}-A_{2})(A_{4}-A_{3})$. From this reason, we see that $\det X_{0}\not\equiv 0$, which implies $x_{0}\neq 0$.

Therefore, from eq. 131, we see that

[TABLE]

which implies that $a_{1}-b_{1}+\dots-b_{3}+a_{4}=T$.