数理解析研究所講究録
第 1566 巻 2007 年 64-78

64

NUMERICAL REAL INVERSIONFORMULAS OF
THE LAPLACE TRANSFORM BY USINGTHE SINC FUNCTIONS
T. MATSUURA (松浦勉), A. AL-SHUAIBI, H. FUJIWARA(藤原宏志), S. SAITOH (齋藤三郎) and M. SUGIHARA(杉原正顕)

Abstract
We shall give a very natural and numerical real inversion formulaof the Laplace transform for generaldata following the ideas of best
approximations, generalized inverses, the Tikhonov regularization andthe theory of reproducing kernels. Furthermore, we shall additionally
use the Sinc functions (the Sinc method) to our general theory, tosolve the related integral equation. However, the new method in thispaper for the real inversion formula will be reduced to the solutionof linear simultaneous equations. This real inversion formula may
be expected to be practical to calculate the inverses of the Laplacetransform by computers when the real data contain noise or errors.
We shall illustrate examples and justify our computational work.
$L_{2}$

1

Introduction

We shall give a very natural and numerical real inversion formula of
the Laplace transform$\underline{(\mathcal{L}F)(p)=f(}p)=\int_{0}^{\infty}e^{-pt}F(t)dt$

,

$P>0$

(1.1)

’2000 Mathematics Subject Classification: Primary $44A15,35K05$ , 30C40 Key words
and phrases: Laplace transform, real inversion, numerical inversion, reproducing kernel,
Tikhonov regularization, Sobolev space, approximate inverse, sinc method, sampling the-

orem

65
for functions $F$ of some natural function space. This integral transformis, of course, very fundamental in mathematical science. The inversionof the Laplace transform is, in general, given by a complex form,
however, we are interested in and are requested to obtain its real
inversion in many practical problems. However, the real inversion
will be very involved and one might think that its real inversion willbe essentialy involved, because we must catch “analyticity” from thereal or discrete data. Note that the image functions of the Laplacetransform are analytic on some half complex plane. For complexityof the real inversion formula of the Laplace transform, we recall, forexample, the following formulas:$\lim_{narrow\infty}\frac{(-1)^{n}}{n!}(\frac{n}{t})^{n+1}f^{(n)}(\frac{n}{t})=F(t)$(Post [15] and Widder [27,28]), and$\lim_{narrow\infty}\Pi_{k=1}^{n}(1+\frac{t}{k}\frac{d}{dt})[\frac{n}{t}f(\frac{n}{t})]=F(t)$,([27,28]).

Furthermore, see [1-8,16,17,21,25,26,28,29] and the recent related
articles [10] and [11]. See ako the great references $[29,30]$ . The problem may be related to analytic extension problems, see [11] and [22].
In this paper, we shall give a new type and very natural real inversion
formula from the viewpoints of best approximations, generalized inverses and the Tikhonov regularization by combining these fundamental ideas and methods by means of the theory of reproducing kernels.
Furthermore, we shall use the sinc functions (the sinc method) as a
new approarch to solve the crucial Fredholm integral equation of thesecond kind on the half space in our general theory. We shall also
propose a new method for the real inversion of the Laplace transform
based essentially on linear simultaneous equations. We may think that
this real inversion formula is practical and natural. Error analysis willbe also considered.

2

Background General Theorems

Let $E$ be an arbitrary set, and let
be a reproducing kernel Hilbert
space (RKHS) admitting the reproducing kernel $K(p, q)$ on $E$ . For$H_{K}$

66
any Hilbert space
we consider a bounded linear operator from
into . We are generally interested in the best approximation$L$$\mathcal{H}$$H_{K}$$\mathcal{H}$problem(2.2)$\inf_{f\in H_{K}}\Vert Lf-d\Vert_{\mathcal{H}}$for a vector in . However, this extremal problem is quite involvedin existence and representation. See [16,19,20] for the details.
Now, for the Tikhonov regularization, we set, for any fixed positive$d$$\mathcal{H}$$\alpha>0$$K_{L}(\cdot,p;\alpha)=(L^{*}L+\alpha I)^{-1}K(\cdot,p)$

wheredenotes the adjoint operator ofinner product$L^{*}$$L$,. Then, by introducing the$(f,g)_{H_{K}(L;\alpha)}=\alpha(f,g)_{H_{K}}+(Lf, Lg)_{\mathcal{H}}$

