solve_ikeda_eqs.c

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

池田写像とヤコビ行列（線形化した池田写像）

Ikeda map and the Jacobian matrix (linearized ikeda map)

solve_ikeda_eqs.c

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

#include "ikeda_lyap.h"

void ikeda_eqs( double (*x), double (*y), double u[DIM][DIM] )

{

int i, j, k;

double xx;

double yy;

double x2;

double y2;

double xy;

double theta;

double d2;

double temp1;

double temp2;

double q1, q2, q3, q4;

double df[DIM][DIM];

double dfu[DIM][DIM];

// 池田写像

// Ikeda map

x2 = (*x) * (*x);

y2 = (*y) * (*y);

theta = A - B/(double)( 1.0 + x2 + y2 );

xx = 1.0 + MYU * ( (*x) * cos(theta) - (*y) * sin(theta) );

yy = MYU * ( (*x) * sin(theta) + (*y) * cos(theta) );

(*x) = xx;

(*y) = yy;

// ヤコビ行列

// Jacobian matrix

x2 = (*x) * (*x);

y2 = (*y) * (*y);

xy = (*x) * (*y);

theta = A - B/(double)( 1.0 + x2 + y2 );

temp1 = 1.0 + x2 + y2;

temp2 = temp1 * temp1;

q1 = 1.0 - ( 2.0 * B * xy )/(double)temp2;

q2 = ( 2.0 * B * x2 )/(double)temp2;

q3 = 1.0 + ( 2.0 * B * xy )/(double)temp2;

q4 = ( 2.0 * B * y2 )/(double)temp2;

df[0][0] = MYU * (  q1 * cos(theta) - q2 * sin(theta) );

df[0][1] = MYU * ( -q3 * sin(theta) - q4 * cos(theta) );

df[1][0] = MYU * (  q1 * sin(theta) + q2 * cos(theta) );

df[1][1] = MYU * (  q3 * cos(theta) - q4 * sin(theta) );

/*

｜ df[0][0]  df[0][1] ｜

df = ｜ ｜

｜ df[1][0]  df[1][1] ｜

*/

// ヤコビ行列を用いて微小変位ベクトルを計算

// calculate the micro-displacement vector using the Jacobian matrix

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

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

d2 = 0.0;

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

d2 += df[i][k] * u[k][j];

}

dfu[i][j] = d2;

}

}

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

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

u[i][j] = dfu[i][j];

}

}

}