数学の最近のブログ記事

2次方程式が整数解をもつ割合に続き、蛇足とは思ったが、3次方程式が整数解をもつ割合についてもみてみた。

単なる直感で、3次方程式が整数解をもつ割合は2次方程式が整数解をもつ割合より少ないかと思ったら、確かに3次方程式が整数解をもつ割合は2次方程式が整数解をもつ割合の0.64倍程度であるが、案外大きい。

ただ、107の試行で、2つ以上を解をもつ場合というのは5回しか観察できなかった。

あまり豊穣ではないかもしれないが、どうもよく分からない世界である。ともかくプログラムのソースコードをリストする。少なくとも、高次になるほど複数の解をもちにくくなるということはいえるかもしれない。

#include <stdio.h>
#include <stdlib.h>
#include <time.h>

long long int QD(int a, int b, int c, int d, int x);

int main()     {
 int i,j,a,b,c,d;
 int e,f;
 srand(time(NULL));
 
 for ( i = 0; i < 10000000; i++ )  {
	a = rand() % 1000 + 1;
	if (( rand() % 2 + 1) == 1 ) a = -a;
 	b = rand() % 1000 + 1;
	if (( rand() % 2 + 1) == 1 ) b = -b;
 	c = rand() % 1000 + 1;
	if (( rand() % 2 + 1) == 1 ) c = -c;
 	d = rand() % 1000 + 1;
	if (( rand() % 2 + 1) == 1 ) d = -d;
 	e = 0; // 解の数
 	if ( d < 0 )
 	  f = -d;
 	else
 	  f = d;
 	for ( j = -f; j <= f; j++ )    {
		if ( j != 0 && (f % j) == 0 ) {
		  if ( QD(a,b,c,d,j) == 0 )       {
			if ( e == 0 )
		      printf("a = %4d,b = %4d,c = %4d, d = %4d, 解1:%d\n",a,b,c,d,j);
		    else if ( e == 1 )
		      printf("a = %4d,b = %4d,c = %4d, d = %4d, 解2:%d\n",a,b,c,d,j);
		    else
		      printf("a = %4d,b = %4d,c = %4d, d = %4d, 解3:%d\n",a,b,c,d,j);
		     e++;
		    }
		}
	}
 }
 
 return 1;

}

long long int QD(int a, int b, int c, int d, int x)  {
 return a*x*x*x + b*x*x + c*x + d;
}

バックナンバーである数学セミナーの2017年10月号の「不可能性の証明」特集の金子元氏の「超越数論における不可能性」の記事のリュービルの定理や、eの超越性の証明を読んでいて、ううん実に難しいとうなっていたら、最後に参考文献のところで、西岡久美子「超越数とはなにか-代数方程式の解にならない数たちの話」講談社と、塩川宇賢「無理数と超越数」森北出版しかリストされていないので、唖然とした。

普通は、参考文献は少なくとも5つくらいはリストされるものなので、どうしても超越数論はマイナーなんだなと感じてしまった。とはいえ、超越数論は、ヒルベルト23の問題の七番目にも採り上げられているので、それなりに重要だと認識されているとはいえる。

私は西岡久美子さん(お笑いの好きな私は、どうしても、ピン芸人のにしおかすみこを思い出してしまう)のこの本はもっているのだが、塩川宇賢さんの本はもってないので、アマゾンの古本で急いで注文した。新刊では出ていないので入手できなくなることを懼れた。

そこで、塩川宇賢「無理数と超越数」森北出版のまえがきだけ読んだ段階で、単なる雑感のようなことを書いてみる。

さて、数というと、複素数も、四元数も、行列も入るのかもしれないが、ここでは素朴な意味での数を考える。すると、最初にくるのは、自然数と整数。これらが解として満たすのは、モニックな整係数の一次方程式である。

次に有理数となると、それが満たすのは、モニックでなくてよい、有理係数の一次方程式である。

次に無理数となると、有理係数の複数次方程式の解である場合と、有理係数の複数次方程式の解にはならない場合とがあり、前者は代数的数、後者は超越数と呼ばれる。

さて、問題は、特定の数が無理数であること、あるいは超越数であることの証明であり、それは場合により、途轍もなく難しい。

超越数であることが証明されている数として、e、π、2の√2乗、log102などがある。

ところが、e + πは、超越数どころか、無理数であるかどうかさえ分からないそうである。そんなもん超越数に決まっとるやんけ、みたいな思い込みは通用しない。

リーマンのゼータ関数に5を代入したζ(5)や、オイラー定数γも、超越数はおろか、無理数かどうかもわからないそうである。ところが、ζ(3)の無理数性を1978年にアペリが証明したというのだから、いよいよわからないものである。

で、ある命題を思いついた。それは、「超越数のうちで、超越数性が証明できるのは、高々加算無限である。」ううん、なんか当たり前で命題にならないか。ともかく、何億年経とうと、その超越数性が証明されない非加算無限の実数たちがひっそりと数の世界で暮らしている、ということだった。これもまた、いわゆる一つのプラトニズムではあるまいか。

2次方程式の整数解の続きで、ax2 + bx + c = 0が、整数解をもつ割合を調べてみようと思い立った。一応、a,b,cも1から999(負の数も含む)の整数と仮定した。

そのプログラミングは自明というか簡単であるが、例えば下記のようなものである。a,b,cを乱数で与えて、2次方程式の整数解 で示唆されたアルゴリズムに基づき、解があるかどうかみる。

すると、少なくとも1つの整数解をもつ確率は、だいたい1.3×10-3くらい、2つの整数解をもつ確率は1.0×10-6より少し小さいくらいだった。これはそんなに小さくないような気もする。