,

(2.3)

we shall construct the Hilbert spacecomprising functions of. This space, of course, admits a reproducing kernel. FMrthermore,
we dirctly obtain$H_{K}(L;\alpha)$$H_{K}$

Proposition 2.1 (/18-20]) The extremal function
regularization$f_{d,\alpha}(\rho)$(2.4)$\inf_{f\in K}\{\alpha||f||_{H_{K}}^{2}+\Vert d-Lf\Vert_{?t}^{2}\}$

exists uniquely and it is represented in terms of the kemel$by$in the Tikhonov$K_{L}(p, q;\alpha)$

:

(2.5)

$f_{d,\alpha}(p)=(d, LK_{L}(\cdot,p;\alpha))_{\mathcal{H}}$

where the kemel $K_{L}(p, q;\alpha)$ is the reproducing kemel for the Hilbert
space $H_{K}(L;\alpha)$ and it is determined as the unique solutionof the equation:$\tilde{K}(p, q;\alpha)$$\tilde{K}(p, q;\alpha)+\frac{1}{\alpha}(L\tilde{K}_{q}, LK_{p})_{\mathcal{H}}=\frac{1}{\alpha}K(p, q)$

(2.6)

with
$\tilde{K}_{q}=\tilde{K}(\cdot, q;\alpha)\in H_{K}$and$K_{p}=K(\cdot,p)\in H_{K}$for$q\in E$

for $p\in E$ .

,

(2.7)

67
In (2.5), when contains errors or noise, we need its error estimate.
For this, we can use the general result:$d$Proposition 2.2 (/14]). In (2.5), we have the estimate$|f_{d,\alpha}(p)| \leq\frac{1}{\sqrt{\alpha}}\sqrt{K(p,p)}\Vert d\Vert_{\mathcal{H}}$.

For the convergence rate or the results for noisy data, see, ([9]).

A Natural Situation for Real Inversion Formulas
3

In order to apply the general theory in Section 2 to the real inversion formula of the Laplace transform, we shall recal the “naturalsituation” based on [17].
We shall introduce the simple reproducing kemel Hilbert space(RKHS)
comprised of absolutely continuous functions $F$ on thepositive real linewith finite norms$H_{K}$$R^{+}$$\{\int_{0}^{\infty}|F’(t)|^{2}\frac{1}{t}e^{t}dt\}^{1/2}$

and satisfying $F(O)=0$ . This Hilbert space admits the reproducing
kernel $K(t, t’)$
(3.8)

$K(t, t’)= \int_{0}^{\min(t,t’)}\xi e^{-\xi}d\xi$$=\{\begin{array}{llll}-te^{-t}- e^{-t}+1 for t\leq t’-t’e^{-t’} -e^{-t}’+1 for t\geq t\end{array}\}$(see [16], pages 55-56). Then we see that$\int_{0}^{\infty}|(\mathcal{L}F)(p)p|^{2}dp\leq\frac{1}{2}||F||_{H_{K}}^{2}$

;

(3.9)

that is, the linear operator$(\mathcal{L}F)(p)p$oninto $L_{2}(R^{+}, dp)=L_{2}(R^{+})$ is bounded. For the reproducing
kernel Hilbert spaces $H_{K}satis\phi ing(3.9)$ , we can find some general$H_{K}$

68
spaces ([17]). Therefore, from the general theory in Section 2, we
obtain
Proposition 3.1 $(l^{15}7)$ . For any $g\in L_{2}(R^{+})$ and for any $\alpha>0$ , thein the $8ense$
best approximation$F_{\alpha,g}^{*}$$\inf_{F\in H_{K}}\{\alpha\int_{0}^{\infty}|F’(t)|^{2}\frac{1}{t}e^{t}dt+\Vert(\mathcal{L}F)(p)p-g||_{L_{2}(R+}^{2})\}$(3.10)$= \alpha\int_{0}^{\infty}|F_{\alpha,g}^{*l}(t)|^{2}\frac{1}{t}e^{t}dt+||(\mathcal{L}F_{\alpha,g}^{*})(p)p-g\Vert_{L_{2}(R^{+})}^{2}$exists uniquely and we obtain the representation$F_{\alpha,g}^{*}(t)= \int_{0}^{\infty}g(\xi)(\mathcal{L}K_{\alpha}(\cdot,t))(\xi)\xi d\xi$

