想象一下有人在一定的角度下跳下阳台 theta 和速度 v0 （阳台的高度表示为 ystar ） . 在2D中考虑这个问题并考虑拖动，你会得到一个微分方程系统，可以用Runge-Kutta方法求解（我选择显式中点，不确定这个的屠夫 table 是什么） . 我实现了它并且它完美地工作，对于一些给定的初始条件，我得到了移动粒子的轨迹 .
我的问题是我想修复两个初始条件（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）实现微分方程的右侧和nessecery参数
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 然后尝试所有可能的速度（在一个有意义的intervall中）但最后不知道如何将得到的数组与所需的结果进行比较......
由于这主要是编码问题，我希望堆栈溢出是寻求帮助的正确位置...
编辑：这里要求的是我尝试从上面解决问题（替换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()
我'm sure that I' m使用 .any() 函数错误，想法是将 z[:,0] 的值舍入为三位数，然后将它们与 xstar 进行比较，如果匹配则循环应终止并重新启动当前的 z ，如果不是，则应将其保存在另一个中数组然后增加 v0 .
3 回答
编辑2018-07-16
在这里，我发布了一个纠正的答案，考虑到空中阻力 .
下面是一个python脚本，用于计算
(v0,theta)值的集合，以便空气拖动的轨迹在某个时间t=tstar通过(x,y) = (xstar,0). 我使用没有空气拖曳的轨迹作为初始猜测，并且还猜测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所需的迭代次数 .这里开始脚本 .
这是图
fig_xdrop_v0_theta.png.这是图
fig_traj_sample.png.编辑2018-07-15
我意识到我忽略了这个问题考虑了空气阻力 . 真可惜我 . 所以，我的答案不正确 . 我担心自己删除我的答案看起来像隐瞒了一个错误，我暂时将其留在下面 . 如果人们认为这是一个令人讨厌的错误答案，我就是O.K.有人删除它 .
微分方程实际上可以手工求解，并且它不需要太多的计算资源来绘制人从地面上的阳台到达的距离作为初始速度
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的初始条件为你正在解决的牛顿公式是，
其中
m是该人的质量，g是我不知道英文名称的常数，但我们都知道它是什么 .来自eq . （3），
与eq . 一起使用（1），
我们也使用了初始条件 .
\int_0^t表示从0到t的积分 .来自eq . （4），
我们在哪里使用了初始条件 . 与eq . 一起使用（3）并且还使用初始条件，
其中
t^2表示t平方 . 从y(t)的表达式中，我们可以得到该人撞击地面的时间tstar. 也就是说，y(tstar) =0. 该等式可以通过二次公式（或类似的名称）来解决我在哪里使用了一个条件
tstar>0. 现在我们知道当他撞到地面时所达到的阳台的距离为x(tstar). 使用上面的x(t)表达式，实际上
x(tstar)取决于v0和theta以及g和ystar. 您将g和ystar作为常量，并且您希望找到所有(v0,theta)，使x(tstar) = xstar为给定的xstar值 .由于方程式的右侧（5）可以廉价计算，你可以在这个2D网格上为
v0和theta设置网格并计算xstar. 然后，您可以看到(v0,theta)谎言的解决方案集的大致位置 . 如果您需要精确的解决方案，您可以从该数据表中选取一个包含解决方案的段 .下面是一个演示这个想法的python脚本 .
我还附上了这个脚本生成的图形 . 黄色曲线是解决方案集
(v0,theta)，这样当你设置xstar = 8.0*1.95和ystar=4.0*1.95时，该人在墙上xstar处于地面 . 蓝色坐标表示x(tstar)，即人们从阳台水平跳跃的距离 . 请注意，在给定的v0（高于阈值aruondv0=9.9）时，有两个theta值（人物投射自己的两个方向）到达目标点(x,y) = (xstar,0).theta值的较小分支可以是负的，这意味着只要初始速度足够高，人就可以向下跳到达目标点 .该脚本还会生成一个数据文件
output.dat，其中包含解决方案封闭的段 .所以经过一些尝试后，我找到了一种方法来实现我想要的......这是我在首发帖子中提到的蛮力方法，但至少现在它有效......
这个想法很简单：定义一个函数
find_v0，它找到给定的thetaav0. 在这个函数中，你取一个v0的起始值（我选择8但这只是我的猜测），然后取起始值并检查difference函数与(xstar,0)有多远的距离 . 在这种情况下的有趣点可以通过将x轴上大于xstar的所有点设置为零（及其相应的y值），然后使用trim_zeros修剪所有零来确定，现在是对应的最后一个元素期望的输出 . 如果差函数的输出小于临界值（在我的情况下为0.1），则传递当前v0，如果不是，则将其增加0.01并再次执行相同的操作 .这个代码看起来像这样（再次替换3）和4））：
这件事的问题在于它需要永远计算（当前设置大约5-6分钟） . 如果有人有一些提示如何改进代码以获得更快或有不同的方法，它仍然会受到赞赏 .
假设x方向的速度永远不会降到零，你可以将x作为独立参数而不是时间 . 然后状态向量是时间，位置，速度，并且该状态空间中的向量场被缩放，使得vx分量始终为1.然后从零到xstar进行积分以计算轨迹符合xstar的状态（近似），如x-值 .
或者使用您自己的集成方法 . 我使用odeint作为显示的文档界面如何在集成中使用此衍生函数 .
得到的时间和y值可能是极端的