为数值求解的微分方程实现初始条件

Question

想象某人在一定角度theta 和速度v0 下从阳台上跳下来（阳台的高度表示为ystar）。在 2D 中查看这个问题并考虑阻力，你会得到一个微分方程组，可以用 Runge-Kutta 方法求解（我选择显式中点，不确定这个屠夫表是什么）。我实现了这个并且它工作得非常好，对于一些给定的初始条件，我得到了移动粒子的轨迹。

我的问题是我想修复两个初始条件（x 轴上的起点为零，y 轴上的起点为ystar）并确保轨迹经过某个点x 轴（我们称之为xstar）。为此当然存在其他两个初始条件的多种组合，在这种情况下是 x 和 y 方向上的速度。问题是我不知道如何实现它。

到目前为止我用来解决问题的代码：

1) Runge-Kutta 方法的实现

import numpy as np 
import matplotlib.pyplot as plt

def integrate(methode_step, rhs, y0, T, N):
    star = (int(N+1),y0.size)
    y= np.empty(star)
    t0, dt = 0, 1.* T/N
    y[0,...] = y0
    for i in range(0,int(N)):
        y[i+1,...]=methode_step(rhs,y[i,...], t0+i*dt, dt)
    t = np.arange(N+1) * dt
    return t,y

def explicit_midpoint_step(rhs, y0, t0, dt):
    return y0 + dt * rhs(t0+0.5*dt,y0+0.5*dt*rhs(t0,y0))

def explicit_midpoint(rhs,y0,T,N):
    return integrate(explicit_midpoint_step,rhs,y0,T,N)

2) 微分方程右边的实现和需要的参数

A = 1.9/2.
cw = 0.78
rho = 1.293
g = 9.81

# Mass and referece length
l = 1.95
m = 118

# Position
xstar = 8*l
ystar = 4*l

def rhs(t,y):
    lam = cw * A * rho /(2 * m)
    return np.array([y[1],-lam*y[1]*np.sqrt(y[1]**2+y[3]**2),y[3],-lam*y[3]*np.sqrt(y[1]**2+y[3]**2)-g])

3) 用它解决问题

# Parametrize the two dimensional velocity with an angle theta and speed v0

v0 = 30
theta = np.pi/6

v0x = v0 * np.cos(theta)
v0y = v0 * np.sin(theta)

# Initial condintions
z0 = np.array([0, v0x, ystar, v0y])

# Calculate solution
t, z = explicit_midpoint(rhs, z0, 5, 1000)

4) 可视化

plt.figure()
plt.plot(0,ystar,"ro")
plt.plot(x,0,"ro")
plt.plot(z[:,0],z[:,1])
plt.grid(True)
plt.xlabel(r"$x$")
plt.ylabel(r"$y$")
plt.show()

为了使问题具体化：考虑到这个设置，我如何找到v0 和theta 的所有可能组合，使得z[some_element,0]==xstar

我当然尝试了一些东西，主要是修复 theta 的蛮力方法，然后尝试所有可能的速度（以有意义的间隔），但最后不知道如何将结果数组与想要的结果……

由于这主要是一个编码问题，我希望堆栈溢出是寻求帮助的正确地方...

编辑： 按照这里的要求，我尝试从上面解决问题（替换 3）和 4）..

theta = np.pi/4.
xy = np.zeros((50,1001,2))
z1 = np.zeros((1001,2))
count=0

for v0 in range(0,50):
    v0x = v0 * np.cos(theta)
    v0y = v0 * np.sin(theta)
    z0 = np.array([0, v0x, ystar, v0y])

    # Calculate solution
    t, z = explicit_midpoint(rhs, z0, 5, 1000)    

    if np.around(z[:,0],3).any() == round(xstar,3):
        z1[:,0] = z[:,0]
        z1[:,1] = z[:,2]
        break
    else:
        xy[count,:,0] = z[:,0]
        xy[count,:,1] = z[:,2]
        count+=1


