/*------------------------------------------------------------------------------
ダブルスクロール方程式とヤコビ行列（線形化したダブルスクロール方程式）
double_scroll equations and the Jacobian (linearized double_scroll equations)
solve_double_scroll_eqs.c
------------------------------------------------------------------------------*/
#include "double_scroll_lyap.h"
void double_scroll_eqs( double v[DIM], double dfu[DIM][DIM], double dvdt[DIM], double df_k[DIM][DIM] )
{
int i;
int j;
int k;
double d2;
double xx;
double yy;
double zz;
double temp;
double df[DIM][DIM];
xx = v[0]; /* xx = Vc1 */
yy = v[1]; /* yy = Vc2 */
zz = v[2]; /* zz = iL */
/*------------------------------------------------------------------*/
// ダブルスクロール方程式
// double scroll equations
temp = G * (yy - xx) - func_g( xx );
dvdt[0] = temp/(double)C1;
temp = G * (xx - yy) + zz;
dvdt[1] = temp/(double)C2;
dvdt[2] = -yy/(double)L;
/*------------------------------------------------------------------*/
// ヤコビ行列
// Jacobian matrix
temp = -G - d_func_g( xx );
df[0][0] = temp/(double)C1;
df[0][1] = G/(double)C1;
df[0][2] = 0.0;
df[1][0] = G/(double)C2;
df[1][1] = -G/(double)C2;
df[1][2] = 1.0/(double)C2;
df[2][0] = 0.0;
df[2][1] = -1.0/(double)L;
df[2][2] = 0.0;
/*
｜ df[0][0] df[0][1] df[0][2] ｜
｜ ｜
df = ｜ df[1][0] df[1][1] df[1][2] ｜
｜ ｜
｜ df[2][0] df[2][1] df[2][2] ｜
*/
// ヤコビ行列を用いて微小変位ベクトルを計算
// 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] * dfu[k][j];
}
df_k[i][j] = d2;
}
}
}
double func_g( double xx )
{
double fx;
fx = M0 * xx + 0.5 * (M1 - M0) * fabs( xx + Bp ) + 0.5 * (M0 - M1) * fabs( xx - Bp );
return fx;
}
double d_func_g( double xx )
{
double temp;
temp = M0 + 0.5 * (M1 - M0) * ( xx + Bp )/(double)fabs( xx + Bp )
+ 0.5 * (M0 - M1) * ( xx - Bp )/(double)fabs( xx - Bp );
return temp;
}
/*-----------------------------------------------------------------------------------------*/
void solve_double_scroll_eqs_using_RK4( double v[DIM], double u[DIM][DIM] )
{
int i;
int j;
double vv[DIM];
double v_k1[DIM], v_k2[DIM], v_k3[DIM], v_k4[DIM];
double dfu[DIM][DIM];
double df_k1[DIM][DIM], df_k2[DIM][DIM], df_k3[DIM][DIM], df_k4[DIM][DIM];
/*******************************************************/
for (i = 0; i < DIM; i++) {
vv[i] = v[i];
for (j = 0; j < DIM; j++) {
dfu[i][j] = u[i][j];
}
}
/*******************************************************/
// vvとdfuを渡し、v_k1とdf_k1の計算値を得る。
// pass vv and dfu, and get the calculated values of v_k1 and df_k1
double_scroll_eqs( vv, dfu, v_k1, df_k1 );
for (i = 0; i < DIM; i++) {
vv[i] = v[i] + DELTA_T * v_k1[i]/(double)2.0;
}
for (i = 0; i < DIM; i++) {
for (j = 0; j < DIM; j++) {
dfu[i][j] = u[i][j] + DELTA_T * df_k1[i][j]/(double)2.0;
}
}
/*******************************************************/
double_scroll_eqs( vv, dfu, v_k2, df_k2 );
for (i = 0; i < DIM; i++) {
vv[i] = v[i] + DELTA_T * v_k2[i]/(double)2.0;
}
for (i = 0; i < DIM; i++) {
for (j = 0; j < DIM; j++) {
dfu[i][j] = u[i][j] + DELTA_T * df_k2[i][j]/(double)2.0;
}
}
/*******************************************************/
double_scroll_eqs( vv, dfu, v_k3, df_k3 );
for (i = 0; i < DIM; i++) {
vv[i] = v[i] + DELTA_T * v_k3[i];
}
for (i = 0; i < DIM; i++) {
for (j = 0; j < DIM; j++) {
dfu[i][j] = u[i][j] + DELTA_T * df_k3[i][j];
}
}
/*******************************************************/
double_scroll_eqs( vv, dfu, v_k4, df_k4 );
for (i = 0; i < DIM; i++) {
v[i] += DELTA_T * (v_k1[i] + 2.0 * v_k2[i] + 2.0 * v_k3[i] + v_k4[i])/(double)6.0;
}
for (i = 0; i < DIM; i++) {
for (j = 0; j < DIM; j++) {
u[i][j] += DELTA_T * (df_k1[i][j] + 2.0 * df_k2[i][j] + 2.0 * df_k3[i][j] + df_k4[i][j])/(double)6.0;
}
}
/*******************************************************/
}