数学の最近のブログ記事

私は数楽者なので、数式そのものがミーハー的に好きである。なので、たまたま外資系に勤務していたとき、数理関係の発明を扱っていたので、特許出願用に作成する必要があるときは、なるべく本格的な数式を書こうとした。

ここで本格的な数式というのは、数式の数学的な意味ではない。たかが数楽者にそんな高尚なことが分かるわけがない。数楽者に思いつくのは上にも書いたようなミーハーなことで要するに、数式の外見が本格的ということである。

本格的というのは数学の専門書みたいだということで、するとTeXで作成した数式、ということである。そこで私は、奥村晴彦著、「[改訂第5版] LaTeX2e 美文書作成入門」技術評論社を買って、付いてきたDVDでTeXを当時使っていたWindows 7のパソコンに導入した。

ちょっと記憶を辿って、私が特許出願で作成した数式のうち一番込み入っているのは、以下のものだと思う。

それで、Windows 10パソコンに移行したとき、LaTeX2e 美文書作成入門の最新版を買わなきゃ、と思った。

でも、Windows 10パソコンにプリインストールされてるWORD 2016をいじっているとき、普通に数式が作れることが分かった。それで作成した数式の例は、こちらである。

もちろんこの数式でもよくて、むしろ特許公報をみているとWORDで作成した数式の方が多いような気がする。そもそもTeXで作成した数式は、線が細くで見にくい。

しかしTeXで作成した数式の方がダンディーなのである。しかも私は、GUIでなく、スクリプトを書いて数式を表現したいのである。

だから私はやはりLaTeX2e 美文書作成入門の最新版を買うとしよう。ただ、細々とTeXのスクリプトを書くのも楽勝というわけではないのだが。

じつは、ストレンジ・アトラクタの微分方程式を数値的に描いてみたいとおもったモチベーションは、「カオスとフラクタル」講談社ブルーバックスでたまたま見かけた煙たなびくみたいな画像に一目ぼれしたからだが、いろいろ工夫してみてもどうしても意図する画像にならないでいろいろ調べているとき、ポアンカレ断面という用語を知った。

そこでポアンカレ断面とは何か調べてみると、どうも微分方程式の解の軌道内に超平面を配置し、その超平面と軌道の交点をプロットする、というものであるらしい。ストレンジ・アトラクタの場合、cos tというtに関する周期項があるので、t = 0, 2π, ... 2nπ,...を超平面として、軌道の交点をプロットしてみることにした。

ところがたまたま、以下のようなダフィンの方程式をウェブ上でみつけた。

ちなみにストレンジ・アトラクタの微分方程式を再掲すると、以下のとおりである。

この2つの微分方程式の最大の違いは、xという項の存在で、これは本質的である。ともかく、なんとなくであるが、ダフィンの方程式に関して、ポアンカレ断面を求めてみることにした。なおこのとき、ウェブ上情報の示唆により、δ=0.20、γ=0.30、ω=1.0ととった。

どうやってポアンカレ断面を求めるか、特に文献もウェブ上の情報もみつからなかったので、独自に考えてみた。

まず、角度θの象限を確認する。

  • 第1象限: cosθ > 0 且つsinθ > 0
  • 第2象限: cosθ < 0 且つsinθ > 0
  • 第3象限: cosθ < 0 且つsinθ < 0
  • 第4象限: cosθ > 0 且つsinθ < 0

そして、第4象限から第1象限に移るとき、あるいは第1象限から第4象限に移るとき、θ = 0の点を通るので、ポアンカレ断面としてプロットすることにした。その一例を示すと、以下のとおりである。

これは、関連の書籍でみるポアンカレ断面とだいたい似ているので、一応正しいとみなすことにした。そこで、ストレンジ・アトラクタの場合も同様にしてポアンカレ断面をプロットしてみると、以下のとおりである。ちなみに、このストレンジ・アトラクタでは、k = 1.0、B = 12.0としてみた。

最後にダフィンの方程式のためのプログラムのソースコードを以下に示す。ストレンジ・アトラクタの場合は、fny()として以下の関数を用い、

double fny(double x, double y, double t)  {
 return -0.1*y - x*x*x + 12.0*cos(t);
}
プロットのところを
 px = 300.0 + 28.0 * x;
 py = 320.0 - 28.0 * y;
とした。

