プログラム

Macwiki様の3次元ベクトルクラスを使わせていただいています。
http://macwiki.sourceforge.jp/wiki/index.php/Vector3%E3%82%AF%E3%83%A9%E3%82%B9

using namespace std;
#include<iostream>
#include<math.h>
#include<stdio.h>
#include"vector3.h"
#define     PI      3.141592653589793238463
#define		G		 9.80665

vector3 sig(vector3 (&R1)[3],vector3 (&R2)[3], double l){
	vector3 temp0,temp1,temp3,temp4;
	temp1=R1[1]-R2[1];
	temp0=R1[0]-R2[0];
	temp3=-k*(1.0/l-1.0/temp1.abs())*temp1;
	temp4=-g*(1-temp0.abs()/temp1.abs())/dT*temp1;
	return(temp3+temp4);
}

void main(void){
	vector3 r1[3],r2[3],r3[3];
	vector3 T12=vector3(0.0,0.0,5.0);//N/m
	int i;
	for (i=0;i<3;i++){
		r1[i]=vector3(-4.0, 0.0, 0.0);
		r2[i]=vector3( 0.0, 0.0, 0.0);
		r3[i]=vector3( 0.0, 4.0, 0.0);
	}

	double l12=(r1[1]-r2[1]).abs();
	double l23=(r2[1]-r3[1]).abs();

	long int it;
	for(it=0;it*dT<=30;it++){
		r1[2]=2.0*r1[1]-r1[0] + (sig(r1,r2,l12)+T12%(r1[1]-r2[1])) *dT*dT/m1;
		r2[2]=2.0*r2[1]-r2[0] + (sig(r2,r1,l12)+sig(r2,r3,l23)
								+T12%(r2[1]-r1[1])-T12%(r2[1]-r3[1])) *dT*dT/m1;
		r3[2]=2.0*r3[1]-r3[0] + (sig(r3,r2,l23)-T12%(r3[1]-r2[1])) *dT*dT/m3;

		for(i=0;i<2;i++){
			r1[i]=r1[i+1];
			r2[i]=r2[i+1];
			r3[i]=r3[i+1];
		}
		if (it%100==0){
			cout <<r1[1].x<<" "<<r1[1].y<<"\n";
			cout <<r2[1].x<<" "<<r2[1].y<<"\n";
			cout <<"\n"<<"\n";
			cout <<r2[1].x<<" "<<r2[1].y<<"\n";
			cout <<r3[1].x<<" "<<r3[1].y<<"\n";
			cout <<"\n"<<"\n";
		}

	}
}