Here,$K_{\alpha}(\cdot, t)$.

(3.11)

is determined by the functional equation$K_{\alpha}(t, t’)= \frac{1}{\alpha}K(t,t’)-\frac{1}{\alpha}((\mathcal{L}K_{\alpha,t’})(p)p, (\mathcal{L}K_{t})(p)p)_{L_{2}(R+})$(3.12)

for$K_{\alpha,t’}=K_{a}(\cdot, t’)$and$K_{t}=K(\cdot,t)$.

4 Sampling Theory and Reproducing
Kernels
In order to solve the integral equation (3.11), numerically, we shall
employ the sinc method. At first we $shaU$ fix notations and basic
results in the sampling theory following the book by F. Stenger[23]and at the same time we shall show the basic relation of the sampling
theory and the theory of reproducing kernels.
We shall consider the integral transform, for a function in$g$$L_{2}(-\pi/h, +\pi/h),$ $(h>0)$69$f(z)= \frac{1}{2\pi}\int_{-\pi/h}^{\pi/h}g(t)e^{-izt}dt$.

(4.13)In order to identify the image space following the theory ofreproducing
kernels [16], we form the reproducing kernel$K_{h}(z, \overline{u})=\frac{1}{2\pi}\int_{-\pi/h}^{\pi/h}e^{-izt}\overline{e^{-iut}}dt$(4.14)$= \frac{1}{\pi(z-\overline{u})}\sin\frac{\pi}{h}(z-\overline{u})$$:= \frac{1}{h}$Sinc$( \frac{z-\overline{u}}{h})$$:= \frac{1}{h}S(k, h)(z-\overline{u}+hk)$,

by the notations in [23]. The image space of (4.13) is called the PaleyWiener spacecomprised of all analytic functions of exponential
type $satis\theta ing$ , for some constant $C$ and as $zarrow\infty$$W( \frac{\pi}{h})$$|f(z)|\leq C$

exp$( \frac{\pi|z|}{h})$and$\int_{R}|f(x)|^{2}dx<\infty$.

From the identity$K_{h}(jh,j’h)= \frac{1}{h}\delta(j,j’)$(the Kronecker’s), since

is the reproducing kemel for the
Hilbert space , from the general theory of integral transforms andthe Parseval’s identity we have the isometric identities in (4.13)$\delta$$\delta(j,j’)$$\ell^{2}$$\frac{1}{2\pi}\int_{-\pi/h}^{\pi/h}|g(t)|^{2}dt$$=h \sum_{j}|f(jh)|^{2}$$= \int_{R}|f(x)|^{2}dx$.

70( , Of) is
That is, the reproducing kernel Hilbert spacewithcharacterized as a space comprising the Paley-Wiener space
and with the norm squares above. Here we used the well-known result; that is,thatis a unique set for the Paley-Wiener space$f(jh)=0$ for all implies $f\equiv 0$ . Then, the reproducing property ofstates that
$K_{h}$$H_{K_{h}}$$z$$W( \frac{\pi}{h})$$\{jh)\}_{j}$$W( \frac{\pi}{h})$$j$$K_{h}(z,\overline{u})$$f(x)= \langle f(\cdot), K_{h}(\cdot, x)\rangle_{H_{K_{h}}}=h\sum_{j}f(jh)K_{h}(jh, x)$$= \int_{R}f(\xi)K_{h}(\xi,x)d\xi$.

, this representation is the sampling
theorem which represents the whole data $f(x)$ in terms of the discrete
data $\{f(jh)\}_{j}$ . For a general theory for the sampling theory and error, see [16].
estimates for some finite points
In particular, on the real line$x$$\{hj\}_{j}$5

New Algorithm

By setting$(\mathcal{L}K_{\alpha}(\cdot, t))(\xi)\xi=H_{\alpha}(\xi,t)$,

which is needed in (3.11), we obtain the Fredholm integral equationof the second type$\alpha H_{\alpha}(\xi, t)+\int_{0}^{\infty}H_{\alpha}(p,t)\frac{1}{(p+\xi+1)^{2}}dp$$=f( \xi, t)=-\frac{e^{-t\xi}e^{-t}}{\xi+1}(t+\frac{1}{\xi+1})+\frac{1}{(\xi+1)^{2}}$.