これを暗号化に使えないかと思案したが、なにしろ計算量が小さすぎて、到底堅牢な暗号にはなりえないような気がした。

ちなみに、解の一例は以下のとおりである。
a = -55,b = 494,c = 560, 解1:10
a = -215,b = 342,c = 557, 解1:-1
a = -789,b = 242,c = 547, 解1:1
a = -387,b = 756,c = 36, 解1:2
a = -387,b = -290,c = 677, 解1:1
a = -211,b = 395,c = 714, 解1:3
a = -4,b = -179,c = 45, 解1:-45
a = -729,b = -67,c = 662, 解1:-1
a = -573,b = 3,c = 570, 解1:1
a = -132,b = -162,c = 294, 解1:1
a = 404,b = 846,c = 76, 解1:-2
a = -269,b = -491,c = 948, 解1:-3
a = 38,b = 795,c = 757, 解1:-1
a = 78,b = -366,c = 288, 解1:1
a = -1,b = -20,c = 629, 解1:-37
a = -1,b = -20,c = 629, 解2:17


#include <stdio.h>
#include <stdlib.h>
#include <time.h>

int QD(int a, int b, int c, int x);

int main()     {
 int i,j,a,b,c;
 int d,f;
 srand(time(NULL));
 
 for ( i = 0; i < 10000000; i++ )  {
	a = rand() % 1000 + 1;
	if (( rand() % 2 + 1) == 1 ) a = -a;
 	b = rand() % 1000 + 1;
	if (( rand() % 2 + 1) == 1 ) b = -b;
 	c = rand() % 1000 + 1;
 	f = c;
 	if (( rand() % 2 + 1) == 1 ) c = -c;	
 	d = 0; // 解の数
 	for ( j = -f; j <= f; j++ )    {
		if ( j != 0 && (f % j) == 0 ) {
		  if ( QD(a,b,c,j) == 0 )       {
			if ( d == 0 )
		      printf("a = %4d,b = %4d,c = %4d, 解1:%d\n",a,b,c,j);
		    else
		      printf("a = %4d,b = %4d,c = %4d, 解2:%d\n",a,b,c,j);
		     d++;
		    }
		}
	}
 }
 
 return 1;

}

int QD(int a, int b, int c, int x)  {
 return a*x*x + b*x + c;
}

精霊の導きにより?レスラー系のアトラクタにポアンカレ断面を入れて、そこに現れた下記のような軌跡があまりにも特異なので驚いたのだが、まだその不思議さの全貌をみていないような気がして、もう少し拡張してみた。

すなわち、前回はa = b = 0.2で、cを2.0から4.3まで小刻みに増分していったのだが、cが4.0以降も綾があるような気がして、cを4.0から5.6まで増分してみることにした。また、折角なので、以下のようにプログラムのソースコードを変えた。

  1. y ≧ 0以上の領域だけをみることにして、表示の縦方向の縮尺を大きくした。
  2. x2 + y2の平方根でなく、xの値だけで縦方向をプロットすることにより計算を簡略化した。

すると、以下のとおり、蜘蛛の巣状からカオスになったあと思うと突然、霧が晴れたような長い期間があり、あれっカオスはもうないのかと思ったら、c = 5.6の直前でまたカオスが始まるという、なんとも不思議な性質を呈するのだった。

そこで次に、c = 5.0から6.0までを見てみたら、また新たな側面が現れた。ここまでくると普通は発散するか、カオス性が消えうせたりするものだが、レスラー系のアトラクタのポアンカレ断面はなかなか非自明な様相を呈して面白いのだった。

さらに、c = 5.8から6.8までを見てみたら、更にまた新たな様相を呈した。ここでよく観察すると、複数の直線と曲線が入り交じった構造がみてとれ、最後の方で、希釈されたカオスとでもいうべき微妙な領域があらわれる。この相次ぐ異なる貌は、まるで怪人二十面相であり、まさに覆面レスラーと呼びたいが、どうなんだろう。

懲りないで、c = 6.7から7.7までみたら、これはカオスではないのか、ともかくまた濃い領域があらわれた。

さらなる変化を求めて、c = 7.7から8.7までみたら、さらに変化があらわれたが、きりがないのでこのくらいにしておく。

ともかく、一応ソースコードを掲げる。

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

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

#include "Main.h"
//---------------------------------------------------------------------------
#pragma package(smart_init)
#pragma resource "*.dfm"
TForm1 *Form1;
double fnx(double y,double z); 
double fny(double a,double x,double y);
double fnz(double b,double c,double x,double z);