plt.figure()
plt.plot(0,ystar,"ro")
plt.plot(xstar,0,"ro")

for k in range(0,50):    
    plt.plot(xy[k,:,0],xy[k,:,1])

plt.plot(z[:,0],z[:,1])
plt.grid(True)
plt.xlabel(r"$x$")
plt.ylabel(r"$y$")
plt.show()

我确定我使用了错误的.any() 函数，我的想法是将z[:,0] 的值四舍五入到三位数字，然后将它们与xstar 进行比较，如果匹配，则循环应该终止并重新运行当前的z，如果不是，则应将其保存在另一个数组中，然后增加v0。

Answer 1

所以经过一番尝试后，我找到了实现我想要的方法......这是我在开始的帖子中提到的蛮力方法，但至少现在它有效......

这个想法很简单：定义一个函数find_v0，它可以找到给定的theta 和v0。在此函数中，您为v0 取一个起始值（我选择 8，但这只是我的猜测），然后取起始值并使用difference 函数检查有趣点与(xstar,0) 的距离。 .在这种情况下，有趣的点可以通过将 x 轴上大于 xstar 的所有点设置为零（及其相应的 y 值）然后用 trim_zeros 修剪所有零点来确定，现在是最后一个的元素对应于所需的输出。如果差异函数的输出小于临界值（在我的情况下为 0.1），则将当前的 v0 传递给打开，如果不是，则将其增加 0.01 并再次执行相同的操作。

此代码如下所示（再次替换 3) 和 4)）：

th = np.linspace(0,np.pi/3,100)

def find_v0(theta):
    v0=8
    while(True):
        v0x = v0 * np.cos(theta)
        v0y = v0 * np.sin(theta)
        z0 = np.array([0, v0x, ystar, v0y])

        # Calculate solution
        t, z = explicit_midpoint(rhs, z0, 5, 1000)    

        for k in range(1001):
            if z[k,0] > xstar: 
                z[k,0] = 0
                z[k,2] = 0

        x = np.trim_zeros(z[:,0])
        y = np.trim_zeros(z[:,2])

        diff = difference(x[-1],y[-1])

        if diff < 0.1:
            break
        else: v0+=0.01
    return v0#,x,y[0:]

v0 = np.zeros_like(th)

from tqdm import tqdm

count=0
for k in tqdm(th):
    v0[count] = find_v0(k)
    count+=1

v0_interp = interpolate.interp1d(th,v0)

plt.figure()
plt.plot(th,v0_interp(th),"g")
plt.grid(True)
plt.xlabel(r"$\theta$")
plt.ylabel(r"$v_0$")
plt.show()

这个东西的问题是计算需要很长时间（当前设置大约需要 5-6 分钟）。如果有人有一些提示如何改进代码以获得更快的速度或有不同的方法，仍然会受到赞赏。

Answer 2

假设x方向的速度永远不会下降到零，您可以将x作为独立参数而不是时间。然后状态向量是时间、位置、速度，并且该状态空间中的向量场被缩放，使得 vx 分量始终为 1。然后从零积分到 xstar 以计算轨迹与 xstar 相交的状态（近似值）为 x-价值。

def derivs(u,x):
    t,x,y,vx,vy = u
    v = hypot(vx,vy)
    ax = -lam*v*vx
    ay = -lam*v*vy - g
    return [ 1/vx, 1, vy/vx, ax/vx, ay/vx ]

odeint(derivs, [0, x0, y0, vx0, vy0], [0, xstar])

或使用您自己的集成方法。我使用 odeint 作为文档接口来展示如何在集成中使用这个衍生函数。

生成的时间和 y 值可能是极端的

Answer 3

编辑 2018-07-16

考虑到空中阻力，我在这里发布了一个更正的答案。