(5.15)We shall use the double exponential transform $f_{0}n_{oW}ing$ the idea [24]$\xi=\phi(x)=\exp$ ($\phi’(x)=\frac{\pi}{2}$$\frac{\pi}{2}$sinh ),

cosh $x\exp$ ($x$$\frac{\pi}{2}$sinh ).
$x$is a monotonically increasing function and $\phi(-\infty)=$Note that thisand $\phi(\infty)=\infty$ . In addition, for examples$\phi(x)$$0$$\phi(-4)=2.416\cross 10^{-19}$

71
and$\phi(16)=5.860\cross 10^{3030999}$.

So there is no need for setting so wide interval of integration from apractical point of view.
Then, we have
$H_{\alpha}(\xi,t)=H_{\alpha}(\phi(x),t)=\tilde{H}_{\alpha}(x, t)$and so,$\tilde{H}_{\alpha}(x, t)=\sum_{j_{=-\infty}}^{j=\infty}\tilde{H}_{\alpha}(jh, t)Sinc(\frac{x}{h}-j)$and$\alpha\tilde{H}_{\alpha}(x,t)$$+ \int_{-\infty}^{\infty}\tilde{H}_{\alpha}(z, t)\frac{1}{(\phi(z)+\phi(x)+1)^{2}}\phi^{j}(z)dz=f(\phi(x), t)$.

(5.16)We shall approximate as follows:$\tilde{H}_{\alpha}(x,t)\simeq\sum_{j=-N}^{j=N}\tilde{H}_{\alpha}(jh, t)Sinc(\frac{x}{h}-j)$

For error estimates for some finite points
have$\{hj\}_{j}$.

, see [16]. Then, we$\alpha\sum_{j}\tilde{H}_{\alpha}(jh,t)Sinc(\frac{x}{h}-j)$$+ \int_{-\infty}^{\infty}\sum_{k}\tilde{H}_{\alpha}(kh, t)Sinc(\frac{z}{h}-k)\frac{1}{(\phi(z)+\phi(x)+1)^{2}}\phi’(z)dz=f(\phi(x),t)$(5.17)

Rom the identities$\int_{-\infty}^{\infty}Sinc(\frac{x}{h}-i)Sinc(\frac{x}{h}-j)dx=h\delta_{ij}$and$\int_{-\infty}^{\infty}\frac{1}{(\phi(z)+\phi(x)+1)^{2}}Sinc(\frac{x}{h}-l)dx=\frac{1}{(\phi(z)+\phi(lh)+1)^{2}}h$by setting$A_{lk} \equiv\frac{1}{(\phi(kh)+\phi(lh)+1)^{2}}\phi’(kh)h$

,

.

72
we obtain the equation$\alpha\tilde{H}_{\alpha}(lh, t)+\sum_{k}\tilde{H}_{\alpha}(kh, t)A_{lk}=f(\phi(lh), t)\equiv f(lh, t)$(5.18)and the representation$F^{*}(t)= \int_{0}^{\infty}g(\xi)H_{\alpha}(\xi,t)d\xi=\int_{-\infty}^{\infty}g(\phi(x))H_{\alpha}(\phi(x),t)\phi’(x)dx$$\simeq h\sum_{i}g(\phi(ih))\tilde{H}_{\alpha}(ih,t)\phi’(ih)$6.

Inverses for More General FUnctions

As one of the main features of our method, we can easily generalize theapproximation function space. By a suitable transform, our inversion
formula in Section 5 is applicable for more general functions as follows:
We assume that $F$ satisfies the properties (P):$F\in C^{1}[0, \infty)$,$F’(t)=o(e^{\alpha t})$,$0< \alpha<k-\frac{1}{2}$$F(t)=o(e^{\beta t})$,$0<\beta<k$

,

and
一

$\frac{1}{2}$Then, the function(6.19)$G(t)=\{F(t)-F(O)-tF’(0)\}e^{-kt}$

belongs to$H_{K}$. Then,$( \mathcal{L}G)(p)=f(p+k)-\frac{F(0)}{p+k}-\frac{F’(0)}{(p+k)^{2}}$.

(6.20)

Therefore, if we know $F(O)$ and $F’(O)$ , then from$g(p)=(\mathcal{L}G)(p)$by the method in Section 5, we obtain $G(t)$ and so, $hom$ the identity$F(t)=G(t)e^{kt}+F(0)+tF’(0)$(6.21)

