제한된 3체 문제(restricted three-body problem)를 풀고 있습니다.
한 질량은 좌표의 원점에 정지해있고, 다른 한 질량은 원점을 중심으로 원운동을 하구요.
궤적을 살펴보고자 하는 질량의 초기 위치와 초기 속도를 대입하게 되면 궤적을 구할 수 있게 프로그램을 짰는데요...

원래 제한된 3체 문제에서는 타겟이 되는 질량의 역학적 에너지가 보존이 안되는 건지 궁금하네요..

.. 고민하다가 소스 코드를 올려드립니다.
초기 값을 150 150 0 0
그리고 140 160 0 0 을 넣어보세요..;; 파일은 엑셀로 출력되구요.

A열부터 차례로 시간, x, y, x방향 속도, y방향 속도, 원운동하는 두 번째 질량의 x, y, 타겟의 역학적 에너지의 합입니다..

#include <stdio.h>
#include <math.h>
#define c 20.0
#define m1 40.0
#define m2 20.0
#define R 100.0
#define h 0.01

double fwx(double t, double x, double y)
{
double res;
double wr=sqrt(c*m1/pow(R,3));
res=c*m1/(pow(x,2)+pow(y,2))*(-x)/sqrt(pow(x,2)+pow(y,2))+c*m2/(pow(x-R*cos(wr*t),2)+pow(y-R*sin(wr*t),2))*(-x+R*cos(wr*t))/sqrt(pow(x-R*cos(wr*t),2)+pow(y-R*sin(wr*t),2));
return res;
}

double fwy(double t, double x, double y)
{
double res;
double wr=sqrt(c*m1/pow(R,3));
res=c*m1/(pow(x,2)+pow(y,2))*(-y)/sqrt(pow(x,2)+pow(y,2))+c*m2/(pow(x-R*cos(wr*t),2)+pow(y-R*sin(wr*t),2))*(-y+R*sin(wr*t))/sqrt(pow(x-R*cos(wr*t),2)+pow(y-R*sin(wr*t),2));
return res;
}

double fenergy(double t, double x, double y, double wx, double wy)
{
double res;
double wr=sqrt(c*m1/pow(R,3));
res=-c*m1/sqrt(pow(x,2)+pow(y,2))-c*m2/sqrt(pow(x-R*cos(wr*t),2)+pow(y-R*sin(wr*t),2))+(pow(wx,2)+pow(wy,2))/2.0;
return res;
}

int fimp1(double x, double y)
{
int res;
if(sqrt(pow(x,2)+pow(y,2))<=pow(10,1))
res=0;
else
res=1;
return res;
}

int fimp2(double t, double x, double y)
{
int res;
double wr=sqrt(c*m1/pow(R,3));
if(sqrt(pow(x-R*cos(wr*t),2)+pow(y-R*sin(wr*t),2))<=pow(10,1))
res=0;
else
res=1;
return res;
}

int fimp(int fimp1, int fimp2)
{
int res;
res=fimp1*fimp2;
return res;
}

int main(void)
{
double t, x, y, wx, wy;
double kx[4], kwx[4], ky[4], kwy[4];
double wr=sqrt(c*m1/pow(R,3));
char fname[]= "test 001.xls";
FILE *f;

f=fopen(fname, "w");

printf("input the initial values of x, y, dx/dt, dy/dt in order\n");
scanf("%lf %lf %lf %lf", &x, &y, &wx, &wy);

fprintf(f, "c=%.2lf m1=%.2lf m2=%.2lf R=%.2lf h=%.2lf\n", c, m1, m2, R, h);
fprintf(f, "x0=%.2lf y0=%.2lf vx0=%.2lf vy0=%.2lf\n\n", x, y, wx, wy);

for(t=0;fimp(fimp1(x,y),fimp2(t,x,y));t=t+h)
{
fprintf(f,"%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\n",t,x,y,wx,wy,R*cos(wr*t),R*sin(wr*t),fenergy(t,x,y,wx,wy));

kx[0]=wx;
ky[0]=wy;
kwx[0]=fwx(t,x,y);
kwy[0]=fwy(t,x,y);
kx[1]=wx+0.5*kwx[0]*h;
ky[1]=wy+0.5*kwy[0]*h;
kwx[1]=fwx(t+0.5*h,x+0.5*kx[0]*h,y+0.5*ky[0]*h);
kwy[1]=fwy(t+0.5*h,x+0.5*kx[0]*h,y+0.5*ky[0]*h);
kx[2]=wx+0.5*kwx[1]*h;
ky[2]=wy+0.5*kwy[1]*h;
kwx[2]=fwx(t+0.5*h,x+0.5*kx[1]*h,y+0.5*ky[1]*h);
kwy[2]=fwy(t+0.5*h,x+0.5*kx[1]*h,y+0.5*ky[1]*h);
kx[3]=wx+kwx[2]*h;
ky[3]=wy+kwy[2]*h;
kwx[3]=fwx(t+h,x+kx[2]*h,y+ky[2]*h);
kwy[3]=fwy(t+h,x+kx[2]*h,y+ky[2]*h);

x=x+(kx[0]+2*kx[1]+2*kx[2]+kx[3])*h/6.0;
y=y+(ky[0]+2*ky[1]+2*ky[2]+ky[3])*h/6.0;
wx=wx+(kwx[0]+2*kwx[1]+2*kwx[2]+kwx[3])*h/6.0;
wy=wy+(kwy[0]+2*kwy[1]+2*kwy[2]+kwy[3])*h/6.0;
}

fclose(f);
}


초기값 140 160 0 0 에서.. 이놈이 중심으로부터 되게 멀리 떨어지는 걸 확인할 수 있는데요.
이게 정상인건지 아니면 수치해석하다가 나타난 오류인건지 그것도 잘 모르겠네요;; 역학적 에너지가 보존되는지부터가 궁금함..