ポアンカレ断面

user-pic
0

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

そこでポアンカレ断面とは何か調べてみると、どうも微分方程式の解の軌道内に超平面を配置し、その超平面と軌道の交点をプロットする、というものであるらしい。ストレンジ・アトラクタの場合、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;
}
//---------------------------------------------------------------------------

トラックバック(0)

トラックバックURL: http://commutative.world.coocan.jp/mt/mt-tb.cgi/4621

コメントする

このブログ記事について

このページは、あやたろうが2017年8月 6日 05:41に書いたブログ記事です。

ひとつ前のブログ記事は「フラクタルの応用」です。

次のブログ記事は「ボーダーライン」です。

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