下面是一个 python 脚本，用于计算 (v0,theta) 值的集合，以便空中拖动的轨迹在某个时间通过 (x,y) = (xstar,0) t=tstar。我使用没有空气阻力的轨迹作为初始猜测，并猜测x(tstar) 对v0 的依赖性进行第一次改进。到达正确的v0 所需的迭代次数通常为 3 到 4 次。脚本在我的笔记本电脑上在 0.99 秒内完成，包括生成数字的时间。

脚本生成两个图形和一个文本文件。

• fig_xdrop_v0_theta.png
  • 黑点表示解集(v0,theta)
  • 黄线表示参考(v0,theta)，如果没有空气阻力，这将是一个解决方案。
• fig_traj_sample.png
  • 当从解集中采样(v0,theta) 时，检查轨迹（蓝色实线）是否通过(x,y)=(xstar,0)。
  • 黑色虚线表示没有空气阻力的轨迹作为参考。
• output.dat
  • 包含(v0,theta)的数值数据以及tstar的着陆时间和找到v0所需的迭代次数。

这里开始脚本。

#!/usr/bin/env python3

import numpy as np
import scipy.integrate
import matplotlib as mpl
import matplotlib.pyplot as plt
import matplotlib.image as img

mpl.rcParams['lines.linewidth'] = 2
mpl.rcParams['lines.markeredgewidth'] = 1.0
mpl.rcParams['axes.formatter.limits'] = (-4,4)
#mpl.rcParams['axes.formatter.limits'] = (-2,2)
mpl.rcParams['axes.labelsize'] = 'large'
mpl.rcParams['xtick.labelsize'] = 'large'
mpl.rcParams['ytick.labelsize'] = 'large'
mpl.rcParams['xtick.direction'] = 'out'
mpl.rcParams['ytick.direction'] = 'out'



############################################
len_ref = 1.95
xstar = 8.0*len_ref
ystar = 4.0*len_ref
g_earth = 9.81
#
mass = 118
area = 1.9/2.
cw = 0.78
rho = 1.293
lam = cw * area * rho /(2.0 * mass)
############################################
ngtheta=51
theta_min = -0.1*np.pi
theta_max =  0.4*np.pi
theta_grid = np.linspace(theta_min, theta_max, ngtheta)
#
ngv0=100
v0min =6.0
v0max =18.0
v0_grid=np.linspace(v0min, v0max, ngv0)
# .. this grid is used for the initial coarse scan by reference trajecotry
############################################
outf=open('output.dat','w')
print('data file generated: output.dat')
###########################################


def calc_tstar_ref_and_x_ref_at_tstar_ref(v0, theta, ystar, g_earth):
    '''return the drop time t* and drop point x(t*) of a reference trajectory
    without air drag.
    '''

    vx = v0*np.cos(theta)
    vy = v0*np.sin(theta)
    ts_ref = (vy+np.sqrt(vy**2+2.0*g_earth*ystar))/g_earth
    x_ref  = vx*ts_ref 
    return (ts_ref, x_ref)


def rhs_drag(yvec, time, g_eath, lamb):
    '''
    dx/dt = v_x
    dy/dt = v_y
    du_x/dt = -lambda v_x sqrt(u_x^2 + u_y^2)
    du_y/dt = -lambda v_y sqrt(u_x^2 + u_y^2) -g

    yvec[0] .. x
    yvec[1] .. y
    yvec[2] .. v_x
    yvec[3] .. v_y
    '''

    vnorm = (yvec[2]**2+yvec[3]**2)**0.5
    return  [ yvec[2], yvec[3], -lamb*yvec[2]*vnorm, -lamb*yvec[3]*vnorm -g_earth]