we have the inverse $F(t)hom$ the data $f(p),$ $F(O)$ and $F’(O)$ throughthe above procedures.

73

7

Numerical Experiments

We used $h=0.05$ and for , we take the span with 0.01. For thesimultaneous equations (5.18), we take from $\ell=-200$ to $\ell=799$ ;that is, 1000 equations. We solved such equations for $[0,4.99]$ withthe span 0.01 for .
We shall give a numerical experiment for the typical example
$t$

$t$

$F_{0}(t)=K(t, 1)=\{\begin{array}{ll}\text{一} te^{-t}-e^{-t}+1 for 0\leq t\leq 11-2e^{-1} for 1\leq t,\end{array}$whose Laplace transform is$( \mathcal{L}F_{0})(p)=\frac{1}{p(p+1)^{2}}[1-(p+2)e^{-(p+1)]}$.

(7.22)

Figure 1: For $F_{0}(t)=K(t, 1)$ and for $\alpha=10^{-1},10^{-4},10^{-8},10^{-12},10^{-16}$ .

74

Figure

2:

For

the

$10^{-1},10^{-4},10^{-8},10^{-12},10^{-16}$

.

step

function

Figure 3: For the characteristic function$10^{-1},10^{-4},10^{-8},10^{-12},10^{-16}$.

$F_{0}(t)$$F_{0}(t)$

and

for

$\alpha$on [1/2, 3/2] and for$\alpha=$75

Figure 4: For the function $F_{0}(t)=1/4-1/4e^{-2t}(1+2t)$ and for $\alpha=$$10^{-1},10^{-2},10^{-3},10^{-4}$.

Acknowledgements

Al-Shuaibi visiting Gunma University was supported by the Japan
Cooperation Center, Petroleum, the Japan Petroleum Institute (JPI)and King Fahd University of Petroleum and Minerals, and he wishes
to express his deep thanks Professor Saburou Saitoh and Mr. Hideki
Konishi of the JPI for their kind hospitality. S. Saitoh is upported inpart by the Grant-in-Aid for Scientific Research (C)(2)(No. 16540137;
No. 19540164) from the Japan Society for the Promotion Science. S.
Saitoh and T. Matsuura are partially supported by the Mitsubishi
Foundation, the 36th, Natural Sciences, No. 20 (2005-2006).
$s$References[1] Abdulaziz Al-Shuaibi, On the inversion of the Laplace transformby use of a regularized displacement opera $tor$, Inverse Problems,13(1997), 1153-1160.
[2] Abdulaziz Al-Shuaibi, The Riemann zeta function used in the inversion of the Laplace transform, Inverse Problems, 14(1998), 1-7.

76
[3] Abdulaziz Al-Shuaibi, Inversion of the Laplace transform via
Post- Widder formula, Integral Transforms and Special Functions,11(2001), 225-232.
[4] K. Amano, S. Saitoh and A. Syarif, A Real Inversion Fomula for
the Laplace $\pi_{ansform}$ in a Sobolev Space, Z. Anal. Anw. 18(1999),1031-1038.
[5] K. Amano, S. Saitoh and M. Yamamoto, Error estimates of the real
$inver8ion$ fomulas of the Laplace transform, Integral Transforms
and Special Functions, 10(2000), 165-178.

[6] A. Boumenir and A. Al-Shuaibi, The Inverse Laplace $\pi_{ansform}$$0perator8$ , J. Math. Anal.
and Analytic Pseudo-Differential
Appl. 228(1998), 16-36.
[7] D.-W. Byun and S. Saitoh, A real inversion fomula for the Laplacetransform, Z. Anal. Anw. 12(1993), 597-603.
[8] G. Doetsch, Handbuch der Laplace
Birkhaeuser Verlag, Basel, 1950.

$?$}$nnsformation$ , Vol. 1.,[9] H. W. Engl, M. Hanke and A. Neubauer, Regularization of Inverse
Problems, Mathematics and Its Applications, 376(2000), Kluwer

Academic Publishers.
[10] V. V. Kryzhniy, Regularized $Inver8ion$ of Integral $\pi ansforma-$, Inverse Problems, 19(2003),
tions of Mellin Convolution
1227-1240.
[11] V. V. Kryzhniy, Numeri cal inversion of the Laplace transform:$\Phi^{pe}$analysis via regularized analytic continuation, Inverse Problems,22(2006), 579-597.

