ストレンジ・アトラクタ

user-pic
0

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


すると、以下のように連立一次の常微分方程式であらわされ、ルンゲ・クッタ法で解ける。なお、この微分方程式にアトラクタを描かせるためには、定数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;
   }
}
//---------------------------------------------------------------------------

トラックバック(0)

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

コメントする

このブログ記事について

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

ひとつ前のブログ記事は「経験知」です。

次のブログ記事は「延命措置」です。

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