ところで全く関係ないかもしれないが、ヒカシューの「二十世紀の終わりに」という曲を聞いていたら、以下のような歌詞があった。

二十世紀の終わりに恋をするなら
惑星のちからと死の魔術が必要
恋の断面図を透かしてみてごらん
そこは地獄かはたまた暗闇か

もしかしてポアンカレ断面は、恋の断面図なのか、なわけないよね、と意味不明のことをつい呟いてしまった。

//---------------------------------------------------------------------------

#include <vcl.h>
#include <stdlib.h>
#include <math.h>
#pragma hdrstop

#include "Main.h"
//---------------------------------------------------------------------------
#pragma package(smart_init)
#pragma resource "*.dfm"
TForm1 *Form1;
double fny(double x, double y, double t);
double ModT(double t);
int MODE(double th);
const double pi = 3.14156;
double d = 0.0;
//---------------------------------------------------------------------------
__fastcall TForm1::TForm1(TComponent* Owner)
	: TForm(Owner)
{
}
//---------------------------------------------------------------------------
void __fastcall TForm1::Close1Click(TObject *Sender)
{
 Close();
}
//---------------------------------------------------------------------------
void __fastcall TForm1::Draw1Click(TObject *Sender)
{
const double H = 0.001;
int i,j;
double x = 1.074,y = 3.021,t = 0.0;
double a1,a2,a3,a4;
double b1,b2,b3,b4;
double px,py;
double th;
int N = 0;
int n = 0;
int mode0 = 0, mode1 = 0;

 while (n < 100000 )  {
  N++;
  th = ModT(t);
  mode0 = MODE(th);
  a1 = y; 
  b1 = fny(x,y,t);
  a2 = y + H*b1/2.0;
  b2 = fny(x+H*a1/2.0,y+H*b1/2.0,t+H/2.0);
  a3 = y+H*b2/2.0;
  b3 = fny(x+H*a2/2.0,y+H*b2/2.0,t+H/2.0);
  a4 = y+H*b3;
  b4 = fny(x+H*a3,y+H*b3,t+H);
  x += H*(a1+2.0*a2+2.0*a3+a4)/6.0;
  y += H*(b1+2.0*b2+2.0*b3+b4)/6.0;
  t += H;
  if ( N > 1000 && (( mode1 == 3 && mode0 == 0 ) || ( mode1 == 0 && mode0 == 3 )) )     {
	  px = 400.0 + 140.0 * x;
	  py = 320.0 - 140.0 * y;
     Canvas->Pixels[px][py] = clBlack;
     n++;
  }
  mode1 = mode0;
 }
}

double fny(double x, double y, double t)  {
 return -0.20*y + x - x*x*x + 0.30*cos(t);
}

double ModT(double t)                     {
 double s = floor(t / (2.0 * pi));
 return t - 2.0*s*pi;
}

int MODE(double th)      {
 if ( cos(th + d) >= 0 && sin(th + d) >= 0 ) return 0;
 if ( cos(th + d) < 0 && sin(th + d) >= 0 ) return 1;
 if ( cos(th + d) < 0 && sin(th + d) < 0 ) return 2;
 if ( cos(th + d) >= 0 && sin(th + d) < 0 ) return 3;
 return 3;
}
//---------------------------------------------------------------------------

山口昌哉、「カオスとフラクタル」講談社ブルーバックスを更に読んでいて、アトラクタの挙動をするダフィンの方程式という常微分方程式が示されていたので、久々にルンゲ・クッタ法で解いてみることにした。その常微分方程式は以下のとおりである。


すると、以下のように連立一次の常微分方程式であらわされ、ルンゲ・クッタ法で解ける。なお、この微分方程式にアトラクタを描かせるためには、定数kとBをある特定の範囲の値に選ばなくてはならないが、ここではそのような値の一例として、k = 0.25, B = 8.0と選んだ。

この方程式は、時間tの項cos(t)を陽にもつのが特徴である。そこで、x,yにcos(t)を第3の軸として加え、三次元として表示させることにした。いろいろな側面の描画を便宜上異なる色で示した。各図からは、不規則なアトラクタになっていることが見て取れる。

先ず、x-yの軌跡は次のとおりである。

x-y-tの3Dの軌跡は次のとおりである。

x-tの軌跡は次のとおりである。

y-tの軌跡は次のとおりである。