//---------------------------------------------------------------------------
__fastcall TForm1::TForm1(TComponent* Owner)
	: TForm(Owner)
{
}
//---------------------------------------------------------------------------
void __fastcall TForm1::Close1Click(TObject *Sender)
{
 Close();
}
//---------------------------------------------------------------------------
void __fastcall TForm1::Draw1Click(TObject *Sender)
{
const double H = 0.0001;
const pi = 3.14156;
double x,y,z;
double a = 0.2, b = 0.2, c;
double a1,a2,a3,a4;
double b1,b2,b3,b4;
double c1,c2,c3,c4;
double px, py;
int p_sgn, c_sgn;
int i,N;

for ( i = 0; i < 800; i++ )      {

 N = 0;
 c = 4.0 + 1.6*i/800.0;
 x = 2.0; y = 2.0; z = 13.0;
 if ( x + y >= 0 )
 p_sgn = c_sgn = 1;
else
 p_sgn = c_sgn = -1;
 
 while (N++ < 50000000 )  {
  a1 = fnx(y,z); 
  b1 = fny(a,x,y);
  c1 = fnz(b,c,x,z);
  a2 = fnx(y+H*b1/2.0,z+H*c1/2.0);
  b2 = fny(a,x+H*a1/2.0,y+H*b1/2.0);
  c2 = fnz(b,c,x+H*a1/2.0,z+H*c1/2.0);
  a3 = fnx(y+H*b2/2.0,z+H*c2/2.0);
  b3 = fny(a,x+H*a2/2.0,y+H*b2/2.0);
  c3 = fnz(b,c,x+H*a2/2.0,z+H*c2/2.0);
  a4 = fnx(y+H*b3,z+H*c3);
  b4 = fny(a,x+H*a3,y+H*b3);
  c4 = fnz(b,c,x+H*a3,z+H*c3);
  x += H*(a1+2.0*a2+2.0*a3+a4)/6.0;
  y += H*(b1+2.0*b2+2.0*b3+b4)/6.0;
  z += H*(c1+2.0*c2+2.0*c3+c4)/6.0;
  if ( x + y >= 0 )
   c_sgn = 1;
   else
   c_sgn = -1;
  if ( N > 10000 )                {
   if ( p_sgn != c_sgn )            {
     px = 30 + i;
     if ( y >= 0 )
	   py = 200.0 - 60.0*x;

     Canvas->Pixels[px][py] = clBlack;
   }
  }
  p_sgn = c_sgn;
 } // end of while
} // end of for 
}
//---------------------------------------------------------------------------
double fnx(double y,double z)  {
 return -y-z;
}

double fny(double a,double x,double y)  {
 return x + a*y;
}

double fnz(double b,double c,double x,double z)  {
 return b + x*z - c*z;
}

レスラー系のアトラクタを描いてみたものの、どうもあまりピンとこない。あまりカオスぽくないのである。それで次第にレスラー系のアトラクタのことは興味が薄れていった。

それからしばらくして、ある休日の朝、近所のスーパーに買い物に行こうと思い立った。意外かもしれないが、私はスーパーに行くのが好きである。いろいろな商品を眺めていると楽しくなる。すると家内が、卵も買ってきて、と頼むので、メモに書いて、胸ポケットに入れた。そうしないと最近、すぐに忘れてしまう。

そしてスーパーから買い物して帰ってきたはいいが、家内が「卵は?」と聞く。「えっ、買ったよ。」と応えた。レシートをみると確かに卵はリストされているのだが、確かに買い物袋には卵が入ってない。そこで記憶を辿ると、どうも卵を、買い物籠から買い物袋に移すときにそこに置き忘れてしまったようである。すなわち、卵は割れやすいので、最後に買い物袋入れようとおもって、いったん別のところに置き、そのままになってしまったようである。

そこで仕方なく、スーパーに戻ることにした。スーパーに行き、レジの女性に、「あのう、卵を忘れたみたいなんですけど。」と言ったら、サービスカウンタに聞いてみて下さい、という。そこでサービスカウンタにも、レシートを提示しつつ、「卵を忘れたみたいなんですけど。」といったら、ちゃんと卵がそこにキープされていて、あっさり渡してもらえてほっとした。私が忘れた卵をラッキーと持ち去る人がいたら仕方ないと思ったが、そういう人は案外いないものである。

で、ほっとして卵をもって帰りながら、空を眺めていたら急に、「レスラー系にポアンカレ断面を入れなさい。」と声が聞こえてきたような気がした。えっ、どういうこと?と考えてみた。もちろん、具体的な示唆はなにも与えられない。

で、自宅でコーヒーを飲みながら、1時間くらい思案して、以下のとおり、なんとか方針が決まった。
  1. レスラー系は、x,y,zの3次元のシステムだが、簡便のため、x,yの2次元でみる。
  2. ポアンカレ断面は、y = -xの、原点を通る、勾配-1の直線とし、この直線と常微分方程式の解の交点をみる。
  3. 解の交点は、原点からのユークリッド距離でプロットする。でも後で考えたら、勾配-1の直線との交点なんだから、単にxまたはyの座標値でもよかった。
  4. 一応、y ≧ 0と、そうでない場合とで分けて、前者は座標の正値に、後者は座標の負値にプロットする。
  5. a = b = 0.2で、cを2.0から4.3まで小刻みに増やして、それぞれでロジスティック写像のようにプロットする。

するとなんとまあ、下記のように、まるで蜘蛛の巣のような軌跡が現れて、なんだこれはとしばらく茫然となるほどだった。

そうか、精霊?はこれを私に見せたくて、「レスラー系にポアンカレ断面を入れなさい。」と言霊を降らしたのだな、と思ってしまった次第である。ともかく例によって、プログラムのソースコードを示す。
//---------------------------------------------------------------------------

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

