HOME > natural science Laboratory > コンピュータ・シミュレーション講座 > 仮想物理実験室

2重振子のC++プログラムソース(直交座標系+ルンゲクッタ)

文責:遠藤 理平 (2010年12月24日) カテゴリ:仮想物理実験室(247)

2重振り子などの拘束条件を伴う数値計算を行う場合、変数変換して極座標形式がよく利用されます。 しかしながら、直交座標形式でラグランジュ未定定数法を用いたほうがシミュレーションを行う上で見通しが良い場合が多いです。 ラグランジュ未定乗数法を用いた2重振子のシミュレーションのプログラムソースを公開します。

2重振子のシミュレーション

2つの2重振子を用意し、 L_{12} の長さを 1% 変えたてシミュレーションします。 詳細はラグランジュ未定乗数法を用いた2重振子のシミュレーションをご覧ください。
シルバー L_{12} = 10.0
ゴールド L_{12} = 10.1

計算結果

中心座標 (0,0,40) から長さ 15 のひもで固定された2つの球の時間発展を計算します。 以下の図は、ステップごとの座標をファイルに書き出しプロットした結果です。
赤は球1、緑は球2の軌跡です。
初期位置:球1 r1=(0,0,55)、球2 r2=(0,0,70)
初期速度:球1 v1=(0,0,0)、 球2 v2=(1,0,0)

球の動きは非常に複雑に見えます。初期値を少しずらしただけで軌道が全く変わってしまう運動、カオス運動となることが知られています。 初期値を少し変えて、計算してみてください。

C++プログラムソース 2重振子(直交座標系+ルンゲクッタ)

C++と言っても、C++的なことは使ってませんが。

//////////////////////////////////////////////////////////////////////////
// 2重振子(直交座標系+ルンゲクッタ)
//////////////////////////////////////////////////////////////////////////
#include <math.h>
#include <fstream>
#include <sstream>
#include <iostream>
using namespace std;
double PI = acos(-1.0);

//////////////////////////////////////////////////////////////////////////
// 変数の定義
//////////////////////////////////////////////////////////////////////////
ofstream file1( "2-furiko.data" );

//--------------------------------------------------------
// 仮想物理実験室変数の定義
//--------------------------------------------------------
double t = 0.0;       //時刻
double t_end = 30.0;//計算を終了する時刻
double dt= 0.001;     //時間刻み
int skipDat = 1;   //データ書出し回数の間引き回数

double l = 15.0;    //半径
double g =9.8;
double radian, radian_phi, radian_theta; 
double theta;
//支柱の定義
double poll_x = 0.0;
double poll_y = 0.0;
double poll_z = 40.0;


// ボールの定義
double ball_r = 3.0;//ボールの半径
struct BALL {
  double x0, y0, z0;
  double x, y, z;
  double vx, vy, vz;
  double l;
  double m;
};
const int N = 4;
BALL ball[N];
void SetUp(void){
  for(int i=0; i<N; i++ ){
    ball[i].l  = l;
    ball[i].m  = 1.0;
    ball[i].x  = 0.0;
    ball[i].y  = 0.0;
    ball[i].z  = 0.0;
    ball[i].vx = 0.0;
    ball[i].vy = 0.0;
    ball[i].vz = 0.0;
  }
  ball[0].z  = poll_z + l;
  ball[1].z  = poll_z + l * 2.0;
  ball[2].z  = poll_z + l;
  ball[3].z  = poll_z + l + l*1.01;

  ball[1].vx  = 1.0;
  ball[1].vy  = 0.0;
  ball[3].vx  = 1.0;
  ball[3].vy  = 0.0;
}
//////////////////////////////////////////
// ルンゲクッタ用
void RungeKutta4(BALL ball1[], BALL ball2[]);
void Calculate();

int main(){
  SetUp();
  for(int tn=0; tn <= int(t_end/dt) ;tn++){
    t =dt* tn;
    if(tn%skipDat == 0){
      // 2つの2重振子の座標を書き出す
      file1 << t << " " <<  ball[0].x << " " << ball[0].z << " " << ball[1].x << " " << ball[1].z << " " << ball[2].x << " " << ball[2].z << " " << ball[3].x << " " << ball[3].z << endl;
    }
    Calculate();
  }
  return 0;
}

