4.6 Eigen ライブラリィを用いた計算プログラム

個人的には、この手の計算には MATLAB を使うのが好みであるが、 ここでは、C++ でベクトル、 行列に関する演算をするための計算ライブラリィ Eigen を利用してみる。

WWWサイトから、 eigen-eigen-bdd17ee3b1b3.tar.gz のようなソース・ファイル一式を 入手して、以下のようにインストールする (/usr/include/ にファイルをインストールする場合)。

MacBook のターミナルで次のようにインストール

tar xzf eigen-eigen-bdd17ee3b1b3.tar.gz
cd eigen-eigen-bdd17ee3b1b3
tar cf - Eigen | (cd /usr/include/; sudo tar xpf -)

$\Omega=D_1=D(0;1)$ , $z_0=1/2$ の場合を解いてみよう。 この場合は、

$$\varphi(z)=\eps\frac{z-z_0}{1-\overline{z_0} z}$$

が解であることが分かっている ($\eps$ は $\vert\eps\vert=1$ を満たす定数)。

conformalmap.cpp

/*
 * conformalmap.cpp --- 等角写像の数値計算例 (基本解の方法の天野要バージョン)
 *
 * g++ -I/usr/X11R6/include conformalmap.cpp -L~/lib -lglscd -L/usr/X11R6/lib -lX11
 */

#include <iostream>
#include <complex>
#include <Eigen/Dense>

extern "C" {
#include <glsc.h>
};

using namespace Eigen;

// MatrixXd は大きさが指定できて、成分が double
// d の変わりに i, f, cf, cd

// φ(z0)=0, φ'(z0)>0 を満たす等角写像
std::complex<double> phi(std::complex<double> z,
                         std::complex<double> z0)
{
  std::complex<double> w;
  w = (z-z0)/(1.0-std::conj(z0)*z);
  return w;
}

// 天野の方法による等角写像
std::complex<double> f(std::complex<double> z,
                       std::complex<double> z0,
                       int N,
                       double Q0, 
                       VectorXd Q, VectorXcd zeta)
{
  int k;
  std::complex<double> s;
  s = Q0;
  for (k = 0; k < N; k++)
    s += Q(k) * log((z-zeta(k))/(z0-zeta(k)));
  return (z-z0) * std::exp(s);
}

int main(void)
{
  int i, j, k, N, m;
  double theta, pi, dt, R, Q0, r, dr;
  // 虚数単位と原点に移る点 z0=1/2
  std::complex<double> I(0,1), z0(0.6,0.0), zp, w, w2;

  R = 2;
  N = 80;
  VectorXcd zeta(N), z(N);
  VectorXd b(N), Q(N);
  MatrixXd a(N,N);

  pi = 4.0 * atan(1.0);
  dt = 2 * pi / N;
  // 円周 |z|=2 上の等分点
  for (k = 0; k < N; k++)
    zeta(k) = R * exp(I * (k * dt));
  //  std::cout << "zeta=" << zeta << std::endl;
  // 単位円周 |z|=1 上の等分点
  for (j = 0; j < N; j++)
    z(j) = exp(I * (j * dt));

  // 係数行列の準備
  for (j = 0; j < N; j++)
    for (k = 0; k < N; k++) {
      a(j,k) = log(abs(z(j)-zeta(k)));
    }
  // 右辺
  for (j = 0; j < N; j++) {
    b(j) = - log(abs(z(j)-z0));
  }
  // 電荷を求める
  Q = a.partialPivLu().solve(b);

  // Q0 を求める
  Q0 = 0;
  for (k = 0; k < N; k++) {
    Q0 += Q(k) * log(abs(z0-zeta(k)));
  }

  g_init("GRAPH", 150.0, 150.0);
  g_device(G_BOTH);
  g_def_scale(0, - 1.2, 1.2, -1.2, 1.2, 10.0, 10.0, 140.0, 140.0);
  g_sel_scale(0);

  m = 20;
  // 同心円の像
  dr = 1.0 / m;
  dt = 2 * pi / (4 * m);
  for (i = 1; i <= m; i++) {
    r = i * dr;
    for (j = 0; j <= 4 * m; j++) {
      zp = r * exp(I * (j * dt));
      w=f(zp,z0,N,Q0,Q,zeta);
      if (j == 0)
        g_move(w.real(), w.imag());
      else
        g_plot(w.real(), w.imag());
    }
  }
  // 放射線の像
  for (j = 0; j <= 4 * m; j++) {
    theta = j * dt;
    for (i = 0; i <= m; i++) {
      r = i * dr;
      zp = r * exp(I * theta);
      w=f(zp,z0,N,Q0,Q,zeta);
      if (i == 0)
        g_move(w.real(), w.imag());
      else
        g_plot(w.real(), w.imag());
    }
  }
  g_sleep(-1.0);
  g_term();
}

g++ -I/usr/X11R6/include conformalmap.cpp -L/usr/local/lib -lglscd -L/usr/X11R6/lib -lX11

図: $w=\dfrac{z-z_0}{1-\overline{z_0} z}$ , $z_0=1/2$