よく分からないが、x-tの軌跡が、「カオスとフラクタル」講談社ブルーバックスにあった図に近いかもしれない。

//---------------------------------------------------------------------------

#include <vcl.h>
#include <stdlib.h>
#include <math.h>
#pragma hdrstop

#include "Main.h"
//---------------------------------------------------------------------------
#pragma package(smart_init)
#pragma resource "*.dfm"
TForm1 *Form1;
double fny(double x, double y, double t);
double ModT(double t);
const double pi = 3.14156;
int DT = 0;

//---------------------------------------------------------------------------
__fastcall TForm1::TForm1(TComponent* Owner)
	: TForm(Owner)
{
}
//---------------------------------------------------------------------------
void __fastcall TForm1::Close1Click(TObject *Sender)
{
  Close();
}
//---------------------------------------------------------------------------
void __fastcall TForm1::Draw1Click(TObject *Sender)
{
const double H = 0.00001;
int i,j;
double x = 1.074,y = 3.021,t = 0.0;
double a1,a2,a3,a4;
double b1,b2,b3,b4;
double px,py;
double TH1, TH2, COS1, COS2, SIN1, SIN2, th;
int N = 0;
int n = 1;

TH1 = pi/4.0;
TH2 = pi/4.0;

COS1 = cos(TH1);
SIN1 = sin(TH1);
COS2 = cos(TH2);
SIN2 = sin(TH2);
for ( i = 1; i < 1000; i++ )
 for ( j = 1; j < 760; j++ )
 Canvas->Pixels[i][j] = clWhite;

 while (N++ < 30000000 )  {
  a1 = y; 
  b1 = fny(x,y,t);
  a2 = y + H*b1/2.0;
  b2 = fny(x+H*a1/2.0,y+H*b1/2.0,t+H/2.0);
  a3 = y+H*b2/2.0;
  b3 = fny(x+H*a2/2.0,y+H*b2/2.0,t+H/2.0);
  a4 = y+H*b3;
  b4 = fny(x+H*a3,y+H*b3,t+H);
  x += H*(a1+2.0*a2+2.0*a3+a4)/6.0;
  y += H*(b1+2.0*b2+2.0*b3+b4)/6.0;
  t += H;
  if ( N > 10000 )               {
	th = ModT(t);
	if ( DT == 1 )                            {
	 px = 300.0 + 40.0*(COS1*x + SIN1*COS2*y + 2.0*SIN1*SIN2*cos(th));
	 py = 320.0 - 40.0*(-SIN1*x + COS1*COS2*y + 2.0*COS1*SIN2*cos(th));
	 Canvas->Pixels[px][py] = clNavy;
	} else if ( DT == 3 )                     {
	 px = 380.0 + 60.0*y;
	 py = 320.0 - 100.0*cos(th);
	 Canvas->Pixels[px][py] = clGreen;
	}  else if ( DT == 2 )                     {
	 px = 350.0 + 60.0*x;
	 py = 320.0 - 100.0*cos(th);
	 Canvas->Pixels[px][py] = clRed;
	} else                                     {
	 px = 300.0 + 40.0*x;
	 py = 320.0 - 40.0*y;
	 Canvas->Pixels[px][py] = clBlack;
	}
	n++;
  }
 }

  return;
}
//---------------------------------------------------------------------------

double fny(double x, double y, double t)  {
 return -0.25*y - x*x*x + 8.0*cos(t);
}

double ModT(double t)                     {
 int M = (int) t / (2.0 * pi);
 return (double) t - M * 2.0 * pi;
}
void __fastcall TForm1::RadioGroup1Click(TObject *Sender)
{
   switch(RadioGroup1->ItemIndex)   {
	 case 0:
	 DT = 0;
	 break;
	 case 1:
	 DT = 1;
	 break;
	 case 2:
	 DT = 2;
	 break;
	 default:
	 DT = 3;
	 break;
   }
}
//---------------------------------------------------------------------------

高木関数

user-pic
0

山口昌哉、「カオスとフラクタル」講談社ブルーバックスを更に眺めていると、高木関数というものがと言及されていた。すなわち、ここで高木とは、類体論で有名な高木貞治のことである。

高木貞治は、ワイエルシュトラスとは独立に、以下のようなアイデアによって、[0,1]の区間の至るところで連続であり、且つ至るところで微分不可能な関数T(x)を構成した。その数式を示すと以下のとおりである。

