/*------------------------------------------------------------------------------------
ダフィング方程式から得られるジャパニーズ・アトラクタ
Japanese attractor derived from the Duffing equations
dx/dt = y
dy/dt = -K * y - x^3 + B * cos(t)
duffing_eqs.c
------------------------------------------------------------------------------------*/
#include "duffing_eqs.h"
int main(void)
{
int i;
int ii;
double t;
double x;
double y;
// 初期値の設定
// set the initial values
x = X0;
y = Y0;
t = 0.0;
// ４次のルンゲ・クッタ法を用いたダフィング方程式の計算
// calculate the Duffing equations using the 4th-order Runge-Kutta method
for( i = 0; i < SKIPS; i++ ){
for( ii = 0; ii < DIV; ii++ ){
duffing_RK4( t, &x, &y );
t += DELTA_T;
// printf("%.10lf %.10lf\n", x, y );
}
}
printf("#----------------------------------------------\n" );
printf("# ITERATION: %d\n", (int)ITERATION );
printf("# SKIPS: %d\n", (int)SKIPS );
printf("# DIV: %d\n", (int)DIV );
printf("# DELTA_T: %lf\n", (double)DELTA_T );
printf("# PHASE_DEGREE: %lf\n", (double)PHASE_DEGREE );
printf("#----------------------------------------------\n" );
printf("# x y\n" );
// ４次のルンゲ・クッタ法を用いたダフィング方程式の計算
// calculate the Duffing equations using the 4th-order Runge-Kutta method
for( i = 0; i < ITERATION; i++ ){
for( ii = 0; ii < DIV; ii++ ){
duffing_RK4( t, &x, &y );
t += DELTA_T;
// 全ての軌跡が欲しいとき
// if we want all trajectories
// printf("%.10lf %.10lf\n", x, y );
}
// ジャパニーズ・アトラクタ(２π毎のデータ)
// Japanese attractor (data for every 2 * PI)
printf("%.10lf %.10lf\n", x, y );
}
return 0;
}
/*-----------------------------------------------------------------------------------------*/
/*------------------------------------------
ダフィング方程式
Duffing equations
------------------------------------------*/
/* dx/dt = y */
double fx( double t, double x, double y )
{
double tt;
tt = t + 2.0 * M_PI * PHASE_DEGREE/(double)360.0;
return y + 0.0 * x + 0.0 * cos( tt );
}
/*------------------------------------------------
dy/dt = -K * y - x^3 + B * cos(t)
dy/dt = -K * y + x - x^3 + B * cos(t)
------------------------------------------------*/
double fy( double t, double x, double y )
{
double x3;
double tt;
x3 = x * x * x;
tt = t + 2.0 * M_PI * PHASE_DEGREE/(double)360.0;
return -K * y - x3 + B * cos( tt );
// return -K * y + x - x3 + B * cos(t);
}
/*-----------------------------------------------------------------------------------------*/
/*------------------------------------------
4次のルンゲ・クッタ法
4th-order Runge-Kutta Mehtod
------------------------------------------*/
void duffing_RK4( double t, double (*x), double (*y) )
{
double h;
double tt;
double kx1, kx2, kx3, kx4;
double ky1, ky2, ky3, ky4;
double xx, yy;
xx = (*x);
yy = (*y);
h = DELTA_T/(double)2.0;
/*--------------------*/
kx1 = DELTA_T * fx( t, xx, yy );
ky1 = DELTA_T * fy( t, xx, yy );
xx = (*x) + kx1/(double)2.0;
yy = (*y) + ky1/(double)2.0;
/*--------------------*/
tt = t + h;
kx2 = DELTA_T * fx( tt, xx, yy );
ky2 = DELTA_T * fy( tt, xx, yy );
xx = (*x) + kx2/(double)2.0;
yy = (*y) + ky2/(double)2.0;
/*--------------------*/
tt = t + h;
kx3 = DELTA_T * fx( tt, xx, yy );
ky3 = DELTA_T * fy( tt, xx, yy );
xx = (*x) + kx3;
yy = (*y) + ky3;
/*--------------------*/
tt = t + DELTA_T;
kx4 = DELTA_T * fx( tt, xx, yy );
ky4 = DELTA_T * fy( tt, xx, yy );
/*--------------------*/
xx = (*x) + ( kx1 + 2.0 * kx2 + 2.0 * kx3 + kx4 )/(double)6.0;
yy = (*y) + ( ky1 + 2.0 * ky2 + 2.0 * ky3 + ky4 )/(double)6.0;
(*x) = xx;
(*y) = yy;
}