def try_tstar_drag(v0, theta, ystar, g_earth, lamb, tstar_search_grid):
    '''one trial run to find the drop point x(t*), y(t*) of a trajectory 
    under the air drag.
    '''

    tinit=0.0
    tgrid = [tinit]+list(tstar_search_grid)
    yvec_list = scipy.integrate.odeint(rhs_drag,
                                       [0.0, ystar, v0*np.cos(theta), v0*np.sin(theta)],
                                       tgrid, args=(g_earth, lam))
    y_drag = [yvec[1] for yvec in yvec_list]
    x_drag = [yvec[0] for yvec in yvec_list]

    if y_drag[0]<0.0:
        ierr=-1
        jtstar=0
        tstar_braket=None
    elif y_drag[-1]>0.0:
        ierr=1
        jtstar=len(y_drag)-1
        tstar_braket=None
    else:
        ierr=0
        for jt in range(len(y_drag)-1):
            if y_drag[jt+1]*y_drag[jt]<=0.0:
                tstar_braket=[tgrid[jt],tgrid[jt+1]]
                if abs(y_drag[jt+1])<abs(y_drag[jt]):
                    jtstar = jt+1
                else:
                    jtstar = jt
                break

    tstar_est = tgrid[jtstar]
    x_drag_at_tstar_est = x_drag[jtstar]
    y_drag_at_tstar_est = y_drag[jtstar]
    return (tstar_est, x_drag_at_tstar_est, y_drag_at_tstar_est, ierr, tstar_braket)


def calc_x_drag_at_tstar(v0, theta, ystar, g_earth, lamb, tstar_est,
                         eps_y=1.0e-3, ngt_search=20,
                         rel_range_lower=0.8, rel_range_upper=1.2,
                         num_try=5):
    '''compute the dop point x(t*) of a trajectory under the air drag.
    '''

    flg_success=False
    tstar_est_lower=tstar_est*rel_range_lower
    tstar_est_upper=tstar_est*rel_range_upper
    for jtry in range(num_try):
        tstar_search_grid = np.linspace(tstar_est_lower, tstar_est_upper, ngt_search)
        tstar_est, x_drag_at_tstar_est, y_drag_at_tstar_est, ierr, tstar_braket \
            = try_tstar_drag(v0, theta, ystar, g_earth, lamb, tstar_search_grid)
        if ierr==-1:
            tstar_est_upper = tstar_est_lower
            tstar_est_lower = tstar_est_lower*rel_range_lower
        elif ierr==1:
            tstar_est_lower = tstar_est_upper
            tstar_est_upper = tstar_est_upper*rel_range_upper
        else:
            if abs(y_drag_at_tstar_est)<eps_y:
                flg_success=True
                break
            else:
                tstar_est_lower=tstar_braket[0]
                tstar_est_upper=tstar_braket[1]

    return (tstar_est, x_drag_at_tstar_est, y_drag_at_tstar_est, flg_success)


def find_v0(xstar, v0_est, theta, ystar, g_earth, lamb, tstar_est,
            eps_x=1.0e-3, num_try=6):
   '''solve for v0 so that x(t*)==x*.
   '''

   flg_success=False
   v0_hist=[]
   x_drag_at_tstar_hist=[]
   jtry_end=None
   for jtry in range(num_try):
       tstar_est, x_drag_at_tstar_est, y_drag_at_tstar_est, flg_success_x_drag \
           = calc_x_drag_at_tstar(v0_est, theta, ystar, g_earth, lamb, tstar_est)
       v0_hist.append(v0_est)
       x_drag_at_tstar_hist.append(x_drag_at_tstar_est)
       if not flg_success_x_drag:
           break
       elif abs(x_drag_at_tstar_est-xstar)<eps_x: 
           flg_success=True
           jtry_end=jtry
           break
       else:
           # adjust v0
           # better if tstar_est is also adjusted, but maybe that is too much.
           if len(v0_hist)<2:
               # This is the first run. Use the analytical expression of
               # dx(tstar)/dv0 of  the refernece trajectory
               dx = xstar - x_drag_at_tstar_est
               dv0 = dx/(tstar_est*np.cos(theta))
               v0_est += dv0
           else:
               # use linear interpolation
               v0_est = v0_hist[-2] \
                        + (v0_hist[-1]-v0_hist[-2]) \
                        *(xstar                 -x_drag_at_tstar_hist[-2])\
                        /(x_drag_at_tstar_hist[-1]-x_drag_at_tstar_hist[-2])

   return (v0_est, tstar_est, flg_success, jtry_end)        