実は、「カオスとフラクタル」講談社ブルーバックスには高木関数の級数表示のところで誤記があり、そのままプログラミングしたら全くそれらしい関数が構成されなかった。そこで、ウェブで調べるとともに、自分でも考えて、以下のようにした。すなわち、φ(22x)のようにすると、もともと[0,1]だった定義域が広がってしまう。そこで、xが1より大きい場合は、その小数部分のみで関数の値を計算することにより、φ(22x)という値域の拡張を達成した。これで計算すると、想定されていた高木関数の表示が得られた。

折角なので、φ(x) = 4x(1-x)というワイエルシュトラス関数の場合も同様に計算して、併せて表示した。高木関数を黒、ワイエルシュトラス関数を赤のラインで示したのが以下のとおりである。

中心点だけ高木関数とワイエルシュトラス関数の値が一致し、その他の場合は高木関数 < ワイエルシュトラス関数となっていることが見て取れる。

でも待てよ、と考えてみた。私の拙い知識では、至るところ連続且つ微分不可能というワイエルシュトラスの関数は、三角関数を駆使したもので、φ(x) = 4x(1-x)のような関数形ではなかったような気がする。これはワイエルシュトラス関数のalternativeなのだろうか。機会があったら調べてみることにする。

また、改めて関数の連続性と微分可能性の定義を考えてみたが、すっかり忘れているようで、それを改めて勉強しなおして、高木関数の微分不可能性の証明も考えてみたい。

ともかく、例によって以下、プログラムのソースコードを示す。

//---------------------------------------------------------------------------

#include <vcl.h>
#include <math.h>
#pragma hdrstop

#include "Main.h"
//---------------------------------------------------------------------------
#pragma package(smart_init)
#pragma resource "*.dfm"
TForm1 *Form1;
double phai(double x);
double W(double x);
//---------------------------------------------------------------------------
__fastcall TForm1::TForm1(TComponent* Owner)
	: TForm(Owner)
{
}
//---------------------------------------------------------------------------
void __fastcall TForm1::Exit1Click(TObject *Sender)
{
 Close();
}
//---------------------------------------------------------------------------
void __fastcall TForm1::Draw1Click(TObject *Sender)
{
  double x,y,v,by,py,u,r;
  int i,j,bx,px;
  bx = 20;
  by = 660.0;
  Canvas->Pen->Color = clBlack;
  for ( i = 0; i < 700; i++ )               {
	x = (double) i / 700;
	u = x;
	r = 2.0;
	v = phai(u)/r;
	y = 0.0;
	for ( j = 0; j < 100; j++ )        {
	  y += v;
	  u *= 2.0;
	  r *= 2.0;
	  v = phai(u)/r;
	}
	px = i + 20;
	py = 660.0 - 700.0*y;
	Canvas->MoveTo(bx,by);
	Canvas->LineTo(px,py);
	bx = px; by = py;
  }
 // W-function 
    bx = 20;
  by = 660.0;
  Canvas->Pen->Color = clRed;
  for ( i = 0; i < 700; i++ )               {
	x = (double) i / 700;
	u = x;
	r = 2.0;
	v = W(u)/r;
	y = 0.0;
	for ( j = 0; j < 100; j++ )        {
	  y += v;
	  u *= 2.0;
	  r *= 2.0;
	  v = W(u)/r;
	}
	px = i + 20;
	py = 660.0 - 700.0*y;
	Canvas->MoveTo(bx,by);
	Canvas->LineTo(px,py);
	bx = px; by = py;
  }
}
//---------------------------------------------------------------------------

double phai(double x)   {
 double y;
 y = x - floor(x);
 if ( y >= 0 && y < 0.5 )  {
	return 2.0 * y;
 }	else if ( y >= 0.5 && y < 1.0 )  {
	return 2.0*(1.0 - y);
 }
 return 0.0;
}

double W(double x)      {
 double y;
 y = x - floor(x);
 return 4.0 * y * (1.0 - y);
}

Hidden Dimention Galleriesでは、以下のような一般化が与えられる。

すなわち、以下のような二元線形変換が定義される。

xn+1 = axn + byn + e
yn+1 = cxn + dyn + f

そして、これに確率pを組み合わせて、(a,b,c,d,e,f,p)という7つの実数の組で1つの変換が与えられる。

