発振器

周期的な外力 F sin(wt)　または F cos(wt) なんてのが加わっている N 次元の非自律系を
sin(wt) → x ，cos(wt) → y と置き換えることで，N+2 次元の自律系に書き直せる．

x, y の時間 t（横軸）に対する関係．（b=1）

例として，次式で与えられる Duffing 方程式を扱う．（N = 2）

これを次のように 4 (= N + 2) 次元の自律系にして計算したものを図に示す．

Duffing 方程式（w=1.0, F=7.5, p=0.05, q=1.0）の軌道．横軸 u，縦軸 v．

以下は，プロットデータを出力するのに用いたプログラム．

#include 

#include 

#include 
void func(double* i_x_, double* i_nx_);

void mul(int i_dim, double* i_x_, double i_scale, double* i_nx_);

void add(int i_dim, double* i_x1_, double* i_x2_, double* i_nx_);

void println(int i_dim, double* i_x_);
//peridic force

double g_w = 1.0;

double g_F = 7.5;
//

double g_p = 0.05;

double g_q = 1.0;


int main(void){

  int i = 0;

  double t_max = 30;

  double t = 0.0;

  double dt = 0.01;

  int dim = 4;

  int N = (int)(t_max/dt);

  double* x_ = new double[dim];

  double* dx_ = new double[dim];
  //oscillator

  x_[0] = 0.0;  

  x_[1] = 1.0;

  //

  x_[2] = 1.0;

  x_[3] = 1.0;
  for(i = N; i >= 0; i--){

    std::cout << t << " ";

    println(dim, x_);

    func(x_, dx_); //dx_ has dx/dt 

    mul(dim, dx_, dt, dx_); //dx = (dx/dt)*dt

    add(dim, x_, dx_, x_); //x = x + dx; 

    t = t + dt;

  }

  delete[] x_;

  delete[] dx_;

}


inline void func(double* i_x_, double* i_nx_){

  double x = i_x_[0];

  double y = i_x_[1];

  double r2 = x*x+y*y; 

  i_nx_[0] = x + g_w*y - x*r2;

  i_nx_[1] = -g_w*x + y - y*r2;
  

  double u = i_x_[2];

  double v = i_x_[3];

  //u' = v

  //v' = -p*v - q*u^3 + F cos(wt)

  i_nx_[2] = v;

  i_nx_[3] = -g_p*v - g_q*u*u*u + g_F * y;

}

inline void mul(int i_dim, double* i_x_, double i_scale, double* i_nx_){

  int i;

  for(i = i_dim - 1; i >= 0; i--){

    i_nx_[i] = i_x_[i] * i_scale;

  }

}

inline void add(int i_dim, double* i_x1_, double* i_x2_, double* i_nx_){

  int i;

  for(i = i_dim - 1; i >= 0; i--){

    i_nx_[i] = i_x1_[i] + i_x2_[i];

  }

}

inline void println(int i_dim, double* i_x_){

  int i;

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

    std::cout << i_x_[i] << " ";

  }

  std::cout << "\n";

}