数学の最近のブログ記事

森肇教授の「カオス 流転する自然」岩波科学ライブラリーに、レスラー系のアトラクタというものが示されていたので、これも描いてみることにした。なお、レスラー系といっても、アントニオ猪木であるわけがなく、創案者はドイツ人のレスラー(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;
}

3体問題

user-pic
0

さて、2体問題に続き、いよいよ懸案の3体問題の記述に着手することにした。私が若い頃、ルンゲクッタ法で2体問題までは記述できていたのに、3体問題で挫折したのは、3体問題だと、記述が面倒になるからなのだった。それは難しいというよりも、たんに面倒というだけなのだが、当時使用していたコンピュータ言語はBASICで、関数を定義して呼び出すということさえできない。それで細かい記入ミスをしてしまい、私の散漫な注意力では書き間違いをみつけきれなかったようである。

2つの質点の間の引力を記述する関数は、fx(x1,y1,x2,y2)とfy(x1,y1,x2,y2)として定義したので、それをそのまま使えばよい。また、面倒な変数の置き換えは、Perlを使うことにした。

そこで例えば、下記のような内容のテキストファイルをRK.txtとして用意し、これに下記ソースコード1のPerlプログラムをかけることで、うまく変換でき、その結果のRKX.txtを利用することで、ソースコード2のC++のプログラムを得た。これは難なくコンパイルできた。ソースコード2をちらっとでも見ていただくと、練達のプログラマは除き、これを手書きだけで正しく書くことが若干煩わしいことは分かっていただけるだろう。こうして3体問題がなんとか記述できたのだが、適当に位置と速度を与えても、理解困難な不規則な軌道をとって、あっという間に発散してしまう。そこで、とりあえず原点を中心とした円に三角形を形成するように同じ質量の質点を配置し、円に対する接線の方向の初期速度を与えることで、下記のような安定解を得た。

これは、当然に昔から知られていて、ラグランジュの正三角形解と呼ばれているそうである。これは安定解といっても、実は針の先端で踊る天使のように?不安定で、一方の質点の質量を1%増やしただけでも、下記のように発散してしまう。

既知の3体問題の解はいろいろ知られているようであるが、ソースコードのせいで長くなったので、それは別稿とする。

RK.txt

 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);

ソースコード1

open FI, "<RK.txt" or die "zzz";
open FO, ">RKX.txt";

$S1{x1} = "x1+H*a1/2.0";
$S1{y1} = "y1+H*b1/2.0";
$S1{u1} = "u1+H*c1/2.0";
$S1{v1} = "v1+H*d1/2.0";
$S1{x2} = "x2+H*e1/2.0";
$S1{y2} = "y2+H*f1/2.0";
$S1{u2} = "u2+H*f1/2.0";
$S1{v2} = "v2+H*h1/2.0";
$S1{x3} = "x3+H*p1/2.0";
$S1{y3} = "y3+H*q1/2.0";
$S1{u3} = "u3+H*r1/2.0";
$S1{v3} = "v3+H*s1/2.0";

$S2{x1} = "x1+H*a2/2.0";
$S2{y1} = "y1+H*b2/2.0";
$S2{u1} = "u1+H*c2/2.0";
$S2{v1} = "v1+H*d2/2.0";
$S2{x2} = "x2+H*e2/2.0";
$S2{y2} = "y2+H*f2/2.0";
$S2{u2} = "u2+H*f2/2.0";
$S2{v2} = "v2+H*h2/2.0";
$S2{x3} = "x3+H*p2/2.0";
$S2{y3} = "y3+H*q2/2.0";
$S2{u3} = "u3+H*r2/2.0";
$S2{v3} = "v3+H*s2/2.0";

$S3{x1} = "x1+H*a3";
$S3{y1} = "y1+H*b3";
$S3{u1} = "u1+H*c3";
$S3{v1} = "v1+H*d3";
$S3{x2} = "x2+H*e3";
$S3{y2} = "y2+H*f3";
$S3{u2} = "u2+H*f3";
$S3{v2} = "v2+H*h3";
$S3{x3} = "x3+H*p3";
$S3{y3} = "y3+H*q3";
$S3{u3} = "u3+H*r3";
$S3{v3} = "v3+H*s3";

while (<FI>)     {

 printf FO "%s",$_;
}

seek FI,0,0;
$N = 2;
while (<FI>)    {
 chomp;
   foreach $x (sort keys %S1)  {
     $y = $S1{$x};
     $_ =~ s/$x/$y/g;
   }
 printf FO "%s%d%s\n",substr($_,0,2),$N,substr($_,3,120);
}

seek FI,0,0;
$N++;
while (<FI>)    {
 chomp;
   foreach $x (sort keys %S2)  {
     $y = $S2{$x};
     $_ =~ s/$x/$y/g;
   }
 printf FO "%s%d%s\n",substr($_,0,2),$N,substr($_,3,120);

}
seek FI,0,0;
$N++;
while (<FI>)    {
 chomp;
   foreach $x (sort keys %S3)  {
     $y = $S3{$x};
     $_ =~ s/$x/$y/g;
   }
 printf FO "%s%d%s\n",substr($_,0,2),$N,substr($_,3,120);

}

ソースコード2

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

#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 pi = 3.14159;
double TH= 2.0*pi/3.0;
int N = 0;           // 繰り返し回数
m1 = 1.0;
m2 = 1.0;
m3 = 1.0;


// 位置と速度の初期値
x1 = 3.0*cos(TH);  y1 = 3.0*sin(TH); u1 = -0.1*sin(TH); v1 = 0.1*cos(TH);
x2 = 3.0*cos(2.0*TH); y2 = 3.0*sin(2.0*TH); u2 = -0.1*sin(2.0*TH); v2 = 0.1*cos(2.0*TH);
x3 = 3.0; y3 = 0.0; u3 = 0.0; v3 = 0.1;

while ( N++ < 100000000 )        {

 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;
 px = 290.0 + 70.0*x1;
 py = 290.0 - 70.0*y1;
 Canvas->Pixels[px][py] = clRed;
 px = 290.0 + 70.0*x2;
 py = 290.0 - 70.0*y2;
 Canvas->Pixels[px][py] = clBlue;
 px = 290.0 + 70.0*x3;
 py = 290.0 - 70.0*y3;
 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;
}

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

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

本格的というのは数学の専門書みたいだということで、すると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;
}
//---------------------------------------------------------------------------

このアーカイブについて

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

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

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

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