# make a reference table of t* and x(t*) of a trajectory without air drag
# as a function of v0 and theta.
tstar_ref=np.empty((ngtheta,ngv0))
xdrop_ref=np.empty((ngtheta,ngv0))
for j1 in range(ngtheta):
    for j2 in range(ngv0):
        tt, xx = calc_tstar_ref_and_x_ref_at_tstar_ref(v0_grid[j2], theta_grid[j1], ystar, g_earth)
        tstar_ref[j1,j2] = tt
        xdrop_ref[j1,j2] = xx

# make an estimate of v0 and t* of a dragged trajectory for each theta 
# based on the reference trajectroy's landing position xdrop_ref.
tstar_est=np.empty((ngtheta,))
v0_est=np.empty((ngtheta,))
v0_est[:]=-1.0
# .. null value
for j1 in range(ngtheta):
    for j2 in range(ngv0-1):
        if (xdrop_ref[j1,j2+1]-xstar)*(xdrop_ref[j1,j2]-xstar)<=0.0:
            tstar_est[j1] = tstar_ref[j1,j2]
            # .. lazy
            v0_est[j1] \
                = v0_grid[j2] \
                + (v0_grid[j2+1]-v0_grid[j2])\
                *(xstar-xdrop_ref[j1,j2])/(xdrop_ref[j1,j2+1]-xdrop_ref[j1,j2])
            # .. linear interpolation
            break

print('compute v0 for each theta under air drag..')
# compute v0 for each theta under air drag
theta_sol_list=[]
tstar_sol_list=[]
v0_sol_list=[]
outf.write('#  theta  v0  tstar  numiter_v0\n')
for j1 in range(ngtheta):
    if v0_est[j1]>0.0:
        v0, tstar, flg_success, jtry_end \
            = find_v0(xstar, v0_est[j1], theta_grid[j1], ystar, g_earth, lam, tstar_est[j1])
        if flg_success:
            theta_sol_list.append(theta_grid[j1])
            v0_sol_list.append(v0)
            tstar_sol_list.append(tstar)
            outf.write('%26.16e  %26.16e %26.16e %10i\n'
                       %(theta_grid[j1], v0, tstar, jtry_end+1))

theta_sol = np.array(theta_sol_list)            
v0_sol    = np.array(v0_sol_list)            
tstar_sol = np.array(tstar_sol_list)            



### Check a sample 
jsample=np.size(v0_sol)//3
theta_sol_sample= theta_sol[jsample]
v0_sol_sample   = v0_sol[jsample]
tstar_sol_sample= tstar_sol[jsample]

ngt_chk = 50
tgrid  = np.linspace(0.0, tstar_sol_sample, ngt_chk)

yvec_list = scipy.integrate.odeint(rhs_drag,
                                   [0.0, ystar,
                                    v0_sol_sample*np.cos(theta_sol_sample),
                                    v0_sol_sample*np.sin(theta_sol_sample)],
                                   tgrid, args=(g_earth, lam))
y_drag_sol_sample = [yvec[1] for yvec in yvec_list]
x_drag_sol_sample = [yvec[0] for yvec in yvec_list]

# compute also the trajectory without drag starting form the same initial
# condiiton by setting lambda=0.
yvec_list = scipy.integrate.odeint(rhs_drag,
                                   [0.0, ystar,
                                    v0_sol_sample*np.cos(theta_sol_sample),
                                    v0_sol_sample*np.sin(theta_sol_sample)],
                                   tgrid, args=(g_earth, 0.0))
y_ref_sample = [yvec[1] for yvec in yvec_list]
x_ref_sample = [yvec[0] for yvec in yvec_list]


