VisualC++ と OpenGL を利用した仮想物理実験室
ラグランジュ未定乗数法を用いた２重振子のシミュレーション
文責：遠藤 理平 （2010年11月23日）
カテゴリ：仮想物理実験室(325)
ここまで、ラグランジアンを用いて以下の運動のシミュレーションを行いました。
・ラグランジュ運動方程式１：極座標を用いた単振子
・ラグランジュ運動方程式２：極座標を用いた球面振子
・ラグランジュ未定乗数法を用いた球面振子のシミュレーション
今回は、２重振子の運動をラグランジュ未定乗数法を用いてシミュレーションを行ないます。
２重振子の定義
重力場中を運動する２重振子を上の図のように定義します。
重さのない、伸びないひもで結ばれた２つの球(1,2)を想定します。
ひもは伸びないので、２点間の距離が絶えず一定となるように運動を行ないます。
つまり、２つの球の座標 r_1, r_2 は独立に動いているのではなくて、２点間の距離が一定という条件を必ず満たしています。
このように運動を拘束する条件は拘束条件と呼ばれます。
拘束条件が存在する運動を記述する際には、ラグランジュ未定乗数法と呼ばれる解析力学の手法を用いることが
非常に有用で有ることが知られています。
運動方程式
λ_1 とλ_2 は未定乗数と呼ばれる量で、この場合、ひもが球を引く張力を表します。 λ_1 とλ_2 の具体的な値は、拘束条件から導くことができます。
ラグランジュの未定定数
拘束条件と、その１階微分、２階微分から連立方程式を解くことで、 ラグランジュの未定定数を決定することができます。
２重振子のシミュレーション
２つの２重振子を用意し、
L_{12} の長さを 1% 変えたてシミュレーションします。
シルバー L_{12} = 10.0
ゴールド L_{12} = 10.1
極座標を用いた単振子＋ルンゲクッタ法のプログラム
【0-2-3】仮想物理実験室の構築 （ver1.2）を参照ください
構造体の定義
//支柱の定義
double poll_x = 0.0;
double poll_y = 0.0;
double poll_z = 40.0;
double poll_r = 2.0; //半径
double poll_h = 6.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 = poll_z + l * double(i+1);
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;
}
ルンゲクッタ法を利用した数値計算
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;
}
