数值分析导论：羽毛球的建模

介绍

羽毛球是最快的球类运动，尽管外表看起来很美。在这个例子中，我们将对羽毛球的飞行进行建模，并对其轨迹进行数值模拟。为此，我们将使用一个常微分方程解算器来追踪加速度回到轨迹。

import IPython
IPython.display.YouTubeVideo("wy5UraBy5co")

模型羽毛球飞行

羽毛球可以看作是一个准时的群众。它的加速度受牛顿第二定律支配：

\[m\vec a（m/r）=-mg\vec y-\frac 1 2 \rho v^2 a c x\vec t\]

进一步阅读：

羽毛球物理：http://www.worldadminton.com/reference/documents/5084354.pdf
阻力系数：https://en.wikipedia.org/wiki/Drag_系数

模拟输入：

%matplotlib notebook
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import integrate
import ipywidgets as ipw

v0  = 493. / 3.6 # Initial velocity [m/s]
A   = 4.e-3      # Shuttlecock cross area [m**2]
cx  = .62        # Drag coefficient []
m   = 4.e-2      # Shuttlecock mass [kg]
rho = 1.225      # Air density [kg/m**3]
g   = 9.81       # Gravity [m/s**2]

数值模拟

def derivative(X, t):
    """
    Target ODE: Newton's second law
    """
    x, y, vx, vy = X
    v = (vx**2 + vy**2)**.5
    Tx, Ty = vx / v, vy / v
    ax = -.5 * rho * v**2 * A * cx * Tx / m
    ay = -.5 * rho * v**2 * A * cx * Ty / m - g
    return np.array([vx, vy, ax, ay])


x0, y0 = 0., 0.
theta0 = 45.
X0 = [x0, y0, v0* np.cos(np.radians(theta0)), v0* np.sin(np.radians(theta0))]
t = np.linspace(0., 10., 200)
sol = integrate.odeint(derivative, X0, t)
out = pd.DataFrame(sol, columns = ["x", "y", "vx", "vy"])
out.head()

	x	y	vx	vy
0	0.000000	0.000000	96.834345	96.834345
1	4.323421	4.311938	76.793616	76.351638
2	7.831625	7.788303	63.660684	62.843454
3	10.785471	10.692467	54.394327	53.238940
4	13.337946	13.178880	47.510084	46.039129

plt.figure()
plt.plot(out.x, out.y)
plt.grid()
plt.ylim(0., 50.)
plt.xlabel("Position, $x$")
plt.ylabel("Position, $y$")
plt.show()

<IPython.core.display.Javascript object>

thetas = [0., 10.,15., 20., 30., 45., 60., 80., 85.]
plt.figure()

for theta0 in thetas:
    x0, y0 = 0., 3.
    X0 = [x0, y0, v0* np.cos(np.radians(theta0)), v0* np.sin(np.radians(theta0))]
    t = np.linspace(0., 10., 1000)
    sol = integrate.odeint(derivative, X0, t)
    out = pd.DataFrame(sol, columns = ["x", "y", "vx", "vy"])
    out["t"] = t
    plt.plot(out.x, out.y, label = r"$\theta_0 = $" + "{0}".format(theta0))
plt.legend()
plt.grid()
plt.ylim(0., 50.)
plt.xlabel("Position, $x$")
plt.ylabel("Position, $y$")
plt.show()

<IPython.core.display.Javascript object>

作为初始角度函数的范围 \(\theta_0\)

%%time
thetas = np.linspace(-180., 180., 300)
xmax = np.zeros_like(thetas)


for i in range(len(thetas)):
    theta0 = thetas[i]
    x0, y0 = 0., 3.
    X0 = [x0, y0, v0* np.sin(np.radians(theta0)), -v0* np.cos(np.radians(theta0))]
    t = np.linspace(0., 10., 10000)
    sol = integrate.odeint(derivative, X0, t)
    out = pd.DataFrame(sol, columns = ["x", "y", "vx", "vy"])
    xmax[i] = out[out.y < 0.].iloc[0].x

CPU times: user 1.17 s, sys: 0 ns, total: 1.17 s
Wall time: 1.24 s

plt.figure()
plt.plot(thetas, xmax)
plt.grid()
plt.xlabel(r"Start angle $\theta_0$")
plt.ylabel(r"Range $x_m$")
plt.show()

<IPython.core.display.Javascript object>

更具互动性的解决方案

N = 100




x = np.zeros(N)
y = x.copy()


plt.figure()
line, = plt.plot(x, y, "r-", label = "trajectory")
plt.grid()
plt.xlim(0., 80.)
plt.ylim(0., 80.)
plt.xlabel("Position, $x$ [m]")
plt.ylabel("Position, $y$ [m]")
plt.legend()

@ipw.interact(theta0 = (0.,90., 1.), v0 = (10., 500., 10))
def simulate(theta0 = 45., v0 = 490.):
    v0 /= 3.6
    x0, y0 = 0., 0.
    X0 = [x0, y0, v0* np.cos(np.radians(theta0)), v0* np.sin(np.radians(theta0))]
    t = np.linspace(0., 10, 200)
    sol = integrate.odeint(derivative, X0, t)
    line.set_xdata(sol[:, 0])
    line.set_ydata(sol[:, 1])

<IPython.core.display.Javascript object>

interactive(children=(FloatSlider(value=45.0, description='theta0', max=90.0, step=1.0), FloatSlider(value=490…