8.5.1 complexjacobi.cpp

/*
 * complexjacobi.cpp --- 与えられた k に対し sn(z;k), cn(z;k), dn(z;k)
 *   g++ -I/opt/local/include complexjacobi.cpp
 */

#include <iostream>
#include <iomanip>
#include <cmath>
#include <boost/math/special_functions.hpp> // 面倒なので一気に

using namespace std;
using namespace boost::math;

// M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions, §16.21
// 複素数版 Jacobi の sn(z;k)
complex<double> jacobi_sn(double k, complex<double> z)
{
  double kp,m,x,y,s,c,d,s1,c1,d1,numerator;
  x = z.real(); y = z.imag(); kp = sqrt(1 - k * k); m = k * k;
  s  = jacobi_sn(k,  x); c  = jacobi_cn(k,  x); d  = jacobi_dn(k,  x);
  s1 = jacobi_sn(kp, y); c1 = jacobi_cn(kp, y); d1 = jacobi_dn(kp, y);
  numerator = c1 * c1 + m * s * s * s1 * s1;
  return complex<double>(s * d1 / numerator, c * d * s1 * c1 / numerator);
}

// 複素数版 Jacobi の sn(), cn(), dn()
complex<double> jacobi_elliptic(double k, complex<double> z,
                                complex<double> *cn, complex<double> *dn)
{
  double kp,m,x,y,s,c,d,s1,c1,d1,numerator;
  complex<double> sn;
  x = z.real(); y = z.imag(); kp = sqrt(1 - k * k); m = k * k;
  s  = jacobi_sn(k,  x); c  = jacobi_cn(k,  x); d  = jacobi_dn(k,  x);
  s1 = jacobi_sn(kp, y); c1 = jacobi_cn(kp, y); d1 = jacobi_dn(kp, y);
  numerator = c1 * c1 + m * s * s * s1 * s1;
  sn  = complex<double>(s * d1 / numerator,       c * d * s1 * c1 / numerator);
  *cn = complex<double>(c * c1 / numerator,     - s * d * s1 * d1 / numerator);
  *dn = complex<double>(d * c1 * d1 / numerator, - m * s * c * s1 / numerator);
  return sn;
}

int main(int argc, char **argv)
{
  double k;
  complex<double> z, sn, cn, dn;

  z = complex<double>(0.5, 0.5);
  cout << "z=" << z << endl;
  cout << "Jacobi sn" << endl;
  // Mathematica で験算
  // z = 1/2 + I/2;
  // Table[JacobiSN[z,k^2], {k,0,0.9,0.1}]
  for (k = 0; k <= 0.9; k += 0.1) {
    sn = jacobi_sn(k, z);
    cout << "k=" << k << ", sn=" << sn << endl;
  }

  cout << "Jacobi sn,cn,dn" << endl;
  // Mathematica で験算
  // z = 1/2 + I/2
  // Table[{JacobiSN[z,k^2], JacobiCN[z,k^2], JacobiDN[z,k^2]}, {k,0,0.9,0.1}]
  for (k = 0; k <= 0.9; k += 0.1) {
    sn = jacobi_elliptic(k, z, &cn, &dn);
    cout << "k=" << k
         << ", sn=" << sn << ", cn=" << cn << ", dn=" << dn << endl;
  }
  return 0;
}

簡単なテスト (Mathematica の結果と照らし合わせて一応パス)

$ ./complexjacobi
z=(0.5,0.5)
Jacobi sn
k=0, sn=(0.540613,0.457304)
k=0.1, sn=(0.540868,0.45676)
k=0.2, sn=(0.54163,0.455127)
k=0.3, sn=(0.542892,0.45241)
k=0.4, sn=(0.54464,0.448615)
k=0.5, sn=(0.546858,0.443751)
k=0.6, sn=(0.549523,0.437829)
k=0.7, sn=(0.552609,0.430864)
k=0.8, sn=(0.556086,0.422873)
k=0.9, sn=(0.559922,0.413876)
Jacobi sn,cn,dn
k=0, sn=(0.540613,0.457304), cn=(0.989585,-0.249826), dn=(1,-0)
k=0.1, sn=(0.540868,0.45676), cn=(0.989175,-0.24975), dn=(0.999583,-0.00247149)
k=0.2, sn=(0.54163,0.455127), cn=(0.987946,-0.249518), dn=(0.998323,-0.00987698)
k=0.3, sn=(0.542892,0.45241), cn=(0.985903,-0.249121), dn=(0.996186,-0.0221895)
k=0.4, sn=(0.54464,0.448615), cn=(0.983055,-0.248545), dn=(0.993121,-0.0393642)
k=0.5, sn=(0.546858,0.443751), cn=(0.979413,-0.247769), dn=(0.989054,-0.0613386)
k=0.6, sn=(0.549523,0.437829), cn=(0.974994,-0.246768), dn=(0.983895,-0.0880327)
k=0.7, sn=(0.552609,0.430864), cn=(0.969816,-0.24551), dn=(0.977535,-0.11935)
k=0.8, sn=(0.556086,0.422873), cn=(0.963902,-0.24396), dn=(0.969854,-0.155176)
k=0.9, sn=(0.559922,0.413876), cn=(0.957279,-0.24208), dn=(0.960717,-0.195384)
$