ピタゴラス3体問題2

Δt 時間後の m3、m4 の位置、速度を求める函数
m5 の位置は、微分方程式を解かなくても m3、m4 の位置から一意に求まるので別に計算する。
f:id:ti-nspire:20160509112726j:plain

Define pyt(x30,vx30,y30,vy30,x40,vx40,y40,vy40,dt,tol)=
Func
:Local r34_3,r35_3,r45_3
:Define r34_3=√((x3-x4)^(2)+(y3-y4)^(2))*((x3-x4)^(2)+(y3-y4)^(2))
:Define r35_3=(((1)/(5))*√((8*x3+4*x4)^(2)+(8*y3+4*y4)^(2)))^(3)
:Define r45_3=(((1)/(5))*√((3*x3+9*x4)^(2)+(3*y3+9*y4)^(2)))^(3)
:Return rk23(system(vx3,((−4*(x3-x4))/(r34_3))-((8*x3+4*x4)/(r35_3)),vy3,((−4*(y3-y4))/(r34_3))-((8*y3+4*y4)/(r35_3)),vx4,((3*(x3-x4))/(r34_3))-((3*x3+9*x4)/(r45_3)),vy4,((3*(y3-y4))/(r34_3))-((3*y3+9*y4)/(r45_3))),t,{x3,vx3,y3,vy3,x4,vx4,y4,vy4},{0,dt},{x30,vx30,y30,vy30,x40,vx40,y40,vy40},dt,tol)
:EndFunc