#######################################################################
# canvas setting 
#######################################################################
f_size   = (8,5)
#
a1_left   = 0.15
a1_bottom  = 0.15
a1_width  = 0.65
a1_height = 0.80
#
hspace=0.02
#
ac_left   = a1_left+a1_width+hspace
ac_bottom = a1_bottom
ac_width  = 0.03
ac_height = a1_height
###########################################



############################################
# plot 
############################################

#------------------------------------------------
print('plotting the solution..')
fig1=plt.figure(figsize=f_size)
ax1 =plt.axes([a1_left, a1_bottom, a1_width, a1_height], axisbg='w')

im1=img.NonUniformImage(ax1,
                            interpolation='bilinear', \
                            cmap=mpl.cm.Blues, \
                            norm=mpl.colors.Normalize(vmin=0.0,
                                                          vmax=np.max(xdrop_ref), clip=True)) 
im1.set_data(v0_grid, theta_grid/np.pi, xdrop_ref )
ax1.images.append(im1)
plt.contour(v0_grid, theta_grid/np.pi, xdrop_ref, [xstar], colors='y')
plt.plot(v0_sol, theta_sol/np.pi, 'ok', lw=4, label='Init Cond with Drag')
plt.legend(loc='lower left')

plt.xlabel(r'Initial Velocity $v_0$', fontsize=18)
plt.ylabel(r'Angle of Projection $\theta/\pi$', fontsize=18)
plt.yticks([-0.50, -0.25, 0.0, 0.25, 0.50])
ax1.set_xlim([v0min, v0max])
ax1.set_ylim([theta_min/np.pi, theta_max/np.pi])

axc =plt.axes([ac_left, ac_bottom, ac_width, ac_height], axisbg='w')
mpl.colorbar.Colorbar(axc,im1)
axc.set_ylabel('Distance from Blacony without Drag')
# 'Distance from Blacony $x(t^*)$'

plt.savefig('fig_xdrop_v0_theta.png')
print('figure file genereated: fig_xdrop_v0_theta.png')
plt.close()


#------------------------------------------------
print('plotting a sample trajectory..')
fig1=plt.figure(figsize=f_size)
ax1 =plt.axes([a1_left, a1_bottom, a1_width, a1_height], axisbg='w')
plt.plot(x_drag_sol_sample, y_drag_sol_sample, '-b', lw=2, label='with drag')
plt.plot(x_ref_sample, y_ref_sample, '--k', lw=2, label='without drag')
plt.axvline(x=xstar, color=[0.3, 0.3, 0.3], lw=1.0)
plt.axhline(y=0.0, color=[0.3, 0.3, 0.3], lw=1.0)
plt.legend()
plt.text(0.1*xstar, 0.6*ystar,
         r'$v_0=%5.2f$'%(v0_sol_sample)+'\n'+r'$\theta=%5.2f \pi$'%(theta_sol_sample/np.pi),
         fontsize=18)
plt.text(xstar, 0.5*ystar, 'xstar', fontsize=18)

plt.xlabel(r'Horizontal Distance $x$', fontsize=18)
plt.ylabel(r'Height $y$', fontsize=18)
ax1.set_xlim([0.0,        1.5*xstar])
ax1.set_ylim([-0.1*ystar, 1.5*ystar])

plt.savefig('fig_traj_sample.png')
print('figure file genereated: fig_traj_sample.png')
plt.close()

outf.close()

这是fig_xdrop_v0_theta.png的数字。

这是fig_traj_sample.png的数字。

---

编辑 2018-07-15

我意识到我忽略了这个问题考虑了空气阻力。我真丢脸。所以，我下面的回答是不正确的。恐怕我自己删除我的答案看起来像是在隐藏一个错误，我现在把它留在下面。如果人们认为不正确的答案在周围徘徊很烦人，我没关系。有人删除它。

---

