solve_mackey-glass_DDE.c

/*------------------------------------------------------------------------

マッキー・グラス遅延微分方程式と線形化した方程式

Mackey-Glass delay differential equation and linearized equation

dx/dt = A * x(t-tau)/( 1 + x(t-tau)^C ) - B * x(t)

solve_mackey-glass_DDE.c

------------------------------------------------------------------------*/

#include "mackey-glass_DDE_lyap.h"

double uu[DIM+1][DIM];

void mackey_glass_DDE( double x, double delay, double (*dxdt) )

{

(*dxdt) = A * delay/(double)( 1.0 + pow( delay, C ) ) - B * x;

}

void solve_mackey_glass_DDE_using_RK4( double v[DIM], double (*dxdt), double u[DIM][DIM] )

{

int i;

int j;

double d1;

double d2;

double dudt;

double vv;

double delay;

double v_k1, v_k2, v_k3, v_k4;

vv = v[0];

delay = v[X_TAU];

/*******************************************************/

// vvとdelayを渡し、v_k1の計算値を得る。

// pass vv and delay, and get the calculated value of v_k1

mackey_glass_DDE( vv, delay, &v_k1 );

vv = v[0] + DELTA_T * v_k1/(double)2.0;

mackey_glass_DDE( vv, delay, &v_k2 );

vv = v[0] + DELTA_T * v_k2/(double)2.0;

mackey_glass_DDE( vv, delay, &v_k3 );

vv = v[0] + DELTA_T * v_k3;

mackey_glass_DDE( vv, delay, &v_k4 );

(*dxdt) = (v_k1 + 2.0 * v_k2 + 2.0 * v_k3 + v_k4)/(double)6.0;

vv = v[0] + DELTA_T * (*dxdt);

for ( i = 1; i < DIM; i++ ) {

v[DIM-i] = v[DIM-i-1];

}

v[0] = vv;

/*******************************************************/

// 線形化したマッキー・グラス遅延微分方程式を用いて微小変位ベクトルを計算

// calculate the micro-displacement vectors using linearized Mackey-Glass DDE

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

d1 = pow( delay, C ); /* x(t-tau)^C */

d2 = (1.0 + d1) * (1.0 + d1);

dudt =  u[DIM-1][i] * ( ( A * ( 1.0 + d1 ) - A * C * d1 )/(double)d2 ) - B * u[0][i];

uu[0][i] = u[0][i] + DELTA_T * dudt;

}

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

for ( j = 0; j < DIM; j++ ) {

uu[i+1][j] = u[i][j];

u[i][j] = uu[i][j];

}

}

}