#include "Main.h"
//---------------------------------------------------------------------------
#pragma package(smart_init)
#pragma resource "*.dfm"
TForm1 *Form1;
double fnx(double y,double z); 
double fny(double a,double x,double y);
double fnz(double b,double c,double x,double z);
double ROTx(double TH1, double TH2, double x, double y, double z);
double ROTy(double TH1, double TH2, double x, double y, double z);
//---------------------------------------------------------------------------
__fastcall TForm1::TForm1(TComponent* Owner)
	: TForm(Owner)
{
}
//---------------------------------------------------------------------------
void __fastcall TForm1::Close1Click(TObject *Sender)
{
 Close();
}
//---------------------------------------------------------------------------
void __fastcall TForm1::Draw1Click(TObject *Sender)
{
const double H = 0.0001;
const pi = 3.14156;
double x,y,z;
double a = 0.2, b = 0.2, c;
double a1,a2,a3,a4;
double b1,b2,b3,b4;
double c1,c2,c3,c4;
double px, py;
int p_sgn, c_sgn;
int i,N;

for ( i = 0; i < 800; i++ )      {

 N = 0;
 c = 2.0 + 2.4*i/800.0;
 x = 2.0; y = 2.0; z = 13.0;
 if ( x + y >= 0 )
 p_sgn = c_sgn = 1;
else
 p_sgn = c_sgn = -1;
 
 while (N++ < 50000000 )  {
  a1 = fnx(y,z); 
  b1 = fny(a,x,y);
  c1 = fnz(b,c,x,z);
  a2 = fnx(y+H*b1/2.0,z+H*c1/2.0);
  b2 = fny(a,x+H*a1/2.0,y+H*b1/2.0);
  c2 = fnz(b,c,x+H*a1/2.0,z+H*c1/2.0);
  a3 = fnx(y+H*b2/2.0,z+H*c2/2.0);
  b3 = fny(a,x+H*a2/2.0,y+H*b2/2.0);
  c3 = fnz(b,c,x+H*a2/2.0,z+H*c2/2.0);
  a4 = fnx(y+H*b3,z+H*c3);
  b4 = fny(a,x+H*a3,y+H*b3);
  c4 = fnz(b,c,x+H*a3,z+H*c3);
  x += H*(a1+2.0*a2+2.0*a3+a4)/6.0;
  y += H*(b1+2.0*b2+2.0*b3+b4)/6.0;
  z += H*(c1+2.0*c2+2.0*c3+c4)/6.0;
  if ( x + y >= 0 )
   c_sgn = 1;
   else
   c_sgn = -1;
  if ( N > 10000 )                {
   if ( p_sgn != c_sgn )            {
     px = 30 + i;
     if ( y >= 0 )
	   py = 300.0 - 30.0*sqrt(x*x + y*y);
     else
	   py = 300.0 + 30.0*sqrt(x*x + y*y);
     Canvas->Pixels[px][py] = clBlack;
   }
  }
  p_sgn = c_sgn;
 } // end of while
} // end of for 
}
//---------------------------------------------------------------------------
double fnx(double y,double z)  {
 return -y-z;
}

double fny(double a,double x,double y)  {
 return x + a*y;
}

double fnz(double b,double c,double x,double z)  {
 return b + x*z - c*z;
}

double ROTx(double TH1, double TH2, double x, double y, double z)  {
 return cos(TH1)*cos(TH2)*x + sin(TH1)*y + cos(TH1)*sin(TH2)*z;
}

double ROTy(double TH1, double TH2, double x, double y, double z)  {
 return -sin(TH1)*cos(TH2)*x + cos(TH1)*y - sin(TH1)*sin(TH2)*z;
}

数学セミナーの2017年10月号の「不可能性の証明」特集の菊池誠氏の記事の最初の方に、ディオファントス方程式の簡単な例として、2次方程式の整数解のことが書かれてあった。

すなわち、a, b, cを0でない整数として、ax2 + bx + c = 0の整数解を求める問題であるが、cの約数を求めて、それをax2 + bx + cのxに代入して、ゼロになっていれば解である、というナイーブな方法でよいという。なお、1という約数も使うので、cは素数であってもよい。

確かに、ある整数n,mがあって、ax2 + bx + c = (ax - n)(x-m) と因数分解されるなら、 c = nmでmが解になっている。もしnがaで割り切れるなら、n' = n/aと置いてn'も解であり、nm = an'mである。

これは2次方程式、すなわち1元2次方程式の整数解を求めるアルゴリズムだったが、ちょっと考えると、1元n次方程式の整数解を求める場合にもそのまま適用できることがわかる。

ところが、もっと一般にn個の未知数を含む整数係数の多項式が整数解を持つか否かを有限に判定する方法をみつけよ、というのがヒルベルトの第10問題である。

そんなの無理に決まってんじゃん、というのが常識的判断であるが、これをきちんと証明しようとすると、とてつもなく難しい。

それでも、50年後の1950年には最終解決には至らないもののデービス、パットナム、ロビンソンによる本質的な貢献があり、それを継承してマチャセビッチが1970年に、ヒルベルトの第10問題を否定的に解決した。

また、ヒルベルトの第10問題は、ゲーデルの不完全性定理とも関係するようであるが、難しすぎて私にはまったく理解できない。

森肇教授の「カオス 流転する自然」岩波科学ライブラリーに、レスラー系のアトラクタというものが示されていたので、これも描いてみることにした。なお、レスラー系といっても、アントニオ猪木であるわけがなく、創案者はドイツ人のレスラー(Roessler)という人らしい。"oe"は本当はウムラウトなのだが、それは往々にして文字化けするのでここでは使わない。その方程式は、Wikiも併せて参照すると、以下のとおりである。

ここで定数a = 0.2、b = 0.2、c = 5.7でルンゲ・クッタで描いた軌跡が以下のとおりである。

なんでも、cが特定の値より大きいときカオスになるらしくて、そこでどのあたりが臨界値か調べてみた。

定数a = 0.2、b = 0.2、c = 4.7でも以下のとおり定常軌道である。

定数a = 0.2、b = 0.2、c = 2.0だと以下のとおり定常軌道である。

定数a = 0.2、b = 0.2、c = 3.0だと以下のとおり定常軌道である。。