微分方程其实可以用手解， 它不需要太多的计算资源 绘制出该人离阳台的距离 在地面上作为初始速度v0 和 角度theta。然后，您可以选择条件(v0,theta) 这样distance_from_balcony_on_the_ground(v0,theta) = xstar 从这个数据表。

让我们写出横坐标和纵坐标 时间t 的人分别是x(t) 和y(t)。 我想你在大楼的墙上拿了x=0 和y=0 作为地面层，我也在这里这样做。让我们说 人在时间的水平和垂直速度t 分别是v_x(t) 和v_y(t)。 t=0 的初始条件为

x(0) = 0
y(0) = ystar
v_x(0) = v0 cos theta
v_y(0) = v0 sin theta

你正在求解的牛顿方程是，

dx/dt = v_x  .. (1)
dy/dt = v_y  .. (2)
m d v_x /dt = 0  .. (3)
m d v_y /dt = -m g  .. (4)

m 是人​​的质量， 而g 是我不知道英文名字的常量， 但我们都知道它是什么。

从等式。 (3)、

v_x(t) = v_x(0) = v0 cos theta.

将其与 eq 一起使用。 (1)、

x(t) = x(0) + \int_0^t dt' v_x(t') = t v0 cos theta,

我们也使用了初始条件。 \int_0^t 表示 从0 到t 的积分。

从等式。 (4)、

v_y(t)
= v_y (0) + \int_0^t dt' (-g) 
= v0 sin theta -g t,

我们使用了初始条件。 将其与 eq 一起使用。 (3) 并且也使用初始条件，

y(t)
= y(0) + \int_0^t dt' v_y(t')
= ystar + t v0 sin theta -t^2 (g/2).

其中t^2 表示t 的平方。 从y(t) 的表达式，我们可以得到时间tstar 该人撞到地面的位置。即y(tstar) =0。 这个方程可以用二次公式求解 （或类似名称）为