前記羊歯の場合、
(0.836,0.044,-0.044,0.836,0,0.169,1/2)
(-0.141,0.302,0.302,0.141,0,0.127,1/6)
(0.141,-0.302,0.302,0.141,0,0.169,1/6)
(0,0,0,0.175337,0,0,1/6)

であった。

上記ウェブサイトには、以下の例が与えられている。
(0.14 ,0.01 ,0.00 ,0.51 ,-0.08 ,-1.31 ,0.10)
(0.43 ,0.52 ,-0.45 ,0.50 ,1.49 ,-0.75 ,0.35)
(0.45 ,-0.49 ,0.47 ,0.47 ,-1.62 ,-0.74 ,0.35)
(0.49 ,0.00 ,0.00 ,0.51 ,0.02 ,1.62 ,0.20)

これで描いた図形が下記のとおりである。

よく分からないが、これはメープルの葉ということらしいが、現物をみたことがないので、ピンとこない。念のためソースコードも示す。

//---------------------------------------------------------------------------

#include <vcl.h>
#include <stdlib.h>
#pragma hdrstop

#include "Main.h"
//---------------------------------------------------------------------------
#pragma package(smart_init)
#pragma resource "*.dfm"
TForm1 *Form1;
//---------------------------------------------------------------------------
__fastcall TForm1::TForm1(TComponent* Owner)
	: TForm(Owner)
{
}
//---------------------------------------------------------------------------
void __fastcall TForm1::Close1Click(TObject *Sender)
{
  Close();
}
//---------------------------------------------------------------------------
void __fastcall TForm1::Draw1Click(TObject *Sender)
{
int i = 0, n, c;
double x,y,px,py,tx,ty;

// x, yに適当な初期値を与える
x = 0.1;
y = 0.1;

while ( i++ < 5000000 )   {

 n = random(136512)%100 + 1;

  if ( n >= 1 && n <= 10 )      {
     tx = 0.14*x + 0.01*y - 0.08;
     ty = 0.51*y - 1.31;
  } else if ( n >= 11 && n <= 45 )  {
     tx = 0.43*x + 0.52*y + 1.49;
     ty = -0.45*x + 0.50*y - 0.75;
  } else if ( n >= 46 && n <= 80 )  {
     tx = 0.45*x - 0.49*y - 1.62;
     ty = 0.47*x + 0.47*y - 0.74;
  } else                            {
     tx = 0.49*x + 0.02;
     ty = 0.51*y + 1.62;
  }

  if ( i > 100 )    {
	px = 350.0 + 70.0*tx;
	py = 300.0 - 70.0*ty;
	if ( i < 4950000 )
	Canvas->Pixels[px][py] = clGreen;
	else
	Canvas->Pixels[px][py] = clOlive;
  }
  x = tx;
  y = ty;
}
  return;
}
//---------------------------------------------------------------------------

以前に、羊歯のグラフィックというのを描いてみたことがあるが、これは要するに本にあるソースコードを打ち込んでみただけなので、アルゴリズムの輪郭もうまく把握できなかった。

最近たまたま、合原一幸編著「カオス カオス理論の基礎と応用」サイエンス社に、IFSのアトラクタとして、上記羊歯のグラフィックの元となると思われるアルゴリズムが説明してあったので、改めてこの説明から、プログラミングしてみることにした。

IFSとは、iterated function systemの略らしくて、日本語では反復関数系と訳される。なんともありきたりなネーミングだと苦笑するのであるが、米国の応用数学者M. Barnsleyが研究したものらしい。

Barnsleyは、羊歯を描くために、下記の4つの写像を用意した。
  1. W1(x,y) = (0.836x + 0.044y,-0.044x + 0.836y + 0.169)
  2. W2(x,y) = (-0.141x + 0.302y,0.302x + 0.141y + 0.127)
  3. W3(x,y) = (0.141x - 0.302y,0.302x + 0.141y + 0.169);
  4. W4(x,y) = (0.0,0.175337y);

そして、(x,y) ← W(x,y)の処理を繰り返すのだが、W(x,y)は実際は、W1(x,y),W2(x,y),W3(x,y),W4(x,y)のどれかである。その呼び出す確率がそれぞれ、1/2, 1/6, 1/6, 1/6である。こうして描いた図は、下記のとおりである。バリエーションとして、点の色を少し変えてみた。