定数a = 0.2、b = 0.2、c = 4.1だと以下のとおり定常軌道である。。

定数a = 0.2、b = 0.2、c = 4.2だと以下のとおり定常軌道かもしれないが微妙である。よって、4.2より少し大きいところに臨界値があると思われる。

プログラムのソースコードは、以下のとおりである。

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

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

#include "Main.h"
//---------------------------------------------------------------------------
#pragma package(smart_init)
#pragma resource "*.dfm"
TForm1 *Form1;
double fnx(double y,double z); 
double fny(double a,double x,double y);
double fnz(double b,double c,double x,double z);
double ROTx(double TH1, double TH2, double x, double y, double z);
double ROTy(double TH1, double TH2, double x, double y, double z);
//---------------------------------------------------------------------------
__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;
const pi = 3.14156;
double x = 2.0,y = 2.0,z = 13.0;
double a = 0.2, b = 0.2, c = 4.25;
double a1,a2,a3,a4;
double b1,b2,b3,b4;
double c1,c2,c3,c4;
double px, py;
int N = 0;

 while (N++ < 50000000 )  {
  a1 = fnx(y,z); 
  b1 = fny(a,x,y);
  c1 = fnz(b,c,x,z);
  a2 = fnx(y+H*b1/2.0,z+H*c1/2.0);
  b2 = fny(a,x+H*a1/2.0,y+H*b1/2.0);
  c2 = fnz(b,c,x+H*a1/2.0,z+H*c1/2.0);
  a3 = fnx(y+H*b2/2.0,z+H*c2/2.0);
  b3 = fny(a,x+H*a2/2.0,y+H*b2/2.0);
  c3 = fnz(b,c,x+H*a2/2.0,z+H*c2/2.0);
  a4 = fnx(y+H*b3,z+H*c3);
  b4 = fny(a,x+H*a3,y+H*b3);
  c4 = fnz(b,c,x+H*a3,z+H*c3);
  x += H*(a1+2.0*a2+2.0*a3+a4)/6.0;
  y += H*(b1+2.0*b2+2.0*b3+b4)/6.0;
  z += H*(c1+2.0*c2+2.0*c3+c4)/6.0;
  if ( N > 100 )                {
   px = 400.0 + 15.0*ROTx(pi/4.0, pi/4.0,x,y,z);
   py = 400.0 + 15.0*ROTy(pi/4.0, pi/4.0,x,y,z);
   Canvas->Pixels[px][py] = clBlack;
  }
 }
}
//---------------------------------------------------------------------------
double fnx(double y,double z)  {
 return -y-z;
}

double fny(double a,double x,double y)  {
 return x + a*y;
}

double fnz(double b,double c,double x,double z)  {
 return b + x*z - c*z;
}

double ROTx(double TH1, double TH2, double x, double y, double z)  {
 return cos(TH1)*cos(TH2)*x + sin(TH1)*y + cos(TH1)*sin(TH2)*z;
}

double ROTy(double TH1, double TH2, double x, double y, double z)  {
 return -sin(TH1)*cos(TH2)*x + cos(TH1)*y - sin(TH1)*sin(TH2)*z;
}

エノンのアトラクタというかエノン写像の続きだが、Wikiによると、ロジスティック写像のように分岐した後カオスになる図が存在するようなので、描いてみることにした。下記のとおりである。

この図は、下記のようなエノン写像において、b = 0.3と固定し、aを1.0から1.42まで次第に増やしていく。
xi+1 = 1 - axi2 + yi
yi+1 = byi

すなわち、図において横軸はaであり、縦軸はxの値である。aが1.42で打ち切られているのは、それよりaを増やすと式が発散してしまうからである。

この図からは、ロジスティック写像よりも複雑な分岐の様子がわかる。なお図において、縦軸と横軸にそれぞれ0.1刻みで目盛りをつけてみた。

なお参考までに、縦軸をyにしたときの図は、こちらである。

見て取れるように、yの方も分岐はあるのだが、値の幅が相対的に小さくなっている。

最後に例によって、プログラムのソースコードを付ける。

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

#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::Close1Click(TObject *Sender)
{
 Close();
}
//---------------------------------------------------------------------------
void __fastcall TForm1::Draw1Click(TObject *Sender)
{
 double x,y,a,b,tx,ty,px,py;
 int N = 0;
 int i;
 b = 0.3;
 for ( i = 0; i < 422; i++ )      {
   a = 1.0 + (double) 0.5*i/500.0; // パラメータaの値のセット
   x = y = 0.1;  // 初期値
   N = 0;
   while ( N++ < 800 )     {
     tx = 1.0 - a*x*x + y;
     ty = b*x;
	 x = tx; y = ty;
	 if ( N > 100 )          {
	   px = 100 + i;
	   py = 250.0 - 150.0 * x;
	   Canvas->Pixels[px][py] = clBlack;
	 }
  }
 }
 // ここで座標線を引く
 Canvas->Pen->Color = clRed;
 Canvas->MoveTo(1,250); Canvas->LineTo(550,250);
 Canvas->MoveTo(100,2); Canvas->LineTo(100,600);
 for ( i = 200; i <= 500; i += 100 )    {
	Canvas->MoveTo(i,25); Canvas->LineTo(i,475);
 }
 for ( i = 1; i <= 15; i++ )             {
   Canvas->MoveTo(80,250+i*15); Canvas->LineTo(120,250+i*15);
 }
 for ( i = 1; i <= 15; i++ )             {
   Canvas->MoveTo(80,250-i*15); Canvas->LineTo(120,250-i*15);
 }
}
//---------------------------------------------------------------------------

