Ooura, Takura

Kyoto University Research Information Repository

TitleAn IMT-type Quadrature Formula with the Same AsymptoticPerformance as the DE Formula(Mathematical Sciences forLarge Scale Numerical Simulations)Author(s)Ooura, TakuraCitation数理解析研究所講究録 (2007), 1573: 13-21Issue Date2007-11URL http://hdl.handle.net/2433/81326RightType Departmental Bulletin PaperTextversionpublisherKyoto UniversityAn IMT-type Quadrature Formula with the SameAsymptotic Performance as the DE Formula(Takuya Ooura)Research Institute for Mathematical Sciences,Kyoto University1 IntroductionIn the family of numerical quadrature formulas obtained by variable transformation, there are the IMT formula $[1, 2]$ and the DE formula [5]. The DE formula hasan asymptotic error estimate $O(\exp(-cN/\log N))$ , where $N$ is the number of thesampling points, and $c$ is a posItive constant independent of $N$ . On the other hand,the asymptotic error of the IMT formula behaves as $O(\exp(-c\sqrt{N}))$ . DE formula’sperformance is therefore better than IMT’s. Though the repeated application ofthe IMT-type transformation [4] is known to give a substantial improvement to theIMT formula, the asymptotic error of the IMT-Double formula including the IMT-type double exponential formula [3], the IMT-Triple formula, the IMT-Quadrupleformula, $\cdots$ behaves as $O(\exp(-cN/(\log N)^{2})),$ $O$ ($\exp(-cN/((\log N)$ (log log $N)^{2}))$ ),$O$ ($\exp(-cN/((\log N)$ (log log $N)(\log\log$ log $N)^{2}$ ))), $\cdots$ , which do not attain the per-formance of the DE formula still. We propose in the present paper an IMT-typequadrature fomula with the asymptotic error estimate $O(\exp(-cN/\log N))$ . Theidea of the proposed formula is to optimally choose the parameters of the IMT-typetransformation depending on the number $N$ of sampling points.2 . IMT FormulaFirst of all, we review the IMT formula. The IMT formula consists of a changeof variable and the trapezoidal rule with equal mesh size. The trapezoidal rulewith equal mesh size is known to be very efficient if applied in the following cases:1. integration of a periodic smooth function;2. integration of a smooth fmction over the entire interval $(-\infty, \infty),\cdot$The IMT formula uses a variable transformation of the first case. On the otherhand, the DE formula uses a variable transformation of the second case.1573 2007 13-21 13A typical IMT formula uses the following variable transformation:$x= \phi(t)=\frac{1}{Q}\int_{0}^{t}$ exp $(- \frac{1}{s}-\frac{1}{1-s})ds$ , $Q= \int_{0}^{1}\exp(-\frac{1}{s}-\frac{1}{1-s})ds$and transforms a given integral$I= \int_{0}^{1}f(x)dx$into$I= \int_{0}^{1}f(\phi(t))\phi’(\cdot t)dt$ ,whose integrand becomes a smooth periodic function. Applying the trapezoidalrule with equal mesh size $h$ results in the IMT formula:$I_{h}=h \sum_{n=1}^{N-1}f(\phi(nh))\phi’(nh)$ , $h=1/N$.The error of the IMT formula is asymptotically estimated [1, 2, 4] as$|I-I_{h}|=O(\exp(-c\sqrt{N}))$ ,where $c$ is independent of $N$.3 Proposed IMT-type QuadratureLet the given integral be$I= \int_{-1}^{1}f(x)dx$ . (1)An IMT-type transformation$x= \phi_{m,k}(t)=erf(\frac{k}{(1-t)^{m}}-\frac{k}{(1+t)^{m}}).$ ’ erf $x= \frac{2}{\sqrt{\pi}}\int_{0}^{x}e^{-t^{2}}dt$ (2)gives the following formula of IMT-type:$I_{h}=h \sum_{n=1}^{N-1}f(\phi_{m,k}(-1+nh))\phi_{m,k}’(-1+nh)$ , $h=2/N$, (3)where $m,$ $k$ are positive parameters which may depend on $N$ . The reason to usethe function erf $x$ is that we can easily compute the integral and that the errorestimation becomes easier. The feature in this proposal is to choose the optimalparameter depending on $N=2/h$.14In order to estimate the error, we assume that $f$ is analytic in some complexdomain which contains $(-1,1)$ . By the error analysis in Appendix 1, the errorterms of the quadrature (3) are given as $|I-I_{h}|\sim E_{1}+E_{2}$ ,$E_{1}$ $=A_{1}$ exp $(-( \alpha’k^{2})^{\frac{1}{2m+1}}(\frac{\pi N}{2m})^{\frac{2m}{2m+1}}\cdot(2m+1)\sin\frac{\pi/2}{2m+1})$ , (4)$E_{2}$ $=A_{2}$ exp $(- \frac{\pi^{3/2}\beta’N}{4mk}I,$ (5)where $\alpha’,$ $\beta^{j}$ are positive constants depending on singularities of $f$ , and $A_{1},$ $A_{2}$ arepositive constants depending on $f$ . We first choose the parameters so that $E_{1}$equals $E_{2}$ , and approximately optimize the order of $E_{1},$ $E_{2}$ in addition. Then, theparameters are$m= \frac{1}{2}$ log $N$, $k= \frac{e\beta’}{\sqrt{\pi}}$ (6)and the error terms are evaluated asymptotically as followslog $E_{1} \sim\log E_{2}\sim-\frac{\pi^{2}N}{2e\log N}$ , $Narrow\infty$ .In fact,$( \frac{\pi N}{2m})^{\frac{2m}{2m+1}}=\frac{\pi N}{2m}(\frac{2m}{\pi N})^{\frac{1}{2m+1}}=\frac{\pi N}{\log N}(\frac{2m}{\pi e^{2m}})^{\frac{1}{2m+1}}=\frac{\pi N}{\log N}(\frac{1}{e}+o(1))$ ,and$(2m+1)$ sin $\frac{\pi/2}{2m+1}=\frac{\pi}{2}+o(1)$ , $marrow\infty$ .Therefore, when the parameters are chosen as in (6), the error of our formula isestimated as$|I-I_{h}|=O(\exp(-cN/\log N))$ . (7)This order is the same as that of the DE Formula.4 Numerical ExamplesWe computed the following integrals$I_{1}$ $= \int_{-1}^{1}\sqrt{1-x^{2}}dx$ ,$I_{2}$ $= \int_{-1}^{1}\frac{dx}{1+x^{2}}$ ,$I_{3}$ $= \int_{-1}^{1}\log(1+x)dx$ ,15$I_{4}$ $= \int_{-1}^{1}\frac{dx}{(2+x)(1-x)^{3/4}(1+x)^{1/4}}$ ,$I_{5}$ $= \int_{-1}^{1}\frac{\cos\pi x}{\sqrt{1-x}}dx$ ,$I_{6}$ $= \int_{-1}^{1}\frac{dx}{\sqrt{100000001-x^{2}}}$by using the transformations:1. Our transformation: $\phi_{m,k}(t)=erf(k(1-t)^{-m}-k(1+t)^{-m}),$ $m=(1/2)\log N$ ,$k=2.2$ ;2. DE transformation: $\phi_{DE}(t)=t\bm{t}\bm{h}$( $(\pi/2)$ sinh $t$).We chose the parameter $k$ empirically: We computed $\int_{-1}^{q}(1-x^{2})^{\alpha}dx,$ $\alpha=0.25,0.5,0.75$with $k$ beginning $hom1.4$ to 3.0 with each 0.2. We found that $k=2.2$ commit-ted the least error. With this empirical fact, we chose $k=2.2$ throughout ourexperiments. The result is shown in Fig. 1, Fig. 2 and Fig. 3.$N_{*val}$ $N_{eva\prime}$Errors in the case of Errors in the case of$I_{1}= \int_{-1}^{1}\sqrt{1-x^{2}}dx$ $I_{2}= \int_{-1}^{1}\frac{dx}{1+x^{2}}$Fig. 1: Comparison between our formula and the DE formula for $I_{1}$ and $I_{2}$The horizontal axis is the number of sampling points actually computed bywhich we mean that we count only those $N$ sampling points such that $|f(\phi_{m,k}(-1+$$nh))\phi_{m,k}’(-1+nh)|$ is greater than $|I-I_{h}|$ . The vertical axis $represents-\log_{10}|I-$$I_{h}|$ . From these figures we see that our formula has almost as high degree of theperformance as the DE formula. In addition, maximal precisions obtained by ourformula are higher than those by the DE formula (the symbol $\triangle$ in the uppermostposition).16$N_{*\nu al}$ $N_{*\nu ael}$Errors in the case of Errors in the case of$I_{3}= \int_{-1}^{1}\log(1+x)dx$ $I_{4}= \int_{-1}^{1}\frac{dx}{(2+x)(1-x)^{3/4}(1+x)^{1/4}}$Fig. 2: Comparison between our formula and the DE formula for $I_{3}$ and $I_{4}$$N_{\nu\cdot/}$ $N_{\nu ol}$Errors in the case of Errors in the case of$I_{5}= \int_{-1\varpi-x}^{1coe\pi x}dx$ $I_{6}$ $\int_{-1}^{1}\frac{dx}{\sqrt 1.00000001-x^{2}}$Fig. 3: Comparison between our formula and the DE formula for $I_{5}$ and $I_{6}$17Appendix 1 Error AnalysisWe assume that $f(z)$ is analytic on a complex domain including the real interval$(-1,1)$ and has singularities which behave as1. $f(z)=a(1-z^{2})^{\alpha}+o(|1-z^{2}|^{\alpha})$ , $zarrow\pm 1,$ $\alpha>-1$ ;2.’$f(z)=b(z-*)^{-1}(z-\overline{*})^{-1}+O(1)$ , $zarrow z_{p}$ .We also assume that ${\rm Im} z_{p}>0$ and that zp is the singularity which is closest to$(-1,1)$ . The error of (3) is estimated [1, 2, 4] by$I-I_{h}=- \sum_{j=1}^{\infty}(-1)^{jN}(C_{jN}+C_{-jN})$ ,where $C_{n}$ are Fourier coefficients:$C_{n}= \int_{-1}^{1}g(x)e^{\pi inx}dx$ , $g(x)=f(\phi_{m,k}(x))\phi_{m,k}’(x)$ .Since $2|{\rm Re} C_{N}|$ is dominant in $|I-I_{h}|$ , it suffices for our purpose to calculate $C_{N}$ .The singularities of $g(z)$ are as follows:1. $g(z)\sim 4am(_{\sqrt{2\pi}^{k}})^{1-\alpha}(1\mp z)^{-m(1-\alpha)-1}\exp(-(1+\alpha)k^{2}((1-z)^{-m}-(1+z)^{-m})^{2})$,$zarrow\pm 1$ ;2. $g(z)\sim b(\phi_{m,k}(z)-z_{p})^{-1}(\phi_{m,k}(z)-\overline{z}_{p})^{-1}\phi_{m,k}’(z)$ , $\phi_{m,k}(z)arrow z_{p}$ ,and we change the path of integral in the following way:$C_{N}= \int_{-1}^{1}g(x)e^{\pi iNx}dx=\int_{\Gamma_{1}+\Gamma_{2}+\Gamma_{3}}g(z)e^{\pi iNz}dz$ .The path of $\Gamma_{1}+\Gamma_{2}+\Gamma_{3}$ is shown in Fig 4, where $\zeta_{p}$ is a pole transformed by$z_{p}=\phi_{m,k}(\zeta_{p})$ and is closest to the origin.lmFig. 4: The path of $\Gamma_{1},$ $\Gamma_{2},$ $\Gamma_{3}$18We consider the integral along $\Gamma_{1}$ . From the saddle point method [1, 2, 4], wehave$\log\int_{\Gamma_{1}}g(z)e^{\pi iNz}dz\sim-((1+\alpha)k^{2})^{\frac{1}{2m+1}}(2m+1)(\frac{\pi N}{2m}e^{-r;}2)^{\frac{2m}{2m+1}}-\pi iN$ .The estimation on $\Gamma_{3}$ is carried out in the same way and we havelog $\int_{\Gamma_{3^{-}}}g(z)e^{\pi iNz}dz\sim-((1+\alpha)k^{2})^{\frac{1}{2m+1}}(2m+1)(\frac{\pi N}{2m}e^{+\frac{\pi}{2}I})^{\frac{2m}{2m+1}}+\pi iN$ .By setting $1+\alpha=\alpha’+\epsilon,$ $\epsilon>0$ , we obtain (4).The integral along $\Gamma_{2}$ is estimated as$\int_{\Gamma_{2}}g(z)e^{\pi iNz}dz\sim\frac{\oint_{|z-z}\pi b}{{\rm Im} z_{p}}\exp(\pi iN\zeta_{p})p|=\epsilon f(z)e^{\pi iN\phi_{m,k}^{-1}(z)}dz$Since $\log(1\pm z)^{-m}=\mp mz+O(m|z|^{2})$ , we have$marrow\infty,$ $k=const.$ ,where $\eta_{p}$ is a point such that $z_{p}=erf\eta_{p}$ and that the quantity(8)is minimized. Here we assume that $|z_{p}|$ is small. Rom Appendix 2, such a min-imizer $\eta_{p}$ exists. We then define $\beta$ as the minimum value of (8). We see that$\beta\approx Im.z_{p},$ $z_{p}\approx(2/\sqrt{\pi})\eta_{p}\approx(4mk/\sqrt{\pi})\zeta_{p}$ . By setting$\beta=\beta’+\epsilon$ ,we obtain (5).Appendix 2 Zeros of erf $u$We consider the domain $\{x+iy|x\geq 0, y\geq 0\}$ , since $erf(-u)=$ -erf $u$ anderf $\overline{u}=\overline{erfu}$. erf $u$ has no zeros in the following domain$D=\{x+iy|x>0, y>0, (\sqrt{\pi}/2)e^{x^{2}}erfx>ye^{y^{2}}\}$ .19Fig. 5: Domain $D$ and zeros of erf $u$In fact, it holds there that$| \int_{0}^{x+iy}e^{-t^{2}}dt|$ $\geq$ $| \int_{0}^{x}e^{-t^{2}}dt|-|\int_{0}^{y}e^{-(x+it)^{2}}dt|$$\geq$ $\frac{\sqrt{\pi}}{2}erfx-ye^{-x^{2}+y^{2}}>0$.The domain $D$ and zeros of erf $u$ are shown in Fig. 5.If $|z_{p}|$ is assumed to be small enough, erf $u-z_{p}$ has zeros at$\eta_{0}=\frac{\sqrt{\pi}}{2}z_{p}+o(|z_{p}|)$ , $\eta_{j}=u_{j}+\frac{\sqrt{\pi}}{2}e^{u_{j}^{2}}z_{p}+o(|z_{p}|)$ , $j=1,2,3,$ $\cdots$ ,where $0,$ $u_{1}\approx 1.45+1\cdot.88i,$ $u_{2}\approx 2.24+2.61i,$ $\cdots$ are zeros of erf $u$ . Since theasymptote of $\partial D$ is $\{x+iy|x>0, y>0, x=y\}$ , we can estimate $|e^{u_{j}^{2}}|\leq M_{1}$ ,${\rm Im}$ log $\eta_{j}\geq M_{2}$ , andwhere $M_{1},$ $M_{2},$ $M_{3}$ are positive constants independent of $j=1,2,3,$ $\cdots$ . If 1 $z_{p}|$ issmall enough, then $\eta_{p}=\eta_{0}$ and $\beta<M_{3}$ .Acknowledgment. Helpful comments by Prof. H. Okamoto of Kyoto Universityand Prof. M. Sugihara of the University of Tokyo are acknowledged with gratitude.20Reference[1] M. Iri, S. Moriguti, Y. Takasawa, On a certain quadrature formula (inJapanese), Kokyuroku, RIMS Kyoto Univ. 91 (1970) 82-118.[2] M. Iri, S. Moriguti, Y. Takasawa, On a certain quadrature formula, J. Comput.Appl. Math. 17 (1987) 3-20 (English translation of [1]).[3] M. Mori, An IMT-type double exponential formula for numerical integration,Publ. RIMS, Kyoto Univ. 14 (1978) 713-729.[4] K. Murota, M. Iri, Parameter tuning and repeated application of the IMT-typetransformation in numerical quadrature, Numer. Math. 38 (1982) 347-363.[5] H. Takahasi, M. Mori, Double exponential formulas for numerical integration,Publ. RIMS, Kyoto Univ. 9 (1974) 721-741.21

An IMT-type Quadrature Formula with the Same Asymptotic Performance as the DE Formula(Mathematical Sciences for Large Scale Numerical Simulations)

http://repository.kulib.kyoto-u.ac.jp/dspace/bitstream/2433/81326/1/1573-02.pdf

An IMT-type Quadrature Formula with the Same Asymptotic Performance as the DE Formula(Mathematical Sciences for Large Scale Numerical Simulations)

Abstract

Similar works

Full text

Available Versions

Kyoto University Research Information Repository