なお、例によって、プログラムのソースコードも付ける。

//---------------------------------------------------------------------------

#include <vcl.h>
#include <stdlib.h>
#pragma hdrstop

#include "Main.h"
//---------------------------------------------------------------------------
#pragma package(smart_init)
#pragma resource "*.dfm"
TForm1 *Form1;
//---------------------------------------------------------------------------
__fastcall TForm1::TForm1(TComponent* Owner)
	: TForm(Owner)
{
}
//---------------------------------------------------------------------------
void __fastcall TForm1::Close1Click(TObject *Sender)
{
  Close();
}
//---------------------------------------------------------------------------
void __fastcall TForm1::Draw1Click(TObject *Sender)
{
int i = 0, n, c;
double x,y,px,py,tx,ty;

// x, yに適当な初期値を与える
x = 0.1;
y = 0.1;

while ( i++ < 30000000 )   {

 n = random(1000)%6 + 1;
 switch(n)  {
  case 4:
   tx = -0.141*x + 0.302*y;
   ty = 0.302*x + 0.141*y + 0.127;
   break;
  case 5:
   tx = 0.141*x - 0.302*y;
   ty = 0.302*x + 0.141*y + 0.169;	
   break;
  case 6:
   tx = 0.0;
   ty = 0.175337*y;	
   break;
  default:
   tx = 0.836*x + 0.044*y;
   ty = -0.044*x + 0.836*y + 0.169;
   break;
 }
  if ( i > 100 )    {
	px = 300.0 + 550.0*tx;
	py = 600.0 - 550.0*ty;
	c = i%6 + 1;
	if ( c == 4 )
    Canvas->Pixels[px][py] = clYellow;
	else if ( c == 5 )
    Canvas->Pixels[px][py] = clRed;
    else if ( c == 6 )
    Canvas->Pixels[px][py] = clBlue;
    else
    Canvas->Pixels[px][py] = clGreen;
  }
  x = tx;
  y = ty;
}
  return;
}
//---------------------------------------------------------------------------

フラクタルを一通り見たというわけではないが、典型的なカオスであるロジスティック写像f(x) = ax(1-x)について見てみる。

すなわち、x ← f(x)の繰り返しで、xの挙動をみる。このとき、初期値は、0と1の間の任意の実数でよいと思われるが、例えば0.5とする。このとき、ロジスティック写像の挙動は、aの値によって全く異なる。

0≦a≦3のとき、xは安定な周期点に収束する。面白いのは3<a<4のときで、2周期から2n周期を経て、一瞬3周期になり、それからカオスとなる。なお、3周期はカオスを意味するというのが、リー・ヨークが示唆したところである。

横軸をaの値、縦軸をxの値とし、3<a<4のグラフは以下のとおりとなる。

これを描いたプログラムのソースコードは、以下のとおりである。このプログラムではなぜか初期値を0.1ととり、100回ループを回して助走をつけてから、プロットするようにした。

//---------------------------------------------------------------------------

#include <vcl.h>
#pragma hdrstop

#include "Main.h"
//---------------------------------------------------------------------------
#pragma package(smart_init)
#pragma resource "*.dfm"
TForm1 *Form1;
//---------------------------------------------------------------------------
__fastcall TForm1::TForm1(TComponent* Owner)
	: TForm(Owner)
{
}
//---------------------------------------------------------------------------
void __fastcall TForm1::Button1Click(TObject *Sender)
{
int n, i;
double x, a, px, py;

 for ( n = 0; n < 500; n++ )   {
  x = 0.1;
  a = 3.0 + (double) n/500;
  for ( i = 0; i < 1000; i++ ) {
    x = a * x * (1.0 - x);
    px = 20 + n;
	py = 500.0 - 400.0 * x;
    if ( i > 100 ) Form1->Canvas->Pixels[px][py] = clBlack;
  }
 }
}
//---------------------------------------------------------------------------
void __fastcall TForm1::Button2Click(TObject *Sender)
{
 Close();
}
//---------------------------------------------------------------------------

完全数に続いて数論ネタだが、数学セミナーの6月号をみていたら、川田浩一氏の記事で、「3つの素数の3乗の和の話」というのがあったので、それがなに?と読んでみると、今年の西暦である2017が3つの素数の3乗の和であらわされる数だというので、そういう数がどれくらいあるのかちょっと計算してみた。