カオスを勉強しているうちに、大学時代に授業に出たことがある森肇教授が書かれた啓蒙書があることを知った。それは、「カオス 流転する自然」岩波科学ライブラリーである。

大学時代、森教授は統計力学の授業を担当していて、大変偉い方だという評判だった。なんでも非平衡の統計力学が専門だと伺ったのだが、当時不真面目な学生だった私はほとんど授業に出なかった。勿体ないというか申し訳ないので、せめてもの償いにこの本を買って勉強させていただくことにした。

その本を読んでいると、エノン写像というのが出てきた。なんでも、以下のような式らしい。
xi+1 = 1 - 1.4xi2 + 0.3yi
yi+1 = yi

ただ、ここにはあまり詳しいことが載ってないので、Wikiで確認して、微妙に違う以下の式でアトラクタを描いてみた。
xi+1 = 1 - 1.4xi2 + yi
yi+1 = 0.3yi

すると、下記のようになかなか格好いい、ナイキの商標のような軌跡になるのだった。

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

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

#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::Close1Click(TObject *Sender)
{
 Close();
}
//---------------------------------------------------------------------------
void __fastcall TForm1::Draw1Click(TObject *Sender)
{
 double x,y,a,b,tx,ty,px,py;
 int N = 0;
 int i,j;
 b = 0.3;
 a = 1.4;
   x = y = 0.1;  // 初期値
   while ( N++ < 10000 )     {
     tx = 1.0 - a*x*x + y;
     ty = b*x;
     x = tx; y = ty;
     if ( N > 100 )         {
       px = 300.0 + 150*x;
       py = 300.0 - 300.0*y;
       Canvas->Pixels[px][py] = clBlack;
    }
  }

}
//---------------------------------------------------------------------------

なんとか3体問題の面白そうな軌道をみてみたいとウェブで検索していて、中村 健という方のこちらのページをみつけて、あまりにもいろいろなほしい情報が開示されているので感激した。そこで、中村さんに感謝しつつ、この情報を使わせていただくことにした。

まず、ここで教わったことは重心座標を使う、ということである。勘の悪い私は、重心座標といわれても、は?としばらく意味が分からなかったが、ああ、3体の重心を座標の原点とみなして描画するのだな、とやっとわかった。確かに、私が試みたある3体の軌道では、3体が一緒になってx座標が増える方向に移動して、とうとう視界からみえなくなってしまった。こういうのも重心座標だと、きちんと追える。

すなわち、下記のように、まず重心(xc,yc)を求め、x1 - xcのような計算で、重心からの相対的な距離をプロットする。ここでLは表示のスケーリングであり、この値を適当に調整して、見やすいように、且つ画面からはみださないように表示される軌道のサイズを決める。

 xc = (m1*x1 + m2*x2 + m3*x3)/(m1+m2+m3);
 yc = (m1*y1 + m2*y2 + m3*y3)/(m1+m2+m3);
 px = 320.0 + L*(x1 - xc);
 py = 290.0 - L*(y1 - yc);
 Canvas->Pixels[px][py] = clRed;
 px = 320.0 + L*(x2 - xc);
 py = 290.0 - L*(y2 - yc);
 Canvas->Pixels[px][py] = clBlue;
 px = 320.0 + L*(x3 - xc);
 py = 290.0 - L*(y3 - yc);
 Canvas->Pixels[px][py] = clPurple;

そこで、中村さんのホームページから気まぐれで、いくつかピックアップして描画してみた。

まず、面白そうだなと思ったのが、周期解の近くとあった、この軌道。単なる気分で、もとの質量がm1 = m2 = m3 = 3だったのを2.5にしてみたら、変な二層構造になってしまった。この場合は、L = 30とした。

質点の番号mxyuv
1(赤)2.5-300-0.5
2(青)2.53000.5
3(紫)2.505.190.50

次は単純にみえる8の字解。この場合、L = 150とした。単純なのは結果であって、これは今世紀に入って、C. Moore, A. Chenciner & R. Montgomery, C. Sim´o らによって発見されたもので、そのパラメータを決定するために相当な努力があった。

質点の番号mxyuv
1(赤)1/3-1000
2(青)1/31000.5
3(紫)1/300-0.60.923

次は、中村さんがBroucke解と呼ぶもの。私が示した軌道は、早めに打ち切ったものである。理由は、途中から計算誤差のせいだと思うが、突然軌道が発散してしまったからである。それに対処するためにHの刻み値を小さくしたら、軌道を描くのが低速になりすぎて、うんざりした。なお、L = 130とした。

質点の番号mxyuv
1(赤)1/3-1000.3407
2(青)1/31000.3407
3(紫)1/3000.1-0.6814

次は、基本に戻って、新周期解の一種。L = 60とした。

質点の番号mxyuv
1(赤)3-300-0.5
2(青)33000.5
3(紫)305.180.4460

ところで、3体問題の軌道を描いてみての感想は、一見安定しているようでいて、わずかな計算誤差などが、軌道の崩壊につながる、ということである。そのため、私が使っている4次のルンゲクッタ法で、刻みの値Hとして極めて小さいものを使ったりしているのだが、するとどんどん軌道の進み速度が低下する。これはやはり、8次のルンゲクッタ法を使う必要があるのか、と思いつつ、何を調べればそれが載っているのかわからず、困っている。