tstar = (v0 sin theta + sqrt((v0 sin theta)^2 + 2g ystar)/g,

我使用了条件tstar>0。现在我们知道了 人击中时离阳台的距离 地面为x(tstar)。使用上面x(t) 的表达式，

x(tstar) = (v0 cos theta) (v0 sin theta + sqrt((v0 sin theta)^2 + 2g ystar))/g.
  .. (5)

实际上x(tstar) 依赖于v0 和theta 以及g 和ystar。 你持有g 和ystar 作为常量，你想找到 所有(v0,theta) 使得x(tstar) = xstar 对于给定的xstar 值。

由于等式的右侧。 (5) 可以廉价地计算， 您可以为v0 和theta 设置网格并计算xstar 在这个二维网格上。然后，您可以看到解决方案集大致在哪里 (v0,theta) 的谎言。如果您需要精确的解决方案，您可以拿起 包含此数据表中的解决方案的段。

下面是一个演示这个想法的python脚本。

我还附上了这个脚本生成的图。 黄色曲线是解集(v0,theta)，使得 人从墙上撞到xstar 当您设置xstar = 8.0*1.95 和ystar=4.0*1.95 时。 蓝色坐标表示x(tstar)，即距离 人从阳台上水平跳下。 请注意，在给定的v0（高于v0=9.9 的阈值）， 有两个 theta 值（人的两个方向 投射自己）到达目标点(x,y) = (xstar,0)。 theta 值的较小分支可以是负数，这意味着只要初始速度足够高，人就可以向下跳跃到达目标点。

该脚本还会生成一个数据文件output.dat，其中包含 解决方案封闭段。

#!/usr/bin/python3

import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import matplotlib.image as img

mpl.rcParams['lines.linewidth'] = 2
mpl.rcParams['lines.markeredgewidth'] = 1.0
mpl.rcParams['axes.formatter.limits'] = (-4,4)
#mpl.rcParams['axes.formatter.limits'] = (-2,2)
mpl.rcParams['axes.labelsize'] = 'large'
mpl.rcParams['xtick.labelsize'] = 'large'
mpl.rcParams['ytick.labelsize'] = 'large'
mpl.rcParams['xtick.direction'] = 'out'
mpl.rcParams['ytick.direction'] = 'out'


############################################
len_ref = 1.95
xstar = 8.0*len_ref
ystar = 4.0*len_ref
g_earth = 9.81
############################################
ngv0=100
v0min =0.0
v0max =20.0
v0_grid=np.linspace(v0min, v0max, ngv0)
############################################
outf=open('output.dat','w')
print('data file generated: output.dat')
###########################################

def x_at_tstar(v0, theta, ystar, g_earth):
    vx = v0*np.cos(theta)
    vy = v0*np.sin(theta)
    return (vy+np.sqrt(vy**2+2.0*g_earth*ystar))*vx/g_earth

ngtheta=100
theta_min = -0.5*np.pi
theta_max =  0.5*np.pi
theta_grid = np.linspace(theta_min, theta_max, ngtheta)
xdrop=np.empty((ngv0,ngtheta))
# x(t*) as a function of v0 and theta.
for j1 in range(ngv0):
    for j2 in range(ngtheta):
        xdrop[j1,j2] = x_at_tstar(v0_grid[j1], theta_grid[j2], ystar, g_earth)

outf.write('# domain [theta_lower, theta_upper] that encloses the solution\n')        
outf.write('# theta such that x_at_tstart(v0,theta, ystart, g_earth)=xstar\n')        
outf.write('# v0  theta_lower  theta_upper  x_lower  x_upper\n')
for j1 in range(ngv0):
    for j2 in range(ngtheta-1):
       if (xdrop[j1,j2+1]-xstar)*(xdrop[j1,j2]-xstar)<=0.0:
           outf.write('%26.16e %26.16e %26.16e %26.16e %26.16e\n'
                          %(v0_grid[j1], theta_grid[j2], theta_grid[j2+1],
                                xdrop[j1,j2], xdrop[j1,j2+1]))
print('See output.dat for the segments enclosing solutions.')
print('You can hunt further for precise solutions using this data.')



#######################################################################
# canvas setting 
#######################################################################
f_size   = (8,5)
#
a1_left   = 0.15
a1_bottom  = 0.15
a1_width  = 0.65
a1_height = 0.80
#
hspace=0.02
#
ac_left   = a1_left+a1_width+hspace
ac_bottom = a1_bottom
ac_width  = 0.03
ac_height = a1_height
###########################################



############################################
# plot 
############################################

print('plotting..')
fig1=plt.figure(figsize=f_size)
ax1 =plt.axes([a1_left, a1_bottom, a1_width, a1_height], axisbg='w')

im1=img.NonUniformImage(ax1,
                            interpolation='bilinear', \
                            cmap=mpl.cm.Blues, \
                            norm=mpl.colors.Normalize(vmin=0.0,
                                                          vmax=np.max(xdrop), clip=True)) 
im1.set_data(v0_grid, theta_grid/np.pi, np.transpose(xdrop))
ax1.images.append(im1)
plt.contour(v0_grid, theta_grid/np.pi, np.transpose(xdrop), [xstar], colors='y')

plt.xlabel(r'Initial Velocity $v_0$', fontsize=18)
plt.ylabel(r'Angle of Projection $\theta/\pi$', fontsize=18)
plt.yticks([-0.50, -0.25, 0.0, 0.25, 0.50])
ax1.set_xlim([v0min, v0max])
ax1.set_ylim([theta_min/np.pi, theta_max/np.pi])

axc =plt.axes([ac_left, ac_bottom, ac_width, ac_height], axisbg='w')
mpl.colorbar.Colorbar(axc,im1)
# 'Distance from Blacony $x(t^*)$'

plt.savefig('fig_xdrop_v0_theta.png')
print('figure file genereated: fig_xdrop_v0_theta.png')
plt.close()



outf.close()