[12] T. Matsuura and S. Saitoh, Analytical and numerical solutionsof linear ordinary differential equations with constant coefficients,
Journal of Analysis and Applications, 3(2005), pp. 1-17.

[13] T. Matsuura and S. Saitoh, Analytical and numerical real inversion formulas of the Laplace transform, The ISAAC Catani
Congress Proceedings (to appear).
[14] T. Matsuura and S. Saitoh, Analytical and numencal inversion
fomulas in the $Gau88ian$ convolution by using the Paley- Wiener
spaces, Applicable Analysis, 85(2006), 901-915.
[15] E. L. Post, Generalized diffentiation, ‘TMrans. Amer. Math. Soc.
32(1930), 723-781.

77[16] S. Saitoh, Integral Iansforms, Reproducing Kemels and their
Applications, Pitman Res. Notes in Math. Series 369, Addison
Wesley Longman Ltd (1997), UK.

[17] S. Saitoh, Approximate real inversion formulas of the Laplacetransform, Far East J. Math. Sci. 11(2003), 53-64.
[18] S. Saitoh, Approximate Real Inversion Formulas of the Gaussian
Convolution, Applicable Analysis, 83(2004), 727-733.
[19] S. Saitoh, Best approximation, Tikhonov rtegularization and reproducing $kernel_{8}$ , Kodai. Math. J. 28(2005), 359-367.
[20] S. Saitoh, Tikhonov regulanzation and the theory of reproducing
kemels, Proceedings of the 12th International Conference on Finite or Infinite Dimensional Complex Analysis and Applications,
Kyushyu University Press (2005), 291-298.
[21] S. Saitoh, Vu Kim Tuan and M. Yamamoto, Conditional Stabilityof a Real Inverse Formula for the Laplace $\pi an8form$, Z. Anal.
Anw. 20(2001), 193-202.
[22] S. Saitoh, N. Hayashi and M. Yamamoto (eds.), Analytic Extension Formulas and their Applications, (2001), Kluwer Academic

Publishers.
[23] F. Stenger, Numerical Methods Based on Sinc and Analytic hnctions, Springer Series in Computational Mathematics, 20, 1993.
[24] H. Ibkahasi and M. Mori, Double exponential formulas for numerical integration, Publ. Res. Inst. Math. Sci. 9(1974), 511-524.
[25] Vu Kim Tuan and Dinh Thanh Duc, Convergence rate of PostWidder appro rzmate inversion of the Laplace transform, Vietnam
J. Math. 28(2000), 93-96.
[26] Vu Kim Tuan and Dinh Thanh Duc, A new real inversion formulafor the Laplace transform and its convergence rate, Frac. Cal. &Appl. Anal. 5(2002), 387-394.
[27] D. V. Widder, The inversion of the Laplace integral and the related moment problem, hans. Amer. Math. Soc. 36(1934), 107-

2000.

[28] D. V. Widder, The Laplace $\pi_{ansf_{07}m}$ , Princeton University
Press, Princeton, 1972.
[29] $http://librury.wolffam.com/inforoenter/MathSource/4738/$[30] $http://www$ . columbia. $edu/ww2\theta 4\theta/abate.html$

78
T. Matsuura
Department of Mechanical System Engineering,
Graduate School of Engineering,
Gunma University, Kiryu, 376-8515, Japan

E-mail: matsuura@me.gunma-u.ac.jp

Abdulaziz Al-Shuaibi
KFUPM Box 449, Dhahran 31261, Saudi Arabia
E-mail: shuaaziz@kfupm.edu.sa
H. Fujiwara
Graduate School of Informatics,
Kyoto University, Japan
E-mail: fujiwara@acs. .kyoto-u.ac.jp$i$S. Saitoh
Department of Mathematics,
Graduate School of Engineering,
Gunma University, Kiryu, 376-8515, Japan
E-mail: ssaitoh@math.sci.gunma-u.ac.jp
and M. Sugihara
Department of Computer Science and Engineering,
Graduate School of Engineering, University of Tokyo
Bunkyou-Ku, Hongou, 113-8656, Japan

E-mail: m-sugihara@mist.i.u-tokyo.ac.jp