別のこととして、パラメータの決定がある。たとえば、レヴィの図形ジュリア集合IFSのアトラクタのどれにおいても、適切なパラメータを使用しないと、描画がほとんど意味をなさなくなるので、私はウェブや文献などから学んだ値をほぼそのまま使っているのだが、こういうのをオリジナルに決定するのって、きっとすごく大変なんだよね、と自分に言い聞かせるだけで、どうするのかまったく想像もつかない。なんか情けないが能力のなさはどうしようもない。ともかく、重複をおそれず、ソースコードもリストしておく。

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

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

#include "Main.h"
//---------------------------------------------------------------------------
#pragma package(smart_init)
#pragma resource "*.dfm"
TForm1 *Form1;
double fx(double x1, double y1, double x2, double y2);
double fy(double x1, double y1, double x2, double y2);
//---------------------------------------------------------------------------
__fastcall TForm1::TForm1(TComponent* Owner)
	: TForm(Owner)
{
}
//---------------------------------------------------------------------------
void __fastcall TForm1::Close1Click(TObject *Sender)
{

 Close();
}
//---------------------------------------------------------------------------
void __fastcall TForm1::Run1Click(TObject *Sender)
{
const double H = 0.000001;  // 刻み
double x1,y1,u1,v1;
double x2,y2,u2,v2;
double x3,y3,u3,v3;
double a1,a2,a3,a4;  // dx1
double b1,b2,b3,b4;  // dy1
double c1,c2,c3,c4;  // du1
double d1,d2,d3,d4;  // dv1
double e1,e2,e3,e4;  // dx2
double f1,f2,f3,f4;  // dy2
double g1,g2,g3,g4;  // du2
double h1,h2,h3,h4;  // dv2
double p1,p2,p3,p4;  // dx3
double q1,q2,q3,q4;  // dy3
double r1,r2,r3,r4;  // du3
double s1,s2,s3,s4;  // dv3
double m1,m2,m3;       // 質量
double px, py;
double xc, yc;  //重心
double pi = 3.14159;
double TH= 2.0*pi/3.0;
double L = 60.0;
int N = 0;           // 繰り返し回数
m1 = 3.0;
m2 = 3.0;
m3 = 3.0;


// 位置と速度の初期値
x1 = -3.0;  y1 = 0.0; u1 = 0.0; v1 = -0.5;
x2 = 3.0; y2 = 0.0; u2 = 0.0; v2 = 0.5;
x3 = 0.0; y3 = 5.18; u3 = 0.446; v3 = 0.0;

while ( N++ < 1000000000 )        {

 a1 = u1;
 b1 = v1;
 c1 = -m2*fx(x1,y1,x2,y2) - m3*fx(x1,y1,x3,y3);
 d1 = -m2*fy(x1,y1,x2,y2) - m3*fy(x1,y1,x3,y3);
 e1 = u2;
 f1 = v2;
 g1 = -m1*fx(x2,y2,x1,y1) - m3*fx(x2,y2,x3,y3);
 h1 = -m1*fy(x2,y2,x1,y1) - m3*fy(x2,y2,x3,y3);
 p1 = u3;
 q1 = v3;
 r1 = -m1*fx(x3,y3,x1,y1) - m2*fx(x3,y3,x2,y2);
 s1 = -m1*fy(x3,y3,x1,y1) - m2*fy(x3,y3,x2,y2);
 a2 = u1+H*c1/2.0;
 b2 = v1+H*d1/2.0;
 c2 = -m2*fx(x1+H*a1/2.0,y1+H*b1/2.0,x2+H*e1/2.0,y2+H*f1/2.0) 
     - m3*fx(x1+H*a1/2.0,y1+H*b1/2.0,x3+H*p1/2.0,y3+H*q1/2.0);
 d2 = -m2*fy(x1+H*a1/2.0,y1+H*b1/2.0,x2+H*e1/2.0,y2+H*f1/2.0) 
      - m3*fy(x1+H*a1/2.0,y1+H*b1/2.0,x3+H*p1/2.0,y3+H*q1/2.0);
 e2 = u2+H*f1/2.0;
 f2 = v2+H*h1/2.0;
 g2 = -m1*fx(x2+H*e1/2.0,y2+H*f1/2.0,x1+H*a1/2.0,y1+H*b1/2.0) 
      - m3*fx(x2+H*e1/2.0,y2+H*f1/2.0,x3+H*p1/2.0,y3+H*q1/2.0);
 h2 = -m1*fy(x2+H*e1/2.0,y2+H*f1/2.0,x1+H*a1/2.0,y1+H*b1/2.0) 
      - m3*fy(x2+H*e1/2.0,y2+H*f1/2.0,x3+H*p1/2.0,y3+H*q1/2.0);
 p2 = u3+H*r1/2.0;
 q2 = v3+H*s1/2.0;
 r2 = -m1*fx(x3+H*p1/2.0,y3+H*q1/2.0,x1+H*a1/2.0,y1+H*b1/2.0) 
      - m2*fx(x3+H*p1/2.0,y3+H*q1/2.0,x2+H*e1/2.0,y2+H*f1/2.0);
 s2 = -m1*fy(x3+H*p1/2.0,y3+H*q1/2.0,x1+H*a1/2.0,y1+H*b1/2.0) 
      - m2*fy(x3+H*p1/2.0,y3+H*q1/2.0,x2+H*e1/2.0,y2+H*f1/2.0);
 a3 = u1+H*c2/2.0;
 b3 = v1+H*d2/2.0;
 c3 = -m2*fx(x1+H*a2/2.0,y1+H*b2/2.0,x2+H*e2/2.0,y2+H*f2/2.0) 
      - m3*fx(x1+H*a2/2.0,y1+H*b2/2.0,x3+H*p2/2.0,y3+H*q2/2.0);
 d3 = -m2*fy(x1+H*a2/2.0,y1+H*b2/2.0,x2+H*e2/2.0,y2+H*f2/2.0) 
      - m3*fy(x1+H*a2/2.0,y1+H*b2/2.0,x3+H*p2/2.0,y3+H*q2/2.0);
 e3 = u2+H*f2/2.0;
 f3 = v2+H*h2/2.0;
 g3 = -m1*fx(x2+H*e2/2.0,y2+H*f2/2.0,x1+H*a2/2.0,y1+H*b2/2.0) 
      - m3*fx(x2+H*e2/2.0,y2+H*f2/2.0,x3+H*p2/2.0,y3+H*q2/2.0);
 h3 = -m1*fy(x2+H*e2/2.0,y2+H*f2/2.0,x1+H*a2/2.0,y1+H*b2/2.0) 
      - m3*fy(x2+H*e2/2.0,y2+H*f2/2.0,x3+H*p2/2.0,y3+H*q2/2.0);
 p3 = u3+H*r2/2.0;
 q3 = v3+H*s2/2.0;
 r3 = -m1*fx(x3+H*p2/2.0,y3+H*q2/2.0,x1+H*a2/2.0,y1+H*b2/2.0) 
      - m2*fx(x3+H*p2/2.0,y3+H*q2/2.0,x2+H*e2/2.0,y2+H*f2/2.0);
 s3 = -m1*fy(x3+H*p2/2.0,y3+H*q2/2.0,x1+H*a2/2.0,y1+H*b2/2.0) 
      - m2*fy(x3+H*p2/2.0,y3+H*q2/2.0,x2+H*e2/2.0,y2+H*f2/2.0);
 a4 = u1+H*c3;
 b4 = v1+H*d3;
 c4 = -m2*fx(x1+H*a3,y1+H*b3,x2+H*e3,y2+H*f3) 
      - m3*fx(x1+H*a3,y1+H*b3,x3+H*p3,y3+H*q3);
 d4 = -m2*fy(x1+H*a3,y1+H*b3,x2+H*e3,y2+H*f3) 
      - m3*fy(x1+H*a3,y1+H*b3,x3+H*p3,y3+H*q3);
 e4 = u2+H*f3;
 f4 = v2+H*h3;
 g4 = -m1*fx(x2+H*e3,y2+H*f3,x1+H*a3,y1+H*b3) 
      - m3*fx(x2+H*e3,y2+H*f3,x3+H*p3,y3+H*q3);
 h4 = -m1*fy(x2+H*e3,y2+H*f3,x1+H*a3,y1+H*b3) 
      - m3*fy(x2+H*e3,y2+H*f3,x3+H*p3,y3+H*q3);
 p4 = u3+H*r3;
 q4 = v3+H*s3;
 r4 = -m1*fx(x3+H*p3,y3+H*q3,x1+H*a3,y1+H*b3) 
      - m2*fx(x3+H*p3,y3+H*q3,x2+H*e3,y2+H*f3);
 s4 = -m1*fy(x3+H*p3,y3+H*q3,x1+H*a3,y1+H*b3) 
      - m2*fy(x3+H*p3,y3+H*q3,x2+H*e3,y2+H*f3);


 x1 += H*(a1 + 2.0*a2 + 2.0*a3 + a4)/6.0;
 y1 += H*(b1 + 2.0*b2 + 2.0*b3 + b4)/6.0;
 u1 += H*(c1 + 2.0*c2 + 2.0*c3 + c4)/6.0;
 v1 += H*(d1 + 2.0*d2 + 2.0*d3 + d4)/6.0;
 x2 += H*(e1 + 2.0*e2 + 2.0*e3 + e4)/6.0;
 y2 += H*(f1 + 2.0*f2 + 2.0*f3 + f4)/6.0;
 u2 += H*(g1 + 2.0*g2 + 2.0*g3 + g4)/6.0;
 v2 += H*(h1 + 2.0*h2 + 2.0*h3 + h4)/6.0;
 x3 += H*(p1 + 2.0*p2 + 2.0*p3 + p4)/6.0;
 y3 += H*(q1 + 2.0*q2 + 2.0*q3 + q4)/6.0;
 u3 += H*(r1 + 2.0*r2 + 2.0*r3 + r4)/6.0;
 v3 += H*(s1 + 2.0*s2 + 2.0*s3 + s4)/6.0;
 
 xc = (m1*x1 + m2*x2 + m3*x3)/(m1+m2+m3);
 yc = (m1*y1 + m2*y2 + m3*y3)/(m1+m2+m3);
 px = 320.0 + L*(x1 - xc);
 py = 290.0 - L*(y1 - yc);
 Canvas->Pixels[px][py] = clRed;
 px = 320.0 + L*(x2 - xc);
 py = 290.0 - L*(y2 - yc);
 Canvas->Pixels[px][py] = clBlue;
 px = 320.0 + L*(x3 - xc);
 py = 290.0 - L*(y3 - yc);
 Canvas->Pixels[px][py] = clPurple;
}
}
//---------------------------------------------------------------------------
double fx(double x1, double y1, double x2, double y2) {
 double r,R;
 r = (x1 - x2)*(x1 - x2) + (y1 - y2)*(y1 - y2);
 R = r * sqrt(r);
 return (x1 - x2)/R;
}

double fy(double x1, double y1, double x2, double y2) {
 double r,R;
 r = (x1 - x2)*(x1 - x2) + (y1 - y2)*(y1 - y2);
 R = r * sqrt(r);
 return (y1 - y2)/R;
}

このアーカイブについて

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

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

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

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