使用odeint查找零值

Question

我该如何使用scipy.integrate.ode找到我的方程的第一个导数等于0的点呢？

我设置了这个函数，可以得到答案，但我不太确定结果的准确性，而且这可能不是最有效的方法。

基本上，我使用这个函数来找出一个初速度的抛体停止运动的时间。对于常微分方程的系统，有没有更好的方法来解决这个问题呢？

import numpy as np
from scipy import integrate

def find_nearest(array,value):
    idx=(np.abs(array-value)).argmin()
    return array[idx], idx

def deriv(x,t):
    # This function sets up the following relations
    # dx/dt = v , dv/dt = -(Cp/m)*(4+v^2)
    return np.array([ x[1], -(0.005/0.1) * (4+ (x[1]**2)) ])

def findzero(start, stop, v0):
    time = np.linspace(start, stop, 100000)
    #xinit are initial vaules of equation
    xinit = np.array([0.0,v0]) 
    x = integrate.odeint(deriv,xinit,time)
    # find nearest velocity value nearest to 0
    value, num = find_nearest(x[:,1],0.0001)
    print 'closest value ',
    print value
    print 'reaches zero at time ',
    print time[num]
    return time[num]
# from 0 to 20 seconds with initial velocity of 100 m/s
b = findzero(0.0,20.0,100.0)

Answer 1

我会使用类似于 scipy.optimize.fsolve() 这样的工具来找出导数的根。通过这个方法，我们可以反向推算出达到这个根所需的时间。

Answer 2

一般来说，解决这类问题的好方法是把你的方程式改写一下，让速度成为自变量，而时间和距离则是因变量。这样，你只需要从初始速度v0积分到速度为0就可以了。

不过，在你给出的例子中，其实根本不需要用到scipy.integrate。这些方程式用纸和笔就能轻松解决（通过变量分离再进行标准积分）。最后的结果是：

t = (m/(2 Cp)) arctan(v0/2)

这里的v0是初始速度，而arctan的结果需要用弧度来计算。

如果初始速度是100米每秒，那么答案是15.5079899282秒。