月別アーカイブ: 2016年12月

数学

ニュートン写像に現れる綺麗な画像

コメントする

ニュートン写像はニュートン=ラフソン法によって、ある初期値から適当な値へ収束するまでの回数でその関数を特徴づけるものです。
フラクタルと似ている考えであることを明記しておきます。

ニュートン写像

\(
\displaystyle N_f(z)=z-\frac{f(z)}{f'(z)}
\)
の繰り返しによって適当な初期値からスタートし、収束させていきます。このページでは画像が主です。説明の方はリーマンのゼータ関数の複素力学系をご覧ください。

基本はこれだけです。簡単な複素関数でも相当複雑になります。
早速適当な関数を用意して、ニュートン写像を見ていきましょう。

[adsense1]

まずは基本的な多項式を見ます。
\(
f(z)=z^2-1
\)

\(
f(z)=z^3-1
\)

\(
f(z)=z^4-1
\)

\(
\displaystyle f(z)=z^5-1, z^6-1, z^7-1, z^10-1
\)

少し複雑な関数の場合
\(
\displaystyle f(z)=\ln z, e^z-1
\)

\(
\displaystyle f(z)=e^{1/z}-0.01
\)

\(
\displaystyle f(z)=ze^{z}\ln z-0.01
\)

\(
\displaystyle f(z)=z^5e^z-0.1
\)

\(
\displaystyle f(z)=z^4e^z\ln z-0.1
\)

\(
\displaystyle f(z)=sin(z)/z
\)

\(
\displaystyle f(z)=z^2 e^{-z^2}-1
\)

\(
\displaystyle f(z)=ln(z) e^{-z^2}-1
\)

最後にリーマンのゼータ関数

複素関数\(f(z)\)が多くの根を持つとき、非常に複雑で綺麗な画像が得られるようです。
恐らく遠く離れたところからたまたま根のあたりに到達することが多くなるからでしょう。

[adsense2]

fortran90によるプログラム

fortran90によるプログラムを置いておきます。入力するべき場所は!=====の間と、関数fと導関数dfです。

program main
  implicit none
  integer,parameter::kmax=49
  double precision,parameter::eps=1d-8
  integer::i,j,k,Nx,Ny
  double precision::xa,xb,ya,yb,x,y,hx,hy
  complex(kind(0d0))::f,df,z,z1
  external::f,df

  !==================
  xa=-5d0; xb=5d0
  Nx=201

  ya=-5d0; yb=5d0
  Ny=201
  !==================
   
  hx=(xb-xa)/dble(Nx-1)
  hy=(yb-ya)/dble(Ny-1)
  do i=0,Nx-1
     x=xa+i*hx
     do j=0,Ny-1
        y=ya+j*hy
        z=dcmplx(x,y)
        do k=0,kmax
           z1=z-f(z)/df(z)
           if(abs((z-z1)/z).le.eps)exit
           z=z1
        enddo
        write(10,'(2e18.8e3,i5)')x,y,k
     enddo
     write(10,*)
     if(mod(i,100).eq.0)write(6,'(A,i6,A,i6)')"  ",i," / ",Nx
  enddo
 
  stop
end program main

function f(z)
  implicit none
  complex(kind(0d0))::f,z
 
  f=z**5*exp(z)-0.1d0
 
  return
end function f

function df(z)
  implicit none
  complex(kind(0d0))::df,z
 
  df=5*z**4*exp(z)+z**5*exp(z)
 
  return
end function df
プログラミングと数値計算, 数学

リーマンのゼータ関数の複素力学系

コメントする

本稿の目的はリーマンのゼータ関数を、ニュートン法による写像を用いて別の視点から見ることです。
綺麗な絵が得られる[1]ので自分のコンピュータ上で見てみたいと思いました。

複素力学系とは？

複素力学系とは、簡単に言えば写像のことです。
実関数による写像の繰り返しで得られる系を力学系と呼び、
複素関数による写像の繰り返しで得られる系を複素力学系と呼びます。
”力学”の言葉の由来はニュートン力学から来ています。

”ニュートン力学”は目に見えるくらいの物体を扱う学問ですが、その正体は時間に関する2階の微分方程式です。
なので、時刻\(t\)の物体の振る舞い\(x(t)\)からちょっとの時間\(\Delta t\)秒後の物体の振る舞い\(x(t+\Delta t)\)は、適当な時間発展をさせる演算子\(f\)によって
\(
x(t)\overset{f}{\longrightarrow} x(t+\Delta t)
\)
と移り変わります。これを何度も繰り返し作用させて任意の時刻\(t’\)まで時間発展させる事が出来ます。

\(
x(t)\overset{f}{\longrightarrow} x(t+\Delta t)\overset{f}{\longrightarrow} x(t+2\Delta t)\overset{f}{\longrightarrow}x(t+3\Delta t)\overset{f}{\longrightarrow} \cdots\overset{f}{\longrightarrow} x(t’)
\)
となると、ニュートン力学を性質を知ろうとするときに重要な情報は、演算子\(f\)と初期値に集約されます。

演算子\(f\)の空間へどんどん写像していくわけですので、ある写像の繰り返しも力学のようだ、だから”力学系”という名前が付けられたわけです。
もう少し詳しく知りたければ、まずは[3]を読むと良いでしょう。

本稿では複素力学系の一つであるニュートン＝ラフソン法によって得られる写像(ここでは便宜的にニュートン写像と呼ぶことにします)を考えます。

複素関数\(f(z)\)の具体的なニュートン写像の数式は以下で与えられます。
\(
\displaystyle N_f(z)=z-\frac{f(z)}{f'(z)}
\)
ある複素平面の初期値からスタートして、この写像を何回繰り返したら収束したか？をその系を特徴づける値とします。
要は、ニュートン＝ラフソン法の初期値を複素平面上の全ての点を取り、収束するまでの回数をカウントするのです。

[adsense1]

ニュートン写像の具体例

本稿の目的はゼータ関数に対して行うことですが、まずは簡単な複素関数のニュートン写像を見てみましょう。
まず、複素関数
\(
\displaystyle f(z)=z^2-1
\)
を考えます。根の周り（\(z=\pm 1\)）では収束回数は少なくなることが予想できます。実際にやってみますとこんな感じに。

非常に単調です。実部の値が正か負かでのみ特徴づけられるようです。

続いて、複素関数
\(
\displaystyle f(z)=z^3-1
\)
を考えます。根は\(z=1, e^{\pm i2\pi/3}\)。

3つは確かに分けられていますが、その境界辺りでは非常に奇怪な模様になります。

さらに、複素関数
\(
\displaystyle f(z)=z^4-1
\)
を考えます。根は\(z=\pm 1, \pm i\)。

ゼータ関数のニュートン写像

では、ゼータ関数に適応してみましょう。
ゼータ関数の導関数は導出方法とプログラムは[2],[9]/リーマンのゼータ関数の数値計算コード(複素平面)に追記したのでそちらをご覧ください。
するとこんなニュートン写像になります。

ゼロ点付近では収束が早くなっています。ゼロ点のある場所がマークしたところです。

カラーを対数にしますとこんな図に。

(2017/01/13 追記)
もっとほかの範囲のニュートン写像を見ていきます。

全体像

拡大その1(全体像の左側の四角)

拡大その2(全体像の右側の四角)

拡大その２, 高解像度版 highres_zetanewton.png (15MB)

元の画像はこちらからダウンロードできます。

[adsense2]

ゼータ関数のニュートン写像を得るfortran90のプログラムはこちらです。

ifort -O3 -qopenmp main.f90

計算範囲はxa,xb,ya,ybで指定されています。

参考文献

[1]Tomoki Kawahira “Riemann’s zeta function and Newton’s method: Numerical experiments from a complex-dynamical viewpoint”

[2]Tom M. Apóstol “Formulas for Higher Derivatives of the Riemann Zeta Function”, Math. of. Comp vol 44, 169, 1985, pp. 223-232.

[3]Chapter 14 力学系（補講１）

[4]ジュリア集合の色付けを工夫して芸術的なフラクタル図形を描く

[5]複素力学系, ニュートン法, ジュリア集合

[6]西沢清子, 藤村雅代 “多項式のニュートン写像における鉢の幅について(カオスをめぐる力学系の諸問題)”

[7]Riemann Zeta Function -wolfram mathworld

[8]ディガンマ関数の数値計算 -シキノート
[9]リーマンのゼータ関数の数値計算コード(複素平面) -シキノート

プログラミングと数値計算, 数学

二重指数関数型数値積分

1件のコメント

二重指数関数型数値積分公式(DE公式またはtanh-sinh公式)とは、変数変換によって被積分関数を別の関数に変換し積分する方法のことです。
計算アイデアは

  1. 変数変換によって台形則が良く働く被積分関数に変換する
  2. 台形則によって数値積分する

という考えに立脚しています。端点に特異点が存在する場合に最適、という特徴を持ちます。
1974年に、高橋 秀俊と森 正武によって提案されました[1]。

二重指数関数型数値積分とは、関数\(f(x)\)の\([-1,1]\)に渡る積分を

\(
\begin{align}
\int_{-1}^{1} f(x) dx &\approx \sum_{i=-N_-}^{N^+} w_i f(x_i),\\
x_i&=\tanh\left(\frac{\pi}{2}\sinh(ih)\right), \\
w_i&=\frac{h\frac{\pi}{2}\cosh(ih)}{\cosh^2\left(\frac{\pi}{2}\sinh(ih)\right)}
\end{align}
\)

として近似します。\(N_-, N_+\)は離散化誤差と打ち切り誤差が等しくなるように決められます。
\(N_-, N_+\)はプログラム上では、\(x_i\)がコンピュータの扱える桁数を超えず、\(w_i\)がアンダーフローを起こさない範囲で決められます。


[adsense1]

計算アイデア

DE公式は変数変換型の数値積分公式と呼ばれます。
DE公式の計算は、

  1. 変数変換によって台形則が良く働く被積分関数に変換する
  2. 台形則によって数値積分する

という考えが元になっています。

大きな特徴として、

・端点特異性に強い
・応用範囲が広い

ことが挙げられます。
例えば、端点で発散してしまう積分
\(
\displaystyle \int_0^1 \frac{1}{\sqrt{x}} dx
\)
や、厄介な積分
\(
\displaystyle \int_0^1 \sin\left(\frac{1}{\sqrt{x}}\right)/\sqrt{x} dx
\)
の計算も少ない分点数で実行できます(後者の厄介な積分はDE公式を用いて4桁程度一致します。しかし、ほかの補間型の積分公式では1桁合うかどうかでしょう)。

なぜ変数変換でうまく計算できるのか、これを理解するにはまず、台形則が良く働く場合とは何かを知らなければなりません。

台形則が良く働く関数とは、関数の端点での微分の値が近い事です。
台形則の刻み幅を\(h\), 区間\([a,b]\)を台形則によって関数\(f(x)\)を数値積分する時、本来の積分値と台形則で求めた値の誤差を表す第1項は、
\(
\displaystyle \frac{h^2}{12}\{f'(b)-f'(a)\}
\)
と表されます。刻み幅を小さくすれば誤差が小さくなるのは直観的に分かりますが、誤差を小さくできる要因は\(\{f'(b)-f'(a)\}\)にもあります。これは、
・両端で微分の値が等しい(周期関数の1周期に渡る積分、\(a,b\)に近づくつれて関数が定数に近づいていく積分)
を考えるとき、高精度の結果を与えると言っています。

上記条件が満たされるとき、他の高次の誤差項\(f^{(n)}(b)-f^{(n)}(a)\)もほとんどゼロになることが期待されるため、数値計算誤差が限りなく小さくなっていきます。

では変数変換によって端点で減衰する被積分関数に変換することを考えましょう。
高橋-森によって以下の変数変換が提案され、その特性が調べられました[1]。
\(
\displaystyle x=\varphi(t)=\tanh\left(\frac{\pi}{2}\sinh(t)\right)
\)
この変数変換をすることによって端点で”急速に“ゼロに近づく被積分関数に変換されるのです。

この”急速に“の速度はすさまじく、tが大きくなる時、被積分関数は
\(
\displaystyle f(\varphi(t))\varphi'(t)\sim A\exp\{-c\exp(-|t|)\}
\)
の形を持ちます。指数の指数の形で減衰するため、二重指数関数型という名前が付けられているのです。

まぁ、なぜ簡単に計算できるのかと言えば、\(1/\sqrt{x}|_{x\to 0}\)で発散していく時、これよりも早く減衰する関数で割れば収束しますよね、という事です。

※三重、四重指数関数型があるのでは、と考えるのは自然です。しかし、このような関数は(多分ですが)正則であるという条件の下では存在しないことが示されています[3,4]。一重指数関数型数値積分公式は存在します(例えばIMT公式)。

ニュートン・コーツ型数値積分、ガウス求積法は補間型の数値積分公式と呼ばれます。
これらの公式は
与えられたデータ点とデータ点を(重み関数)×(多項式)として補間し、その補間された関数を積分する
というアイデアの元考えられています。
しかし、重み関数は方程式から人間が知っていないとならず、上の形以外の関数を無理やり積分しようとしても積分精度は悪いです。
そのため、問題ごとに合わせたアルゴリズムの変更が必要不可欠になります。

数値計算で気を付けるべき点

端点特異点がある場合の収束

端点特異点がある場合、変数変換後の空間で計算区間が足らなくなり、大きな打切り誤差が発生します。
端点で発散している時に顕著です。
例えば、積分
\(
\displaystyle \int_0^1 \sqrt{x}dx
\)
は高精度(16桁近く一致)の結果を与えますが、
積分
\(
\displaystyle \int_0^1 \frac{1}{\sqrt{x}}dx
\)
は8桁程度の一致しかしません。この原因は打切り誤差です。桁落ちが起きないように処理が必要な問題となります。

積分区間の変更

積分
\(
\displaystyle \int_a^b f(x) dx
\)
を行いたい場合、変数変換
\(
\displaystyle x=\frac{b-a}{2}t+\frac{b+a}{2}
\)
を行うことで区間\([-1,1]\)の積分に置き換えることが出来ます。

広義積分に対する二重指数関数型数値積分

半無限区間\([0,\infty]\)の場合、

\(
\begin{align}
\int_{0}^{\infty} f(x) dx &\approx \sum_{i=-N_-}^{N^+} w_i f(x_i),\\
x_i&=\exp\left(\frac{\pi}{2}\sinh(t)\right), \\
w_i&=h\frac{\pi}{2}\cosh(ih)\exp\left(\frac{\pi}{2}\sinh(ih)\right)
\end{align}
\)

無限区間\([-\infty,\infty]\)の場合、

\(
\begin{align}
\int_{-\infty}^{\infty} f(x) dx &\approx \sum_{i=-N_-}^{N^+} w_i f(x_i),\\
x_i&=\sinh\left(\frac{\pi}{2}\sinh(t)\right), \\
w_i&=h\frac{\pi}{2}\cosh(ih)\cosh\left(\frac{\pi}{2}\sinh(ih)\right)
\end{align}
\)

のように変数変換を行うことで二重指数で減衰する関数に変わります。

また、半無限区間の積分の多くの場合は\(exp(-x)\)で減衰します。
これを考慮すると
\(
x=\exp(t-\exp(-t))
\)
という変数変換が有効だということが分かります[1]。本稿では数値計算のアルゴリズムの都合上載せません。

計算の工夫

被積分関数の評価回数を可能な限り減らすため、アルゴリズムを工夫します。

変数変換後の空間で刻み幅\(h\)で求めた数値積分値\(I_h\)は
\(
\displaystyle I_h=h\sum_{i=-N}^N f(\varphi(ih))\varphi'(ih)
\)
と表されます。刻み幅を半分\(h/2\)にすると、分点数は\(4N+1\)になり、
\(
\displaystyle I_{h/2}={h/2}\sum_{i=-2N}^{2N} f(\varphi(ih/2))\varphi'(ih/2)
\)
と表されます。
\(I_h\)を利用して、\(I_{h/2}\)を表現すると、
\(
\displaystyle I_{h/2}=\frac{1}{2}\left\{I_h + h\sum_{i=-N}^{N-1} f(\varphi(ih+h/2))\varphi'(ih+h/2)\right\}
\)
と表されます。

この表式の利点は以下の通りです。
本来、\(I_{h/2}\)を計算するためには\(4N+1\)回右辺を計算しなければなりません。
しかし、\(I_{h}\)と\(I_{h/2}\)の分点の位置が重なっているときがちょうど2N+1点有ります。よって重なっていない残りの\(2N\)点だけを計算してやれば\(I_{h/2}\)を計算できるのです。評価回数が\(4N+1\)回から\(2N\)回に減ったことで単純に計算時間が半分になります。

[adsense2]

fortranプログラム

関数fの数値積分を実行します。
対応しているのは、

・1,2,3次元
・実関数/複素関数(実/複素引数)
・有限区間/半無限区間/無限区間

と

・極座標、球面座標での全空間積分

となっています。
数値計算精度は有効桁数の問題から、8桁程度と思うのが良いと思います。

複素関数の積分は、複素数点と複素数点を直線で結んだ線積分です。
また、半無限区間、無限区間の複素関数の積分では始点\(ca,wa,va\)と角度\(thx,thy,thz\)によって経路を指定することが出来ます。

厳密に全てのルーチンが正しく動作するかは確認していません。
利用する際は連絡先とサイトについてをご覧下さい。

モジュール↓(2000行近くあります)

具体的なプログラム例

積分
\(
\displaystyle \int_0^{5}\sin(\sqrt(x))
\)
を計算。

厳密値
4.3340264879445362
数値計算結果
4.33402648794427

プログラム

全具体例

実装してある計算を全ルーチンを利用したプログラムを書きます。
詳細はコメントを確認してください。

基本的にaを含む変数は積分の始点、bは積分の終点です。

program main
  use DE_integration
  implicit none
  integer::info
  double precision::xa,xb,ya,yb,za,zb,eps,s,thx,thy,thz
  complex(kind(0d0))::ca,cb,wa,wb,va,vb
  double precision,parameter::pi=dacos(-1d0)
  complex(kind(0d0))::cs
 
  double precision::f1,f2,f3,f4,p2,p3,fh1,fi1,fi2,fp2,fi3,fp3
  complex(kind(0d0))::g1,g2,g3,g4,g5,g6,h1,h2,h3,gh1,hh1,&
       gi1,hi1,gi2,hi2,gp2,gi3,hi3,gp3
  external::f1,f2,f3,f4,g1,g2,g3,g4,g5,g6,p2,p3,h1,h2,h3,&
       fh1,gh1,hh1,fi1,gi1,hi1,fi2,gi2,hi2,fp2,gp2,gi3,hi3,gp3,fi3,fp3
 
 
  ! +------- 1D Integration -------+

  ! real f, real x, 1D, finite-range
  xa=0d0;  xb=1d0;  
  eps=1d-12; s=0d0  
  call dde1d(f1,xa,xb,eps,s,info)
  write(6,*)info,s

  ! complex f, real x, 1D, finite-range
  xa=0d0;  xb=1d0;  
  eps=1d-12; s=0d0  
  call cde1d(g1,xa,xb,eps,cs,info)
  write(6,*)info,cs

  ! complex f, complex x, 1D, finite-range
  ca=dcmplx(0d0,1d0);  cb=dcmplx(1d0,3d0)  
  eps=1d-12; s=0d0  
  call zde1d(h1,ca,cb,eps,cs,info)
  write(6,*)info,cs

  ! real f, real x, 1D, semi-infinite
  xa=2d0;  
  eps=1d-12; s=0d0  
  call dde1d_hinf(fh1,xa,eps,s,info)
  write(6,*)info,s

  ! complex f, real x, 1D, semi-infinite
  xa=2d0;
  eps=1d-12; s=0d0
  call cde1d_hinf(gh1,xa,eps,cs,info)
  write(6,*)info,cs

  ! complex f, complex x, 1D, semi-infinite
  ca=dcmplx(0d0,0d0);  thx=pi/4d0
  eps=1d-12; s=0d0
  call zde1d_hinf(hh1,ca,thx,eps,cs,info)
  write(6,*)info,cs

  ! real f, real x, 1D, infinite
  xa=2d0;  
  eps=1d-12; s=0d0  
  call dde1d_inf(fi1,eps,s,info)
  write(6,*)info,s

  ! complex f, real x, 1D, infinite
  xa=2d0;
  eps=1d-12; s=0d0
  call cde1d_inf(gi1,eps,cs,info)
  write(6,*)info,cs

  ! complex f, complex x, 1D, infinite
  ca=dcmplx(0d0,0d0);  thx=pi/4d0
  eps=1d-12; s=0d0
  call zde1d_inf(hi1,ca,thx,eps,cs,info)
  write(6,*)info,cs
 
 
  ! +------- 2D Integration -------+
 
  ! real f, real x,y, 2D, finite-range
  xa=0d0; xb=1d0; ya=0d0; yb=2d0
  eps=1d-12; s=0d0  
  call dde2d(f2,xa,xb,ya,yb,eps,s,info)
  write(6,*)info,s

  ! complex f, real x,y, 2D, finite-range
  xa=0d0; xb=1d0; ya=0d0; yb=2d0
  eps=1d-12; s=0d0  
  call cde2d(g2,xa,xb,ya,yb,eps,cs,info)
  write(6,*)info,cs

  ! complex f, complex x,y, 2D, finite-range
  ca=dcmplx(0d0,1d0);  cb=dcmplx(1d0,3d0)  
  wa=dcmplx(1d0,-1d0);  wb=dcmplx(2d0,0d0)  
  eps=1d-12; s=0d0  
  call zde2d(h2,ca,cb,wa,wb,eps,cs,info)
  write(6,*)info,cs

  ! real f, real x,y, 2D, infinite-range
  eps=1d-8; s=0d0
  call dde2d_inf(fi2,eps,s,info)
  write(6,*)info,s

  ! complex f, real x,y, 2D, infinite-range
  eps=1d-12; s=0d0  
  call cde2d_inf(gi2,eps,cs,info)
  write(6,*)info,cs

  ! complex f, complex x,y, 2D, infinite-range
  ca=dcmplx(0d0,0d0);  thx=0d0
  wa=dcmplx(0d0,0d0);  thy=0d0
  eps=1d-12; s=0d0
  call zde2d_inf(hi2,ca,thy,wa,thx,eps,cs,info)
  write(6,*)info,cs

  ! real f, real x,y, 2D polar, infinite-range
  eps=1d-12; s=0d0
  call dde_polar(fp2,eps,s,info)
  write(6,*)info,s

  ! complex f, real x,y, 2D polar, infinite-range
  eps=1d-12; s=0d0  
  call cde_polar(gp2,eps,cs,info)
  write(6,*)info,cs
 
 
  ! +------- 3D Integration -------+

  ! real f, real x,y,z, 3D, finite-range
  xa=0d0; xb=1d0; ya=0d0; yb=2d0; za=2d0; zb=3d0
  eps=1d-12; s=0d0
  call dde3d(f3,xa,xb,ya,yb,za,zb,eps,s,info)
  write(6,*)info,s

  ! complex f, real x,y,z, 3D, finite-range
  xa=0d0; xb=1d0; ya=0d0; yb=2d0; za=2d0; zb=3d0
  eps=1d-12; s=0d0  
  call cde3d(g3,xa,xb,ya,yb,za,zb,eps,cs,info)
  write(6,*)info,cs

  ! complex f, complex x,y,z, 3D, finite-range  
  ca=dcmplx(0d0,1d0);  cb=dcmplx(1d0,3d0)  
  wa=dcmplx(1d0,-1d0);  wb=dcmplx(2d0,0d0)  
  va=dcmplx(1d0,1d0);  vb=dcmplx(0d0,1d0)
  eps=1d-12; s=0d0  
  call zde3d(h3,ca,cb,wa,wb,va,vb,eps,cs,info)
  write(6,*)info,cs

  ! real f, real x,y,z, 3D, finite-range
  eps=1d-12; s=0d0
  call dde3d_inf(fi3,eps,s,info)
  write(6,*)info,s

  ! complex f, real x,y,z, 3D, finite-range
  eps=1d-12; s=0d0  
  call cde3d_inf(gi3,eps,cs,info)
  write(6,*)info,cs

  ! complex f, complex x,y,z, 3D, finite-range
  ca=dcmplx(0d0,1d0);  thx=0d0
  wa=dcmplx(1d0,-1d0); thy=pi/2d0
  va=dcmplx(1d0,1d0);  thz=0d0
  eps=1d-12; s=0d0  
  call zde3d_inf(hi3,ca,thx,wa,thy,va,thz,eps,cs,info)
  write(6,*)info,cs

  ! real f, real x,y, 3D spherical, infinite-range
  eps=1d-12; s=0d0
  call dde_polar(fp3,eps,s,info)
  write(6,*)info,s

  ! complex f, real x,y, 3D spherical, infinite-range
  eps=1d-12; s=0d0  
  call cde_polar(gp3,eps,cs,info)
  write(6,*)info,cs
 
  stop
end program main

function f1(x)
  implicit none
  double precision,intent(in)::x
  double precision::f1
 
  f1=sin(sqrt(x))
 
  return
end function f1

function g1(x)
  implicit none
  double precision,intent(in)::x
  complex(kind(0d0))::g1
 
  g1=dcmplx(sin(sqrt(x)),exp(-x))
 
  return
end function g1

function h1(z)
  implicit none
  complex(kind(0d0)),intent(in)::z
  complex(kind(0d0))::h1
 
  h1=sin(z)
 
  return
end function h1

function fh1(x)
  implicit none
  double precision,intent(in)::x
  double precision::fh1
 
  fh1=exp(-x)
 
  return
end function fh1

function gh1(x)
  implicit none
  double precision,intent(in)::x
  complex(kind(0d0))::gh1
 
  gh1=dcmplx(sqrt(x)*exp(-x),exp(-x))
 
  return
end function gh1

function hh1(z)
  implicit none
  complex(kind(0d0)),intent(in)::z
  complex(kind(0d0))::hh1
  double precision,parameter::pi=dacos(-1d0)
 
  hh1=exp(dcmplx(0d0,1d0)*0.5d0*pi*z**2)
 
  return
end function hh1

function fi1(x)
  implicit none
  double precision,intent(in)::x
  double precision::fi1
 
  fi1=exp(-x*x)
 
  return
end function fi1

function gi1(x)
  implicit none
  double precision,intent(in)::x
  complex(kind(0d0))::gi1
 
  gi1=dcmplx(1d0/(1d0+x*x),exp(-x*x))
 
  return
end function gi1

function hi1(z)
  implicit none
  complex(kind(0d0)),intent(in)::z
  complex(kind(0d0))::hi1
  double precision,parameter::pi=dacos(-1d0)
 
  hi1=1d0/(1d0+z*z)
 
  return
end function hi1


function f2(x,y)
  implicit none
  double precision,intent(in)::x,y
  double precision::f2
 
  f2=sin(sqrt(x))*exp(-y)
 
  return
end function f2

function g2(x,y)
  implicit none
  double precision,intent(in)::x,y
  complex(kind(0d0))::g2
 
  g2=dcmplx(sin(sqrt(x))/sqrt(y),exp(-x)*y)
 
  return
end function g2

function h2(z,w)
  implicit none
  complex(kind(0d0)),intent(in)::z,w
  complex(kind(0d0))::h2
 
  h2=sin(z)*w
 
  return
end function h2

function fi2(x,y)
  implicit none
  double precision,intent(in)::x,y
  double precision::fi2
 
  fi2=1d0/(1d0+x**4+y**4)
 
  return
end function fi2

function gi2(x,y)
  implicit none
  double precision,intent(in)::x,y
  complex(kind(0d0))::gi2
 
  gi2=dcmplx(exp(-y*y-x*x),exp(-y*y-x*x))
 
  return
end function gi2

function hi2(z,w)
  implicit none
  complex(kind(0d0)),intent(in)::z,w
  complex(kind(0d0))::hi2
 
  hi2=exp(-z*z-w*w)
 
  return
end function hi2

function fp2(r,theta)
  implicit none
  double precision,intent(in)::r,theta
  double precision::fp2
 
  fp2=r*exp(-r)*sin(theta)**2
  fp2=r*fp2 ! Jacobian
 
  return
end function fp2

function gp2(r,theta)
  implicit none
  double precision,intent(in)::r,theta
  complex(kind(0d0))::gp2
 
  gp2=dcmplx(0d0,1d0)*r*exp(-r)*sin(theta)**2
  gp2=r*gp2 ! Jacobian
 
  return
end function gp2

function f3(x,y,z)
  implicit none
  double precision,intent(in)::x,y,z
  double precision::f3
 
  f3=sin(x*z)*exp(-y)
 
  return
end function f3

function g3(x,y,z)
  implicit none
  double precision,intent(in)::x,y,z
  complex(kind(0d0))::g3
 
  g3=dcmplx(sin(sqrt(x))/sqrt(y),exp(-x)*y*sqrt(z))
 
  return
end function g3

function h3(z,w,v)
  implicit none
  complex(kind(0d0)),intent(in)::z,w,v
  complex(kind(0d0))::h3
 
  h3=sin(z)*w*v**2
 
  return
end function h3

function fi3(x,y,z)
  implicit none
  double precision,intent(in)::x,y,z
  double precision::fi3
 
  fi3=exp(-x*x-y*y-z*z)
 
  return
end function fi3

function gi3(x,y,z)
  implicit none
  double precision,intent(in)::x,y,z
  complex(kind(0d0))::gi3
 
  gi3=dcmplx(exp(-y*y-x*x-z*z),2d0*exp(-y*y-x*x-z*z))
 
  return
end function gi3

function hi3(z,w,v)
  implicit none
  complex(kind(0d0)),intent(in)::z,w,v
  complex(kind(0d0))::hi3
 
  hi3=exp(-abs(z)**2-abs(w)**2-abs(v)**2)
 
  return
end function hi3

function fp3(r,theta,phi)
  implicit none
  double precision,intent(in)::r,theta,phi
  double precision::fp3
 
  fp3=r*exp(-r)*sin(theta)**2
  fp3=r*r*sin(theta)*fp3 ! Jacobian
 
  return
end function fp3

function gp3(r,theta,phi)
  implicit none
  double precision,intent(in)::r,theta,phi
  complex(kind(0d0))::gp3
 
  gp3=dcmplx(0d0,1d0)*r*exp(-r)*sin(theta)**2
  gp3=r*r*sin(theta)*gp3 ! Jacobian
 
  return
end function gp3

参考文献

[1] H. Takahasi and M. Masatake, “Double Exponential Formulas for Numerical Integration” (1974)
[2]渡辺二太, 二重指数関数型数値積分公式について(1990)

[3] 森 正武、「数値解析における二重指数関数型変換の最適性」
[4] M. Sugihara, Optimality of the double exponential formula-functional analysis approach,
Numer. Math. 75 (1997) 379-395.

数値積分法の入門とDE公式の説明として、
[5] 戸田 英雄, 小野 令美 数値解析における一つの話題

[6] 森 正武「二重指数関数型変換のすすめ」

[7]森 正武, “FORTRAN 77 数値計算プログラミング”、p.174-176岩波書店、(1986)第一刷