python - 实现数值解微分方程的初始条件

想象有人在一定角度和速度下从阳台上跳下（阳台的高度表示为theta）。从二维的角度看这个问题，考虑阻力，得到一个微分方程组，可以用Runge-Kutta方法来求解（我选择了显式中点，不确定这个问题的butcher表是什么）我实现了这个，它工作得很好，对于某些给定的初始条件，我得到了运动粒子的轨迹。
我的问题是，我要修复两个初始条件（x轴上的起点是零，y轴是v0），并确保轨迹在x轴上的某个点上（我们称之为ystar）。为此，当然存在两个初始条件的多重组合，在这种情况下，在x和y方向上是速度。问题是我不知道如何实现它。
我用来解决这个问题的代码：
1）龙格库塔法的实施

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）微分方程右端和Nessery参数的实现

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()

要使问题具体化：考虑到这一点，我如何找到ystar和xstar的所有可能组合，以便v0
当然，我尝试了一些东西，主要是蛮力法固定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()

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

最佳答案

编辑2018-07-16
在这里，我张贴了一个考虑到空气阻力的正确答案。
下面是一个python脚本，用于计算一组(v0,theta)值，以便空气拖动的轨迹在某个时间通过(x,y) = (xstar,0)。我使用无空气阻力的轨迹作为初始猜测，并猜测第一次优化时t=tstar对x(tstar)的依赖性。达到正确v0所需的迭代次数通常为3到4次脚本在我的笔记本电脑上0.99秒内完成，包括生成图形的时间。
脚本生成两个图形和一个文本文件。
v0
黑点表示解决方案集fig_xdrop_v0_theta.png
黄线表示参考值(v0,theta)，如果没有空气阻力，这将是一个解决方案。
(v0,theta)
从解决方案集中采样时，检查轨迹（蓝色实线）是否通过fig_traj_sample.png。
黑色虚线显示的轨迹没有被空气拖动作为参考。
(x,y)=(xstar,0)
包含(v0,theta)的数值数据、着陆时间output.dat和找到(v0,theta)所需的迭代次数。
剧本开始了。

#!/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()

这是数字tstar。

这是数字v0。

编辑2018-07-15
我意识到我忽略了这个问题考虑空气阻力我真丢人。所以，我下面的回答是不正确的。我担心自己删除我的答案看起来像是隐藏了一个错误，我暂时把它放在下面。如果有人觉得一个错误的答案让人恼火，我同意。有人删除它。
微分方程实际上可以用手求解，
而且不需要太多的计算资源
画出这个人从阳台上走多远
作为初始速度fig_xdrop_v0_theta.png和
角度fig_traj_sample.png然后，您可以选择条件v0
以至于theta
从这个数据表。
让我们写出
时间(v0,theta)的人分别是distance_from_balcony_on_the_ground(v0,theta) = xstar和t。
我想你在大楼的墙上
作为地面，我也是这样做的。假设
当时人的水平和垂直速度x(t)
分别是y(t)和x=0。
y=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)

其中t是人的质量，
v_x(t)是一个常数，我不知道它的英文名，
但我们都知道这是什么。
从式（3）中，

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

将其与式（1）结合使用，

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

我们也用了初始条件。v_y(t)表示
从t=0到m的积分。
从式（4）中，

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

我们使用初始条件的地方。
利用公式（3）和初始条件，

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

其中g表示\int_0^t平方。
从0的表达式中，我们可以得到时间t
在那人撞到地上。也就是说，t^2。
这个方程可以用二次公式来求解
（或类似的名字）如

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

我用了一个条件。现在我们知道了
他撞到阳台时所达到的距离
地面为t。使用上述y(t)的表达式，

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

实际上tstar取决于y(tstar) =0和tstar>0以及x(tstar)和x(t)。
您将x(tstar)和v0作为常量，并希望找到
对于给定的theta值，所有值都g。
由于公式（5）的右侧可以便宜地计算，
可以为ystar和g设置网格并计算ystar
在这个二维网格上。然后，您可以看到解决方案集的大致位置
的谎言。如果你需要精确的解决方案，你可以
将此数据表中的解决方案括起来的段。
下面是一个python脚本，演示了这个想法。
我还附上了一个由这个脚本生成的图形。
黄色曲线是解决方案集
有人在距离墙(v0,theta)处摔倒
当设置为x(tstar) = xstar和xstar时。
蓝色坐标表示v0，即
有人从阳台上横跳下来。
注意，在给定的theta（高于阈值aroundxstar），
有两个(v0,theta)值（人的两个方向
投射自己）到达目标点。
(v0,theta)值的较小分支可以是负的，这意味着只要初始速度足够高，人就可以向下跳到目标点。
该脚本还生成一个数据文件xstar，其中
封闭段的解决方案。

#!/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()

关于python - 实现数值解微分方程的初始条件，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/51333319/