zzの逆関数(3)

前回、グラフで解の存在だけ示したが、
それだけではさみしいので、具体的数値を求めてみよう。
グラフによると、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