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);
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];
}
}
}