//--------------------------------------------------------
// 計算と物体の描画
//--------------------------------------------------------
void Calculate(){
    //位置の算出
    RungeKutta4( &ball[0], &ball[1]);
    RungeKutta4( &ball[2], &ball[3]);
}
/////////////////////////////////////////////
// ルンゲクッタ
double fx1(double t, double m1, double x1, double y1, double z1, double vx1, double vy1, double vz1, double m2, double x2, double y2, double z2, double vx2, double vy2, double vz2){
  return vx1;
}
double fy1(double t, double m1, double x1, double y1, double z1, double vx1, double vy1, double vz1, double m2, double x2, double y2, double z2, double vx2, double vy2, double vz2){
  return vy1;
}
double fz1(double t, double m1, double x1, double y1, double z1, double vx1, double vy1, double vz1, double m2, double x2, double y2, double z2, double vx2, double vy2, double vz2){
  return vz1;
}
double fvx1(double t, double m1, double x1, double y1, double z1, double vx1, double vy1, double vz1, double m2, double x2, double y2, double z2, double vx2, double vy2, double vz2){
  double L01 = sqrt(pow(x1-poll_x,2)+pow(y1-poll_y,2)+pow(z1-poll_z,2));
  double L12 = sqrt(pow(x2-x1,2)+pow(y2-y1,2)+pow(z2-z1,2));
  double cos_12 = ((x1-poll_x)*(x2-x1)+(y1-poll_y)*(y2-y1)+(z1-poll_z)*(z2-z1))/(L01*L12);
  double lamda1 = (L12 * (m1+m2) * ( (pow(vx1,2) + pow(vy1,2) + pow(vz1,2)) - g * (z1 -poll_z) ) + m2 * L01 * cos_12 * (pow(vx2-vx1,2)+pow(vy2-vy1,2)+pow(vz2-vz1,2)) ) / (pow(L01,2)*L12*(m1+m2*(1.0-pow(cos_12,2))));
  double lamda2 = (L12 * cos_12  * ( (pow(vx1,2) + pow(vy1,2) + pow(vz1,2)) - g * (z1 -poll_z) )      + L01          * (pow(vx2-vx1,2)+pow(vy2-vy1,2)+pow(vz2-vz1,2)) ) / (L01*pow(L12,2)*(m1+m2*(1.0-pow(cos_12,2))));
  return - lamda1*(x1-poll_x) -m2*lamda2*(x1-x2);
}
double fvy1(double t, double m1, double x1, double y1, double z1, double vx1, double vy1, double vz1, double m2, double x2, double y2, double z2, double vx2, double vy2, double vz2){
  double L01 = sqrt(pow(x1-poll_x,2)+pow(y1-poll_y,2)+pow(z1-poll_z,2));
  double L12 = sqrt(pow(x2-x1,2)+pow(y2-y1,2)+pow(z2-z1,2));
  double cos_12 = ((x1-poll_x)*(x2-x1)+(y1-poll_y)*(y2-y1)+(z1-poll_z)*(z2-z1))/(L01*L12);
  double lamda1 = (L12 * (m1+m2) * ( (pow(vx1,2) + pow(vy1,2) + pow(vz1,2)) - g * (z1 -poll_z) ) + m2 * L01 * cos_12 * (pow(vx2-vx1,2)+pow(vy2-vy1,2)+pow(vz2-vz1,2)) ) / (pow(L01,2)*L12*(m1+m2*(1.0-pow(cos_12,2))));
  double lamda2 = (L12 * cos_12  * ( (pow(vx1,2) + pow(vy1,2) + pow(vz1,2)) - g * (z1 -poll_z) )      + L01          * (pow(vx2-vx1,2)+pow(vy2-vy1,2)+pow(vz2-vz1,2)) ) / (L01*pow(L12,2)*(m1+m2*(1.0-pow(cos_12,2))));
  return - lamda1*(y1-poll_y) -m2*lamda2*(y1-y2);
}
double fvz1(double t, double m1, double x1, double y1, double z1, double vx1, double vy1, double vz1, double m2, double x2, double y2, double z2, double vx2, double vy2, double vz2){
  double L01 = sqrt(pow(x1-poll_x,2)+pow(y1-poll_y,2)+pow(z1-poll_z,2));
  double L12 = sqrt(pow(x2-x1,2)+pow(y2-y1,2)+pow(z2-z1,2));
  double cos_12 = ((x1-poll_x)*(x2-x1)+(y1-poll_y)*(y2-y1)+(z1-poll_z)*(z2-z1))/(L01*L12);
  double lamda1 = (L12 * (m1+m2) * ( (pow(vx1,2) + pow(vy1,2) + pow(vz1,2)) - g * (z1 -poll_z) ) + m2 * L01 * cos_12 * (pow(vx2-vx1,2)+pow(vy2-vy1,2)+pow(vz2-vz1,2)) ) / (pow(L01,2)*L12*(m1+m2*(1.0-pow(cos_12,2))));
  double lamda2 = (L12 * cos_12  * ( (pow(vx1,2) + pow(vy1,2) + pow(vz1,2)) - g * (z1 -poll_z) )      + L01          * (pow(vx2-vx1,2)+pow(vy2-vy1,2)+pow(vz2-vz1,2)) ) / (L01*pow(L12,2)*(m1+m2*(1.0-pow(cos_12,2))));
  return - g - lamda1*(z1-poll_z) - m2*lamda2*(z1-z2);
}
double fx2(double t, double m1, double x1, double y1, double z1, double vx1, double vy1, double vz1, double m2, double x2, double y2, double z2, double vx2, double vy2, double vz2){
  return vx2;
}
double fy2(double t, double m1, double x1, double y1, double z1, double vx1, double vy1, double vz1, double m2, double x2, double y2, double z2, double vx2, double vy2, double vz2){
  return vy2;
}
double fz2(double t, double m1, double x1, double y1, double z1, double vx1, double vy1, double vz1, double m2, double x2, double y2, double z2, double vx2, double vy2, double vz2){
  return vz2;
}
double fvx2(double t, double m1, double x1, double y1, double z1, double vx1, double vy1, double vz1, double m2, double x2, double y2, double z2, double vx2, double vy2, double vz2){
  double L01 = sqrt(pow(x1-poll_x,2)+pow(y1-poll_y,2)+pow(z1-poll_z,2));
  double L12 = sqrt(pow(x2-x1,2)+pow(y2-y1,2)+pow(z2-z1,2));
  double cos_12 = ((x1-poll_x)*(x2-x1)+(y1-poll_y)*(y2-y1)+(z1-poll_z)*(z2-z1))/(L01*L12);
  double lamda2 = (L12 * cos_12  * ( (pow(vx1,2) + pow(vy1,2) + pow(vz1,2)) - g * (z1 -poll_z) )      + L01          * (pow(vx2-vx1,2)+pow(vy2-vy1,2)+pow(vz2-vz1,2)) ) / (L01*pow(L12,2)*(m1+m2*(1.0-pow(cos_12,2))));
  return  -m1 * lamda2*(x2-x1);
}
double fvy2(double t, double m1, double x1, double y1, double z1, double vx1, double vy1, double vz1, double m2, double x2, double y2, double z2, double vx2, double vy2, double vz2){
  double L01 = sqrt(pow(x1-poll_x,2)+pow(y1-poll_y,2)+pow(z1-poll_z,2));
  double L12 = sqrt(pow(x2-x1,2)+pow(y2-y1,2)+pow(z2-z1,2));
  double cos_12 = ((x1-poll_x)*(x2-x1)+(y1-poll_y)*(y2-y1)+(z1-poll_z)*(z2-z1))/(L01*L12);
  double lamda2 = (L12 * cos_12  * ( (pow(vx1,2) + pow(vy1,2) + pow(vz1,2)) - g * (z1 -poll_z) )      + L01          * (pow(vx2-vx1,2)+pow(vy2-vy1,2)+pow(vz2-vz1,2)) ) / (L01*pow(L12,2)*(m1+m2*(1.0-pow(cos_12,2))));
  return  -m1 * lamda2*(y2-y1);
}
double fvz2(double t, double m1, double x1, double y1, double z1, double vx1, double vy1, double vz1, double m2, double x2, double y2, double z2, double vx2, double vy2, double vz2){
  double L01 = sqrt(pow(x1-poll_x,2)+pow(y1-poll_y,2)+pow(z1-poll_z,2));
  double L12 = sqrt(pow(x2-x1,2)+pow(y2-y1,2)+pow(z2-z1,2));
  double cos_12 = ((x1-poll_x)*(x2-x1)+(y1-poll_y)*(y2-y1)+(z1-poll_z)*(z2-z1))/(L01*L12);
  double lamda2 = (L12 * cos_12  * ( (pow(vx1,2) + pow(vy1,2) + pow(vz1,2)) - g * (z1 -poll_z) )      + L01          * (pow(vx2-vx1,2)+pow(vy2-vy1,2)+pow(vz2-vz1,2)) ) / (L01*pow(L12,2)*(m1+m2*(1.0-pow(cos_12,2))));
  return  -g - m1 * lamda2*(z2-z1);
}
void RungeKutta4(BALL *ball1, BALL *ball2){
  double k1[2][6],k2[2][6],k3[2][6],k4[2][6];
  double x1  = ball1 -> x;
  double vx1 = ball1 -> vx;
  double y1  = ball1 -> y;
  double vy1 = ball1 -> vy;
  double z1  = ball1 -> z;
  double vz1 = ball1 -> vz;
  double m1  = ball1 -> m;
  double x2  = ball2 -> x;
  double vx2 = ball2 -> vx;
  double y2  = ball2 -> y;
  double vy2 = ball2 -> vy;
  double z2  = ball2 -> z;
  double vz2 = ball2 -> vz;
  double m2  = ball2 -> m;
  //cout << sqrt(pow(x1-poll_x,2)+pow(y1-poll_y,2)+pow(z1-poll_z,2))<< " " <<  sqrt(pow(x2-x1,2)+pow(y2-y1,2)+pow(z2-z1,2)) << endl;

  k1[0][0]=dt*fx1( t,m1,x1,y1,z1,vx1,vy1,vz1,m2,x2,y2,z2,vx2,vy2,vz2);
  k1[0][1]=dt*fvx1(t,m1,x1,y1,z1,vx1,vy1,vz1,m2,x2,y2,z2,vx2,vy2,vz2);
  k1[0][2]=dt*fy1( t,m1,x1,y1,z1,vx1,vy1,vz1,m2,x2,y2,z2,vx2,vy2,vz2);
  k1[0][3]=dt*fvy1(t,m1,x1,y1,z1,vx1,vy1,vz1,m2,x2,y2,z2,vx2,vy2,vz2);
  k1[0][4]=dt*fz1( t,m1,x1,y1,z1,vx1,vy1,vz1,m2,x2,y2,z2,vx2,vy2,vz2);
  k1[0][5]=dt*fvz1(t,m1,x1,y1,z1,vx1,vy1,vz1,m2,x2,y2,z2,vx2,vy2,vz2);
  k1[1][0]=dt*fx2( t,m1,x1,y1,z1,vx1,vy1,vz1,m2,x2,y2,z2,vx2,vy2,vz2);
  k1[1][1]=dt*fvx2(t,m1,x1,y1,z1,vx1,vy1,vz1,m2,x2,y2,z2,vx2,vy2,vz2);
  k1[1][2]=dt*fy2( t,m1,x1,y1,z1,vx1,vy1,vz1,m2,x2,y2,z2,vx2,vy2,vz2);
  k1[1][3]=dt*fvy2(t,m1,x1,y1,z1,vx1,vy1,vz1,m2,x2,y2,z2,vx2,vy2,vz2);
  k1[1][4]=dt*fz2( t,m1,x1,y1,z1,vx1,vy1,vz1,m2,x2,y2,z2,vx2,vy2,vz2);
  k1[1][5]=dt*fvz2(t,m1,x1,y1,z1,vx1,vy1,vz1,m2,x2,y2,z2,vx2,vy2,vz2);

  k2[0][0]=dt*fx1( t+dt/2.0,m1,x1+k1[0][0]/2.0,y1+k1[0][2]/2.0,z1+k1[0][4]/2.0,vx1+k1[0][1]/2.0,vy1+k1[0][3]/2.0,vz1+k1[0][5]/2.0,  m2,x2+k1[1][0]/2.0,y2+k1[1][2]/2.0,z2+k1[1][4]/2.0,vx2+k1[1][1]/2.0,vy2+k1[1][3]/2.0,vz2+k1[1][5]/2.0);
  k2[0][1]=dt*fvx1(t+dt/2.0,m1,x1+k1[0][0]/2.0,y1+k1[0][2]/2.0,z1+k1[0][4]/2.0,vx1+k1[0][1]/2.0,vy1+k1[0][3]/2.0,vz1+k1[0][5]/2.0,  m2,x2+k1[1][0]/2.0,y2+k1[1][2]/2.0,z2+k1[1][4]/2.0,vx2+k1[1][1]/2.0,vy2+k1[1][3]/2.0,vz2+k1[1][5]/2.0);
  k2[0][2]=dt*fy1( t+dt/2.0,m1,x1+k1[0][0]/2.0,y1+k1[0][2]/2.0,z1+k1[0][4]/2.0,vx1+k1[0][1]/2.0,vy1+k1[0][3]/2.0,vz1+k1[0][5]/2.0,  m2,x2+k1[1][0]/2.0,y2+k1[1][2]/2.0,z2+k1[1][4]/2.0,vx2+k1[1][1]/2.0,vy2+k1[1][3]/2.0,vz2+k1[1][5]/2.0);
  k2[0][3]=dt*fvy1(t+dt/2.0,m1,x1+k1[0][0]/2.0,y1+k1[0][2]/2.0,z1+k1[0][4]/2.0,vx1+k1[0][1]/2.0,vy1+k1[0][3]/2.0,vz1+k1[0][5]/2.0,  m2,x2+k1[1][0]/2.0,y2+k1[1][2]/2.0,z2+k1[1][4]/2.0,vx2+k1[1][1]/2.0,vy2+k1[1][3]/2.0,vz2+k1[1][5]/2.0);
  k2[0][4]=dt*fz1( t+dt/2.0,m1,x1+k1[0][0]/2.0,y1+k1[0][2]/2.0,z1+k1[0][4]/2.0,vx1+k1[0][1]/2.0,vy1+k1[0][3]/2.0,vz1+k1[0][5]/2.0,  m2,x2+k1[1][0]/2.0,y2+k1[1][2]/2.0,z2+k1[1][4]/2.0,vx2+k1[1][1]/2.0,vy2+k1[1][3]/2.0,vz2+k1[1][5]/2.0);
  k2[0][5]=dt*fvz1(t+dt/2.0,m1,x1+k1[0][0]/2.0,y1+k1[0][2]/2.0,z1+k1[0][4]/2.0,vx1+k1[0][1]/2.0,vy1+k1[0][3]/2.0,vz1+k1[0][5]/2.0,  m2,x2+k1[1][0]/2.0,y2+k1[1][2]/2.0,z2+k1[1][4]/2.0,vx2+k1[1][1]/2.0,vy2+k1[1][3]/2.0,vz2+k1[1][5]/2.0);
  k2[1][0]=dt*fx2( t+dt/2.0,m1,x1+k1[0][0]/2.0,y1+k1[0][2]/2.0,z1+k1[0][4]/2.0,vx1+k1[0][1]/2.0,vy1+k1[0][3]/2.0,vz1+k1[0][5]/2.0,  m2,x2+k1[1][0]/2.0,y2+k1[1][2]/2.0,z2+k1[1][4]/2.0,vx2+k1[1][1]/2.0,vy2+k1[1][3]/2.0,vz2+k1[1][5]/2.0);
  k2[1][1]=dt*fvx2(t+dt/2.0,m1,x1+k1[0][0]/2.0,y1+k1[0][2]/2.0,z1+k1[0][4]/2.0,vx1+k1[0][1]/2.0,vy1+k1[0][3]/2.0,vz1+k1[0][5]/2.0,  m2,x2+k1[1][0]/2.0,y2+k1[1][2]/2.0,z2+k1[1][4]/2.0,vx2+k1[1][1]/2.0,vy2+k1[1][3]/2.0,vz2+k1[1][5]/2.0);
  k2[1][2]=dt*fy2( t+dt/2.0,m1,x1+k1[0][0]/2.0,y1+k1[0][2]/2.0,z1+k1[0][4]/2.0,vx1+k1[0][1]/2.0,vy1+k1[0][3]/2.0,vz1+k1[0][5]/2.0,  m2,x2+k1[1][0]/2.0,y2+k1[1][2]/2.0,z2+k1[1][4]/2.0,vx2+k1[1][1]/2.0,vy2+k1[1][3]/2.0,vz2+k1[1][5]/2.0);
  k2[1][3]=dt*fvy2(t+dt/2.0,m1,x1+k1[0][0]/2.0,y1+k1[0][2]/2.0,z1+k1[0][4]/2.0,vx1+k1[0][1]/2.0,vy1+k1[0][3]/2.0,vz1+k1[0][5]/2.0,  m2,x2+k1[1][0]/2.0,y2+k1[1][2]/2.0,z2+k1[1][4]/2.0,vx2+k1[1][1]/2.0,vy2+k1[1][3]/2.0,vz2+k1[1][5]/2.0);
  k2[1][4]=dt*fz2( t+dt/2.0,m1,x1+k1[0][0]/2.0,y1+k1[0][2]/2.0,z1+k1[0][4]/2.0,vx1+k1[0][1]/2.0,vy1+k1[0][3]/2.0,vz1+k1[0][5]/2.0,  m2,x2+k1[1][0]/2.0,y2+k1[1][2]/2.0,z2+k1[1][4]/2.0,vx2+k1[1][1]/2.0,vy2+k1[1][3]/2.0,vz2+k1[1][5]/2.0);
  k2[1][5]=dt*fvz2(t+dt/2.0,m1,x1+k1[0][0]/2.0,y1+k1[0][2]/2.0,z1+k1[0][4]/2.0,vx1+k1[0][1]/2.0,vy1+k1[0][3]/2.0,vz1+k1[0][5]/2.0,  m2,x2+k1[1][0]/2.0,y2+k1[1][2]/2.0,z2+k1[1][4]/2.0,vx2+k1[1][1]/2.0,vy2+k1[1][3]/2.0,vz2+k1[1][5]/2.0);

  k3[0][0]=dt*fx1( t+dt/2.0,m1,x1+k2[0][0]/2.0,y1+k2[0][2]/2.0,z1+k2[0][4]/2.0,vx1+k2[0][1]/2.0,vy1+k2[0][3]/2.0,vz1+k2[0][5]/2.0   ,m2,x2+k2[1][0]/2.0,y2+k2[1][2]/2.0,z2+k2[1][4]/2.0,vx2+k2[1][1]/2.0,vy2+k2[1][3]/2.0,vz2+k2[1][5]/2.0);
  k3[0][1]=dt*fvx1(t+dt/2.0,m1,x1+k2[0][0]/2.0,y1+k2[0][2]/2.0,z1+k2[0][4]/2.0,vx1+k2[0][1]/2.0,vy1+k2[0][3]/2.0,vz1+k2[0][5]/2.0   ,m2,x2+k2[1][0]/2.0,y2+k2[1][2]/2.0,z2+k2[1][4]/2.0,vx2+k2[1][1]/2.0,vy2+k2[1][3]/2.0,vz2+k2[1][5]/2.0);
  k3[0][2]=dt*fy1( t+dt/2.0,m1,x1+k2[0][0]/2.0,y1+k2[0][2]/2.0,z1+k2[0][4]/2.0,vx1+k2[0][1]/2.0,vy1+k2[0][3]/2.0,vz1+k2[0][5]/2.0   ,m2,x2+k2[1][0]/2.0,y2+k2[1][2]/2.0,z2+k2[1][4]/2.0,vx2+k2[1][1]/2.0,vy2+k2[1][3]/2.0,vz2+k2[1][5]/2.0);
  k3[0][3]=dt*fvy1(t+dt/2.0,m1,x1+k2[0][0]/2.0,y1+k2[0][2]/2.0,z1+k2[0][4]/2.0,vx1+k2[0][1]/2.0,vy1+k2[0][3]/2.0,vz1+k2[0][5]/2.0   ,m2,x2+k2[1][0]/2.0,y2+k2[1][2]/2.0,z2+k2[1][4]/2.0,vx2+k2[1][1]/2.0,vy2+k2[1][3]/2.0,vz2+k2[1][5]/2.0);
  k3[0][4]=dt*fz1( t+dt/2.0,m1,x1+k2[0][0]/2.0,y1+k2[0][2]/2.0,z1+k2[0][4]/2.0,vx1+k2[0][1]/2.0,vy1+k2[0][3]/2.0,vz1+k2[0][5]/2.0   ,m2,x2+k2[1][0]/2.0,y2+k2[1][2]/2.0,z2+k2[1][4]/2.0,vx2+k2[1][1]/2.0,vy2+k2[1][3]/2.0,vz2+k2[1][5]/2.0);
  k3[0][5]=dt*fvz1(t+dt/2.0,m1,x1+k2[0][0]/2.0,y1+k2[0][2]/2.0,z1+k2[0][4]/2.0,vx1+k2[0][1]/2.0,vy1+k2[0][3]/2.0,vz1+k2[0][5]/2.0   ,m2,x2+k2[1][0]/2.0,y2+k2[1][2]/2.0,z2+k2[1][4]/2.0,vx2+k2[1][1]/2.0,vy2+k2[1][3]/2.0,vz2+k2[1][5]/2.0);
  k3[1][0]=dt*fx2( t+dt/2.0,m1,x1+k2[0][0]/2.0,y1+k2[0][2]/2.0,z1+k2[0][4]/2.0,vx1+k2[0][1]/2.0,vy1+k2[0][3]/2.0,vz1+k2[0][5]/2.0   ,m2,x2+k2[1][0]/2.0,y2+k2[1][2]/2.0,z2+k2[1][4]/2.0,vx2+k2[1][1]/2.0,vy2+k2[1][3]/2.0,vz2+k2[1][5]/2.0);
  k3[1][1]=dt*fvx2(t+dt/2.0,m1,x1+k2[0][0]/2.0,y1+k2[0][2]/2.0,z1+k2[0][4]/2.0,vx1+k2[0][1]/2.0,vy1+k2[0][3]/2.0,vz1+k2[0][5]/2.0   ,m2,x2+k2[1][0]/2.0,y2+k2[1][2]/2.0,z2+k2[1][4]/2.0,vx2+k2[1][1]/2.0,vy2+k2[1][3]/2.0,vz2+k2[1][5]/2.0);
  k3[1][2]=dt*fy2( t+dt/2.0,m1,x1+k2[0][0]/2.0,y1+k2[0][2]/2.0,z1+k2[0][4]/2.0,vx1+k2[0][1]/2.0,vy1+k2[0][3]/2.0,vz1+k2[0][5]/2.0   ,m2,x2+k2[1][0]/2.0,y2+k2[1][2]/2.0,z2+k2[1][4]/2.0,vx2+k2[1][1]/2.0,vy2+k2[1][3]/2.0,vz2+k2[1][5]/2.0);
  k3[1][3]=dt*fvy2(t+dt/2.0,m1,x1+k2[0][0]/2.0,y1+k2[0][2]/2.0,z1+k2[0][4]/2.0,vx1+k2[0][1]/2.0,vy1+k2[0][3]/2.0,vz1+k2[0][5]/2.0   ,m2,x2+k2[1][0]/2.0,y2+k2[1][2]/2.0,z2+k2[1][4]/2.0,vx2+k2[1][1]/2.0,vy2+k2[1][3]/2.0,vz2+k2[1][5]/2.0);
  k3[1][4]=dt*fz2( t+dt/2.0,m1,x1+k2[0][0]/2.0,y1+k2[0][2]/2.0,z1+k2[0][4]/2.0,vx1+k2[0][1]/2.0,vy1+k2[0][3]/2.0,vz1+k2[0][5]/2.0   ,m2,x2+k2[1][0]/2.0,y2+k2[1][2]/2.0,z2+k2[1][4]/2.0,vx2+k2[1][1]/2.0,vy2+k2[1][3]/2.0,vz2+k2[1][5]/2.0);
  k3[1][5]=dt*fvz2(t+dt/2.0,m1,x1+k2[0][0]/2.0,y1+k2[0][2]/2.0,z1+k2[0][4]/2.0,vx1+k2[0][1]/2.0,vy1+k2[0][3]/2.0,vz1+k2[0][5]/2.0   ,m2,x2+k2[1][0]/2.0,y2+k2[1][2]/2.0,z2+k2[1][4]/2.0,vx2+k2[1][1]/2.0,vy2+k2[1][3]/2.0,vz2+k2[1][5]/2.0);

  k4[0][0]=dt*fx1( t+dt/2.0,m1,x1+k3[0][0],y1+k3[0][2],z1+k3[0][4],vx1+k3[0][1],vy1+k3[0][3],vz1+k3[0][5]  ,m2,x2+k3[1][0],y2+k3[1][2],z2+k3[1][4],vx2+k3[1][1],vy2+k3[1][3],vz2+k3[1][5]);
  k4[0][1]=dt*fvx1(t+dt/2.0,m1,x1+k3[0][0],y1+k3[0][2],z1+k3[0][4],vx1+k3[0][1],vy1+k3[0][3],vz1+k3[0][5]  ,m2,x2+k3[1][0],y2+k3[1][2],z2+k3[1][4],vx2+k3[1][1],vy2+k3[1][3],vz2+k3[1][5]);
  k4[0][2]=dt*fy1( t+dt/2.0,m1,x1+k3[0][0],y1+k3[0][2],z1+k3[0][4],vx1+k3[0][1],vy1+k3[0][3],vz1+k3[0][5]  ,m2,x2+k3[1][0],y2+k3[1][2],z2+k3[1][4],vx2+k3[1][1],vy2+k3[1][3],vz2+k3[1][5]);
  k4[0][3]=dt*fvy1(t+dt/2.0,m1,x1+k3[0][0],y1+k3[0][2],z1+k3[0][4],vx1+k3[0][1],vy1+k3[0][3],vz1+k3[0][5]  ,m2,x2+k3[1][0],y2+k3[1][2],z2+k3[1][4],vx2+k3[1][1],vy2+k3[1][3],vz2+k3[1][5]);
  k4[0][4]=dt*fz1( t+dt/2.0,m1,x1+k3[0][0],y1+k3[0][2],z1+k3[0][4],vx1+k3[0][1],vy1+k3[0][3],vz1+k3[0][5]  ,m2,x2+k3[1][0],y2+k3[1][2],z2+k3[1][4],vx2+k3[1][1],vy2+k3[1][3],vz2+k3[1][5]);
  k4[0][5]=dt*fvz1(t+dt/2.0,m1,x1+k3[0][0],y1+k3[0][2],z1+k3[0][4],vx1+k3[0][1],vy1+k3[0][3],vz1+k3[0][5]  ,m2,x2+k3[1][0],y2+k3[1][2],z2+k3[1][4],vx2+k3[1][1],vy2+k3[1][3],vz2+k3[1][5]);
  k4[1][0]=dt*fx2( t+dt/2.0,m1,x1+k3[0][0],y1+k3[0][2],z1+k3[0][4],vx1+k3[0][1],vy1+k3[0][3],vz1+k3[0][5]  ,m2,x2+k3[1][0],y2+k3[1][2],z2+k3[1][4],vx2+k3[1][1],vy2+k3[1][3],vz2+k3[1][5]);
  k4[1][1]=dt*fvx2(t+dt/2.0,m1,x1+k3[0][0],y1+k3[0][2],z1+k3[0][4],vx1+k3[0][1],vy1+k3[0][3],vz1+k3[0][5]  ,m2,x2+k3[1][0],y2+k3[1][2],z2+k3[1][4],vx2+k3[1][1],vy2+k3[1][3],vz2+k3[1][5]);
  k4[1][2]=dt*fy2( t+dt/2.0,m1,x1+k3[0][0],y1+k3[0][2],z1+k3[0][4],vx1+k3[0][1],vy1+k3[0][3],vz1+k3[0][5]  ,m2,x2+k3[1][0],y2+k3[1][2],z2+k3[1][4],vx2+k3[1][1],vy2+k3[1][3],vz2+k3[1][5]);
  k4[1][3]=dt*fvy2(t+dt/2.0,m1,x1+k3[0][0],y1+k3[0][2],z1+k3[0][4],vx1+k3[0][1],vy1+k3[0][3],vz1+k3[0][5]  ,m2,x2+k3[1][0],y2+k3[1][2],z2+k3[1][4],vx2+k3[1][1],vy2+k3[1][3],vz2+k3[1][5]);
  k4[1][4]=dt*fz2( t+dt/2.0,m1,x1+k3[0][0],y1+k3[0][2],z1+k3[0][4],vx1+k3[0][1],vy1+k3[0][3],vz1+k3[0][5]  ,m2,x2+k3[1][0],y2+k3[1][2],z2+k3[1][4],vx2+k3[1][1],vy2+k3[1][3],vz2+k3[1][5]);
  k4[1][5]=dt*fvz2(t+dt/2.0,m1,x1+k3[0][0],y1+k3[0][2],z1+k3[0][4],vx1+k3[0][1],vy1+k3[0][3],vz1+k3[0][5]  ,m2,x2+k3[1][0],y2+k3[1][2],z2+k3[1][4],vx2+k3[1][1],vy2+k3[1][3],vz2+k3[1][5]);

  ball1->x  = ball1->x  + (k1[0][0]+2.0*k2[0][0]+2.0*k3[0][0]+k4[0][0])/6.0;
  ball1->vx = ball1->vx + (k1[0][1]+2.0*k2[0][1]+2.0*k3[0][1]+k4[0][1])/6.0;
  ball1->y  = ball1->y  + (k1[0][2]+2.0*k2[0][2]+2.0*k3[0][2]+k4[0][2])/6.0;
  ball1->vy = ball1->vy + (k1[0][3]+2.0*k2[0][3]+2.0*k3[0][3]+k4[0][3])/6.0;
  ball1->z  = ball1->z  + (k1[0][4]+2.0*k2[0][4]+2.0*k3[0][4]+k4[0][4])/6.0;
  ball1->vz = ball1->vz + (k1[0][5]+2.0*k2[0][5]+2.0*k3[0][5]+k4[0][5])/6.0;
  ball2->x  = ball2->x  + (k1[1][0]+2.0*k2[1][0]+2.0*k3[1][0]+k4[1][0])/6.0;
  ball2->vx = ball2->vx + (k1[1][1]+2.0*k2[1][1]+2.0*k3[1][1]+k4[1][1])/6.0;
  ball2->y  = ball2->y  + (k1[1][2]+2.0*k2[1][2]+2.0*k3[1][2]+k4[1][2])/6.0;
  ball2->vy = ball2->vy + (k1[1][3]+2.0*k2[1][3]+2.0*k3[1][3]+k4[1][3])/6.0;
  ball2->z  = ball2->z  + (k1[1][4]+2.0*k2[1][4]+2.0*k3[1][4]+k4[1][4])/6.0;
  ball2->vz = ball2->vz + (k1[1][5]+2.0*k2[1][5]+2.0*k3[1][5]+k4[1][5])/6.0;
}

gnuplotのテンプレート

ついでに数値計算の結果を gnuplot でプロットするテンプレートも記載しておきます。

set terminal gif optimize size 600, 600
set tics font 'Times,18'
set size ratio 1.0
set nokey

set output "2-furiko.gif"
plot "2-furiko.data" u 2:3 w d ,  "2-furiko.data" u 4:5 w d


▲このページのトップNPO法人 natural science トップ

関連記事

仮想物理実験室







▲このページのトップNPO法人 natural science トップ