卓球とルンゲ＝クッタ法

ルンゲ＝クッタ法を用いて、卓球の速い打球の軌道をシミュレーションしてみる。

流体力学的効果を加味したモデルである。ただし、空気抵抗、マグヌス効果についてはこちらの文献を参考にした。

このシミュレーションは、対角線上のコースに、真の上回転をかけて、最大の速度で打ち出すモデルである。

コードはこちら。

#%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import math

#initial settings
N=300;dt=0.01;T=N*dt
x=[0]*N;y=[0]*N;z=[0]*N
dx=dy=dz=0.01
vx=[0]*N;vy=[0]*N;vz=[0]*N
#parameter
g=9.8;r=0.02;A=math.pi*(r**2);m=0.0027
h=0.1525;lx=1.525;ly=2.74
v=26
theta=math.pi/7.2
vh=v*math.cos(theta)
vv=v*math.sin(theta)
CD=0.5;CL=0.35;rho=1.293
x[0]=0;y[0]=0;z[0]=h
vx[0]=vh*lx/(lx**2+ly**2)**(1/2)
vy[0]=vh*ly/(lx**2+ly**2)**(1/2)
vz[0]=vv

def vnorm(a,b,c):
  return (a**2+b**2+c**2)

def Fd(a,b,c,d):
  return -1/2*CD*rho*A*vnorm(a,b,c)**(1/2)*d/m

def Flh(a,b,c,d):
  return 1/2*CL*rho*A*(vnorm(a,b,c)**(1/2))*c*(d/(a**2+b**2)**(1/2))/m
def Flv(a,b,c):
  return 1/2*CL*rho*A*(vnorm(a,b,c)**(1/2))*(a**2+b**2)**(1/2)/m


for i in range(N-1):
  xk1=vx[i]*dt;yk1=vy[i]*dt;zk1=vz[i]*dt
  vxk1=Fd(vx[i],vy[i],vz[i],vx[i])*dt+Flh(vx[i],vy[i],vz[i],vx[i])*dt
  vyk1=Fd(vx[i],vy[i],vz[i],vy[i])*dt+Flh(vx[i],vy[i],vz[i],vy[i])*dt
  vzk1=Fd(vx[i],vy[i],vz[i],vz[i])*dt-g*dt-Flv(vx[i],vy[i],vz[i])*dt
  xk2=(vx[i]+vxk1/2)*dt;yk2=(vy[i]+vyk1/2)*dt;zk2=(vz[i]+vzk1/2)*dt
  vxk2=Fd(vx[i]+vxk1/2,vy[i]+vyk1/2,vz[i]+vzk1/2,vx[i]+vxk1/2)*dt+Flh(vx[i]+vxk1/2,vy[i]+vyk1/2,vz[i]+vzk1/2,vx[i]+vxk1/2)*dt
  vyk2=Fd(vx[i]+vxk1/2,vy[i]+vyk1/2,vz[i]+vzk1/2,vy[i]+vyk1/2)*dt+Flh(vx[i]+vxk1/2,vy[i]+vyk1/2,vz[i]+vzk1/2,vy[i]+vyk1/2)*dt
  vzk2=Fd(vx[i]+vxk1/2,vy[i]+vyk1/2,vz[i]+vzk1/2,vz[i]+vzk1/2)*dt-g*dt-Flv(vx[i]+vxk1/2,vy[i]+vyk1/2,vz[i]+vzk1/2)*dt
  xk3=(vx[i]+vxk2/2)*dt;yk3=(vy[i]+vyk2/2)*dt;zk3=(vz[i]+vzk2/2)*dt
  vxk3=Fd(vx[i]+vxk2/2,vy[i]+vyk2/2,vz[i]+vzk2/2,vx[i]+vxk2/2)*dt+Flh(vx[i]+vxk2/2,vy[i]+vyk2/2,vz[i]+vzk2/2,vx[i]+vxk2/2)*dt
  vyk3=Fd(vx[i]+vxk2/2,vy[i]+vyk2/2,vz[i]+vzk2/2,vy[i]+vyk2/2)*dt+Flh(vx[i]+vxk2/2,vy[i]+vyk2/2,vz[i]+vzk2/2,vy[i]+vyk2/2)*dt
  vzk3=Fd(vx[i]+vxk2/2,vy[i]+vyk2/2,vz[i]+vzk2/2,vz[i]+vzk2/2)*dt-g*dt-Flv(vx[i]+vxk2/2,vy[i]+vyk2/2,vz[i]+vzk2/2)*dt
  xk4=(vx[i]+vxk3)*dt;yk4=(vy[i]+vyk3)*dt;zk4=(vz[i]+vzk3)*dt
  vxk4=Fd(vx[i]+vxk3,vy[i]+vyk3,vz[i]+vzk3,vx[i]+vxk3)*dt+Flh(vx[i]+vxk3,vy[i]+vyk3,vz[i]+vzk3,vx[i]+vxk3)*dt
  vyk4=Fd(vx[i]+vxk3,vy[i]+vyk3,vz[i]+vzk3,vy[i]+vyk3)*dt+Flh(vx[i]+vxk3,vy[i]+vyk3,vz[i]+vzk3,vy[i]+vyk3)*dt
  vzk4=Fd(vx[i]+vxk3,vy[i]+vyk3,vz[i]+vzk3,vz[i]+vzk3)*dt-g*dt-Flv(vx[i]+vxk3,vy[i]+vyk3,vz[i]+vzk3)*dt
  x[i+1]=x[i]+(xk1+xk2*2+xk3*2+xk4)/6
  y[i+1]=y[i]+(yk1+yk2*2+yk3*2+yk4)/6
  z[i+1]=z[i]+(zk1+zk2*2+zk3*2+zk4)/6
  vx[i+1]=vx[i]+(vxk1+vxk2*2+vxk3*2+vxk4)/6
  vy[i+1]=vy[i]+(vyk1+vyk2*2+vyk3*2+vyk4)/6
  vz[i+1]=vz[i]+(vzk1+vzk2*2+vzk3*2+vzk4)/6


ground=[0]*N
net=[ly/2]*30
high=[h*i/30 for i in range(30)]
target=[ly/2*i/N for i in range(N)]
distance=[ly*i/N for i in range(N)]
# xに対応するyの値を用意
#plt.plot(x,y,"r-..")
plt.plot(y,z,"g")
plt.plot(distance,ground,"b+")
plt.plot(target,ground,"r.")
plt.plot(net,high,"r.")
#plt.plot(time,y3,"b--")
#plt.plot(time,y4,"purple")
plt.xlabel("horizontal")
plt.ylabel("vertical")
plt.title("pingpong attack")
plt.xlim(-0.2,3.13)
plt.ylim(-0.76,1)
plt.axes().set_aspect('equal')
# plt.show()で画面に表示
plt.show()