その結果は、西暦3000年までで、以下のとおりである。なお、この計算に使ったPerlプログラムを末尾に示す。下記で、赤字は結果が素数であることを示す。

  1. 0024 = 23 + 23 + 23
  2. 0043 = 23 + 23 + 33
  3. 0062 = 23 + 33 + 33
  4. 0081 = 33 + 33 + 33
  5. 0141 = 23 + 23 + 53
  6. 0160 = 23 + 33 + 53
  7. 0179 = 33 + 33 + 53
  8. 0258 = 23 + 53 + 53
  9. 0277 = 33 + 53 + 53
  10. 0359 = 23 + 23 + 73
  11. 0375 = 53 + 53 + 53
  12. 0378 = 23 + 33 + 73
  13. 0397 = 33 + 33 + 73
  14. 0476 = 23 + 53 + 73
  15. 0495 = 33 + 53 + 73
  16. 0593 = 53 + 53 + 73
  17. 0694 = 23 + 73 + 73
  18. 0713 = 33 + 73 + 73
  19. 0811 = 53 + 73 + 73
  20. 1029 = 73 + 73 + 73
  21. 1347 = 23 + 23 + 113
  22. 1366 = 23 + 33 + 113
  23. 1385 = 33 + 33 + 113
  24. 1464 = 23 + 53 + 113
  25. 1483 = 33 + 53 + 113
  26. 1581 = 53 + 53 + 113
  27. 1682 = 23 + 73 + 113
  28. 1701 = 33 + 73 + 113
  29. 1799 = 53 + 73 + 113
  30. 2017 = 73 + 73 + 113
  31. 2213 = 23 + 23 + 133
  32. 2232 = 23 + 33 + 133
  33. 2251 = 33 + 33 + 133
  34. 2330 = 23 + 53 + 133
  35. 2349 = 33 + 53 + 133
  36. 2447 = 53 + 53 + 133
  37. 2548 = 23 + 73 + 133
  38. 2567 = 33 + 73 + 133
  39. 2665 = 53 + 73 + 133
  40. 2670 = 23 + 113 + 113
  41. 2689 = 33 + 113 + 113
  42. 2787 = 53 + 113 + 113
  43. 2883 = 73 + 73 + 133

これからわかることは、3つの素数の3乗の和であらわされる数は、西暦が始まって2017年までで30回しかなく、年が素数であるという条件を加えると西暦が始まって2017年までで9回しかない、ということで相当稀である。

特に、2017年の直前で、3つの素数の3乗の和であらわされ且つ素数である年は1483年で、この年になにがあったかというと、世界史的には薔薇戦争と、マルチン・ルターの誕生くらいで、日本史的には室町時代で、特に特記すべき事件はないようである。

まあ、結果が素数という条件を外しても、2017年の直前で、3つの素数の3乗の和であらわされる年は、1799年だから、かなり昔である。この年は、エジプトにてロゼッタ・ストーンが発見され、ナポレオンが、ブリュメールのクーデタにてフランス政権を掌握した、というのだから、まあまあのエポック性である。日本史的には、寛政11年だが、特筆すべき出来事はないようである。こうみてくると、3つの素数の3乗の和であらわされる今年は、世界はともかく国内的にはあまり事件はないのだろう、ということになる。

ところで、川田浩一氏はこの記事で、3つの素数の3乗の和であらわされる数と、3つの素数の3乗の和であらわされ且つ素数である数がどれくらいの頻度なのか、π(x) ~ x/log(x)の素数定理を使って計算しているが、私がきちんと理解できていないので、ここでは述べない。

さらに記事に言及すると、素数の和によって自然数をあらわすことをGoldbach問題、自然数のべき乗の和によって自然数をあらわすことをWaring問題と呼び、ここから、素数のべき乗の和に関する問題は、Goldbach-Waring問題と呼ばれる、という。Goldbach問題だけでも天下の難問なのだから、これは相当難しそうである。

ちなみに、Goldbach-Waring問題の一つの到達点として、次の定理があるという。

十分大きい奇数は、9個の素数の3乗の和である。


push @PRIME, 2;

for ( $i = 3; $i < 20; $i++ )  {
 if (IsPrime($i) == 1 ) {
  push @PRIME, ,$i;
 }
}

