zzの逆関数（3） - inamori’s diary

前回、グラフで解の存在だけ示したが、
それだけではさみしいので、具体的数値を求めてみよう。
グラフによると、wが0に近いところに解があるので、
それを求めてみる。

$v + e^v\cos{w} = log(2\pi n)$
$w + e^v\sin{w} = \frac{\pi}{2}$

下より、wは小さいので、

$e^v \approx \frac{\frac{\pi}{2} - w}{w} \equiv X$

これを上に代入して、

$\log{X} + X \approx log(2\pi n)$

Xはそれなりに大きいので、log(2πn)として、
あとは、ニュートン法あたりで数値的に求めればよいだろう。



import std.cstream;

import std.math;
void main() {

    for(int n = 1; n <= 20; n++)

        dout.writefln("%d,%f", n, rev_zz(n));

}
cdouble rev_zz(int n) {

    double  X = 2 * std.math.PI * n;

    double  w = std.math.PI / 2 / (1 + log(X));

    double  v = log( (std.math.PI / 2 - w) / sin(w));

    cdouble u = exp(v + 1i * w);

    

    cdouble du;

    do {

        du = u * exp(u) - X * 1i;

        u -= du / (u + 1) / exp(u);

    } while(abs(du) > 1e-10);

    

    return exp(u);

}
cdouble exp(cdouble c) {

    double  im = c.im;

    return std.math.exp(c.re) * (cos(im) + 1i * sin(im));

}

nが1から20まで求めてみた。



1,2.213583+3.114000i

2,3.034197+5.224328i

3,3.717612+7.108549i

4,4.326578+8.867915i

5,4.886568+10.543436i

6,5.411000+12.156942i

7,5.908002+13.721814i

8,6.382954+15.247012i

9,6.839641+16.738916i

10,7.280851+18.202272i

11,7.708711+19.640738i

12,8.124888+21.057206i

13,8.530719+22.454017i

14,8.927297+23.833098i

15,9.315529+25.196062i

16,9.696178+26.544277i

17,10.069895+27.878912i

18,10.437239+29.200982i

19,10.798695+30.511372i

20,11.154687+31.810862i