foreach $x (@PRIME)    {
 foreach $y (@PRIME)    {
   foreach $z (@PRIME)    {
   $V = $x * $x * $x + $y * $y * $y + $z * $z * $z ;
    if ( $V < 3000 && $x > 0 && $x <= $y && $y <= $z )     {
     $K = sprintf "%04d",$V;
     $YY{$K} = $K;
     $YY{$K}->{S} = sprintf "%d<sup>3</sup> + %d<sup>3</sup> + %d<sup>3</sup>",$x,$y,$z;
    }
   }
 }
}

printf "<ol>\n";

foreach $x (sort keys %YY )    {
 if ( IsPrime($x) == 1 )    {
   printf "<li><font color=red>%04d = %s</font>\n",$x,$YY{$x}->{S};
 } else                     {
   printf "<li>%04d = %s\n",$x,$YY{$x}->{S};
 }
}

printf "</ol>\n";

sub IsPrime  {
my $n = shift;
my $i = 2;
  
  while ( $i * $i <= $n )  {
    if (($n % $i ) == 0 )  {
      return 0;
    }
    $i++;
  }

 return 1;
}

完全数

user-pic
0

藤井四段の29連勝の新記録について、今まで28連勝の記録をもっていた神谷八段が次のように発言したという。

28という完全数は一番好きな数字ですので、それが1位でなくなることは個人的に少々寂しいのですが、凡人がほぼ運だけで作った記録を天才が実力で抜いたというのは、将棋界にとってとてもいいことだと思います。藤井さんがこれからの数十年でどんな世界を見せてくれるのか、ファンの皆様とともに寿命の限り見続けていきたいです。

私は数論オタクでもあるが、素数はともかく完全数なんて意識したことはない。なので、神谷八段から「28という完全数は一番好きな数字ですので」と聞くのは、いかにも面白いと思った。そこで、完全数について改めて調べてみると、どうもメルセンヌ素数と関係あるらしくて、そこで公式に従い、改めて最初の数個の完全数を計算してみると以下のとおりである。

6
28
496
8128
33550336
8589869056
2305843008139952128

これによると、確かに28が二番目の完全数であることがわかる。また、三番目以降は現実的な連勝の数には該当しないこともすぐわかる。

なお、冗長なのでここでは計算しないが、メルセンヌ素数2p-1に対応する素数pは1000万桁以上まで求められているので、そのようなpと、多倍長計算ライブラリを使えば、1000桁以上の完全数を計算するのも容易である。

ちなみに計算に使ったPerlプログラムは以下のとおりである。

@P = (2,3,5,7,13,17,31);

foreach $n (@P)    {
 $x = 1;
 for ( $i = 0; $i < $n-1; $i++ ) {
   $x *= 2;
 }
 $y = $x*2 - 1;
 printf "%d<br>\n",$x*$y;
}

山口昌哉、「カオスとフラクタル」講談社ブルーバックスをきちんと読み直してみると、個別の図形の計算はポール・ピエール・レヴィだとしても、複数個の縮小写像による自己相似集合という表し方を発見したのは、日本では、京都大学理学部の畑 政義氏が最初ということである。現在、京都大学理学研究科/数学・数理解析専攻解析学講座/准教授という肩書である。

さて、レヴィの図形と同じ式で、定数a, b, c, dだけ違う他の図形も示しておく。

  1. a = 0, b = 0.4+0.5i, c = 0, d = 0.4+0.5i
  2. a = 0.4614+0.4614i, b = 0, c = 0.622-0.192i, d = 0
  3. a = 0, b = 0.5+0.2887i, c = 0, d = 0.6667
  4. a = 0.707i, b = 0, c = 0.5, d = 0
  5. a = 0.4614+0.4614i, b = 0, c = 0, d = 0.2896-0.585i

これを見ていると、1. は悪魔系、2.は雲、3.は木の枝、4.は氷柱、5.は海藻または細菌のコロニーと、なかなか表現力豊かである。そういえば、レヴィの図形は蟹か。

このアーカイブについて

このページには、過去に書かれたブログ記事のうち数学カテゴリに属しているものが含まれています。

前のカテゴリは教育です。

次のカテゴリは日常です。

最近のコンテンツはインデックスページで見られます。過去に書かれたものはアーカイブのページで見られます。